解常微分方程的多步法（有admas外插法，admas内插法，一般多步法（Milne法））

最新推荐文章于 2022-12-22 11:10:47 发布

原创最新推荐文章于 2022-12-22 11:10:47 发布 · 5.4k 阅读

43 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #matlab

本文介绍了常微分方程数值解法中的多步法，包括Admas外插法、Admas内插法及Milne法，并通过具体实例展示了如何运用这些方法进行求解。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

解常微分方程的多步法（有admas外插法，admas内插法，一般多步法（Milne法））

admas外插法

%admas外插法，具有k-1阶的精度,下面编写k等于3的外插公式
clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
y1=0.001;
y2=0.008;
y3=0.027;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
yy=zeros(n,1);
yy(1)=y0;
yy(2)=y1;
yy(3)=y2;
yy(4)=y3;
for i=1:n
    xx(i)=a+h*(i-1); 
end
for i=5:n
   k1=55/24*h*(subs(f,{x,y},{xx(i-1),yy(i-1)}));
   k2=-59/24*h*(subs(f,{x,y},{xx(i-2),yy(i-2)})); 
   k3=37/24*h*(subs(f,{x,y},{xx(i-3),yy(i-3)}));
   k4=-9/24*h*(subs(f,{x,y},{xx(i-4),yy(i-4)}));
   yy(i)=k1+k2+k3+k4+yy(i-1);
end
ff1=[xx,yy];

admas内插法

%admas内插法，具有k阶的精度,但是为隐式方程，所以需要求解方程，下面编写当p等于1，k=4的admas内插法
clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
y1=0.001;
y2=0.008;
y3=0.027;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
yy=zeros(n,1);
yy(1)=y0;
yy(2)=y1;
yy(3)=y2;
yy(4)=y3;
for i=1:n
    xx(i)=a+h*(i-1); 
end
for i=5:n
   syms p
   k1=subs(f,{x,y},{xx(i),p})*h*251;
   k2=subs(f,{x,y},{xx(i-1),yy(i-1)})*h*646;
   k3=subs(f,{x,y},{xx(i-2),yy(i-2)})*h*(-264);
   k4=subs(f,{x,y},{xx(i-3),yy(i-3)})*h*106;
   k5=subs(f,{x,y},{xx(i-4),yy(i-4)})*h*(-19);
   yy(i)=double(solve(p==yy(i-1)+k1+k2+k3+k4+k5,p));
end
ff1=[xx,yy];

一般多步法（Milne法）

%一般线性多步法，下面介绍Milne法
clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
y1=0.001;
y2=0.008;
y3=0.027;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
yy=zeros(n,1);
yy(1)=y0;
yy(2)=y1;
for i=1:n
    xx(i)=a+h*(i-1); 
end
for i=3:n
   syms p
   k1=h/3*(subs(f,{x,y},{xx(i),p})+4*subs(f,{x,y},{xx(i-1),yy(i-1)})+subs(f,{x,y},{xx(i-2),yy(i-2)}));
   yy(i)=double(solve(p==yy(i-2)+k1,p));
end
ff1=[xx,yy];