rkf45v.asv

自己编写的变步长龙格库塔方法解强非线性微分方程

function [t,y]=rkf45v(tspan,y0,e,h)
if nargin<4
    h=(tspan(2)-tspan(1))/100;
end
hmin=h/64;
hmax=h*64;
n=1;
t(n)=tspan(1);
y(:,n)=y0;


while t(n)<tspan(2)
    if t(n)+h>br
        h=b-t(n);
    end
    br=tspan(2)-0.00001*abs(tspan(2));
    [y1,y2]=rkf45(t(n),y(:,n),h);
    err=abs(y2(2)-y1(2))+abs(y2(1)-y1(1));
    if ((err<e)|(h<2*hmin))
        y(n+1)=y1;
        
        if((t(n)+h>br)
            t(n+1)=tspan(2);
        else
            t(n+1)=t(n)+h;
        end
        
        n=n+1;
        
        if (err==0)
            s=0;
        else
            s=0.84*(e*h/err)^0.25;
        end
        
        if((s<0.75)&(h>2*hmin))
            h=h/2;
        end
        
        if((s>1.5)&(2*h<hmax))
            h=2*h;
        end
        
        if
        while (abs(y2(2)-y1(2))+abs(y2(1)-y1(1)))/10>e
            h=h/2;
            [y1,y2]=comput(t(n),y(:,n),h);
        end
    else
        while (abs(y2(2)-y1(2))+abs(y2(1)-y1(1)))/10<=e
            h=2*h;
            [y1,y2]=comput(t(n),y(:,n),h);
        end
        h=h/2;
        h=min(h,tspan(2)-t(n));
        [y1,y2]=comput(t(n),y(:,n),h);
    end
    n=n+1;
    t(n)=t(n-1)+h;
    y(:,n)=y2;
    [y1,y2]=comput(t(n),y(:,n),h);
end
t=t'; y=y';