多步法求解微分方程数值解

Posted 2020-07-27 视野

       对于求解一个较大的系统来说，由于单步法所需要的信息较小导致精度不高，从而需要使用多步法来填补单步法的不足之处。常见的多步法有：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