Adams预估计-校正法解微分方程。
function T=Adams_PreEs(f,ab,y0,h)
x=ab(1):h:ab(2);
n=length(x);
T=zeros(2,n);
T(2,:)=x;
T(1,1)=y0;
for t=2:4
x0=T(2,t-1);
k1=feval(f,x0,T(1,t-1));
k2=feval(f,x0+h/2,k1*(h/2)+T(1,t-1));
k3=feval(f,x0+h/2,k2*(h/2)+T(1,t-1));
k4=feval(f,x0+h,k3*h+T(1,t-1));
T(1,t)=T(1,t-1)+(k1+k2*2+k3*2+k4)*(h/6);
end
fn=zeros(1,4);
for t=5:n
x0=T(2,t);
for d=1:4
fn(d)=feval(f,x0-d*h,T(1,t-d));
end
yp=T(1,t-1)+(h*fn*[55,-59,37,-9]')/24;
fn=[feval(f,x0,yp),fn(1:3)];
T(1,t)=T(1,t-1)+(h*fn*[9,19,-5,1]')/24;
end
例: 求解初值问题:
解:(1)建立函数文件:
function dy=dfun(x,y)
dy=3*exp(-x)*cos(y)-x*y;
(2)在command window中调用:
>> T=Adams_PreEs('dfun',[0 3],0,0.02);
>> plot(T(2,:),T(1,:))
(3)与Runge-Kutta方法的比较:
评论