四阶龙格库塔法matlab实现

% 废话不多说,直接复制就可以运行,文章底部有参

考例子,不明白就看看

function [ yy ] = R_K_4( f , y0 , x0 , xn , hh )

% f 是inline function的句柄

% y0 是初始值

% [x0 xn]是计算范围

% hh 是步长

%% ===输入判断===

if (4 == nargin) % == 当没有给步长的时候使用默认步长==

h = 0.1 ;

elseif (5 == nargin)

h = hh ;

end

%% ===参数初始化(可以修改,可从网上找相关参数的值),以下是通用参数===

c = [1/6 1/3 1/3 1/6] ;

r = [0.5 0.5 1] ;

u = [0.5 0 0.5 0 0 1] ;

%% ===

x = x0:h:xn ; % == 计算节点==

y = y0 ;

N = length(x) ; % == 节点数==

for k = 1:N-1 % == 以下就直接按书本公式写的,应该容易理解== kk(1) = f( x(k) , y(k) ) ;

kk(2) = f( x(k)+r(1)*h , y(k) + u(1)*h*kk(1) ) ;

kk(3) = f( x(k)+r(2)*h , y(k) + u(2)*h*kk(1) + u(3)*h*kk(2) ) ;

kk(4) = f( x(k)+r(3)*h , y(k) + u(4)*h*kk(1) + u(5)*h*kk(2) + u(6)*h*kk(3) ) ;

y(k+1) = y(k) + h* ( c(1)*kk(1) + c(2)*kk(2) + c(3)*kk(3) + c(4)*kk(4) ) ;

end

%% ===输出处理===

if (0 == nargout)

for k = 1:N

fprintf( 'x: %f \\t y: %f \\n' , x(k) , y(k) ) ;

end

elseif (1 == nargout)

yy = y ;

else

disp('error : please check your input') ; return ; end end %% === 例子======= % % === input ========= % f = inline('x+y' , 'x' , 'y') ; % R_K_4(f , 1 , 0 , 1 , 0.2) % % === output ======== % x: 0.000000 y: 1.000000 % x: 0.200000 y: 1.242800 % x: 0.400000 y: 1.583636