多步法求解微分方程数值解
对于求解一个较大的系统来说,由于单步法所需要的信息较小导致精度不高,从而需要使用多步法来填补单步法的不足之处。常见的多步法有:Adams显示格式与Adams预测校正格式等。
Matalb代码:
1、Adams方法
function [x,y]=Adams4(f1,tspan,h,x0,y0) x=(tspan(1):h:tspan(2))';%等分区间 n=length(x); y=zeros(n,1);%存储y值 x(1)=x0; y(1)=y0; %使用Rungekutta4方法,为Adams4方法的启动赋初值。 for j=1:3 K1=feval(f1,x(j),y(j)); K2=feval(f1,x(j)+h/2,y(j)+h/2*K1); K3=feval(f1,x(j)+h/2,y(j)+h/2*K2); K4=feval(f1,x(j)+h,y(j)+h*K3); y(j+1)=y(j)+(h/6)*(K1+2*K2+2*K3+K4); end fn=feval(f1,x(1:4),y(1:4)); beta=[-9;37;-59;55]./24; for j=1:n-4 y(j+4)=y(j+3)+h.*(fn')*beta; fn=[fn(2:4);feval(f1,x(j+4),y(j+4))]; end
2、Adams校正格式
function [x,y]=Adams4_jiaozheng(f1,tspan,h,x0,y0) x=(tspan(1):h:tspan(2))';%等分区间 n=length(x); y=zeros(n,1);%存储y值 x(1)=x0; y(1)=y0; %使用Rungekutta4方法,为Adams4_jiaozheng方法的启动赋初值。 for j=1:3 K1=feval(f1,x(j),y(j)); K2=feval(f1,x(j)+h/2,y(j)+h/2*K1); K3=feval(f1,x(j)+h/2,y(j)+h/2*K2); K4=feval(f1,x(j)+h,y(j)+h*K3); y(j+1)=y(j)+(h/6)*(K1+2*K2+2*K3+K4); end fn=feval(f1,x(1:4),y(1:4)); beta=[-9;37;-59;55]/24; beta1=[1;-5;19;9]/24; %Adams4_jiaozheng方法 for j=1:n-4 y(j+4)=y(j+3)+h*fn'*beta; fn=[fn(2:4);feval(f1,x(j+4),y(j+4))]; y(j+4)=y(j+3)+h*fn'*beta1; fn=[fn(1:3);feval(f1,x(j+4),y(j+4))]; end