子函数,四阶经典的龙格-库塔函数:
function [x,y]=Classical_RK4(odefun,xspan,y0,h,varargin)
% 经典Runge-Kutta法求解常微分方程
% 输入参数:
% ---odefun:微分方程的函数描述
% ---xspan:求解区间[x0,xn]
% ---y0:初始条件
% ---h:迭代步长
% ---p1,p2,…:odefun函数的附加参数
% 输出参数:
% ---x:返回的节点,即x=xspan(1):h:xspan(2)
% ---y:微分方程的数值解
x=xspan(1):h:xspan(2);
y(1)=y0;
for k=1:length(x)-1
K1=feval(odefun,x(k),y(k),varargin{:});
K2=feval(odefun,x(k)+h/2,y(k)+h/2*K1,varargin{:});
K3=feval(odefun,x(k)+h/2,y(k)+h/2*K2,varargin{:});
K4=feval(odefun,x(k)+h,y(k)+h*K3,varargin{:});
y(k+1)=y(k)+h/6*(K1+2*K2+2*K3+K4);
end
x=x';y=y';
2.主函数:
>> f=@(t,y)8-3*y;
>> xspan=[0 4];
>> h=0.2;
>> y0=2;
>> [t,y]=Classical_RK4(f,xspan,y0,h)
t =
0
0.2000
0.4000
0.6000
0.8000
1.0000
1.2000
1.4000
1.6000
1.8000
2.0000
2.2000
2.4000
2.6000
2.8000
3.0000
3.2000
3.4000
3.6000
3.8000
4.0000
y =
2.0000
2.3004
2.4654
2.5561
2.6059
2.6333
2.6483
2.6566
2.6611
2.6636
2.6650
2.6657
2.6662
2.6664
2.6665
2.6666
2.6666
2.6666
2.6667
2.6667
2.6667