admas4.m

来自「用Admas四阶预估-校正算法的PECE模式求下列初值问题的数值解」· M 代码 · 共 31 行

function [EX,T,U,E]=adams4(a,b,c,h);
n=round((b-a)/h)+1;
t=zeros(n,1);
u=zeros(n,1);
e=zeros(n,1);
ex=zeros(n,1);
t(1,1)=a;
u(1,1)=c;
ex(1)=exat(0);
for i=1:3
    t(i+1)=t(i)+h;
    u(i+1)=u(i)+h*f(t(i),u(i));
    ex(i+1)=exat(t(i+1));
    e(i+1)=abs(u(i+1)-ex(i+1));
    str=sprintf('%d   t=%g   u=%g   ex=%g   e=%g ',i,t(i+1),u(i+1),ex(i+1),e(i+1));
    disp(str);
    i=i+1;
end
for i=4:n-1
    t(i+1)=t(i)+h;
    u(i+1)=u(i)+h/24*(55*f(t(i),u(i))-59*f(t(i-1),u(i-1))+37*f(t(i-2),u(i-2))-9*f(t(i-3),u(i-3)));
    ff=f(t(i+1),u(i+1));
    u(i+1)=u(i)+h/24*(9*ff+19*f(t(i),u(i))-5*f(t(i-1),u(i-1))+f(t(i-2),u(i-2)));
    ex(i+1)=exat(t(i+1));
    e(i+1)=abs(u(i+1)-ex(i+1));
    str=sprintf('%d   t=%g   u=%g   ex=%g   e=%g ',i,t(i+1),u(i+1),ex(i+1),e(i+1));
    disp(str);
      i=i+1;
end
T=t;U=u;E=e;EX=ex;