adamshh.m

来自「常微分方程组adams数值解法以及非线性方程的二分法解法」· M 代码 · 共 23 行

Adamshh.m
function [x,y]=Adamshh(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    if(n<4)
    k1=feval(fun,x(n),y(n));
    k2=feval(fun,x(n)+h/2,y(n)+k1*h/2);
    k3=feval(fun,x(n)+h/2,y(n)+k2*h/2);
    k4=feval(fun,x(n)+h,y(n)+k3*h);
    y(n+1)=y(n)+1/6*h*(k1+2*(k2+k3)+k4);
    else
       k1=feval(fun,x(n),y(n));
       k2=feval(fun,x(n-1),y(n-1));
       k3=feval(fun,x(n-2),y(n-2));
       k4=feval(fun,x(n-3),y(n-3));
       p4=y(n)+1/24*h*(k1*55-59*k2+37*k3-9*k4);
       dm4=feval(fun,x(n+1),p4);
       y(n+1)=y(n)+1/24*h*(9*dm4+k1*19-5*k2+k3);
    end
end