rkf45v.m

来自「自己编写的变步长龙格库塔方法解强非线性微分方程」· M 代码 · 共 46 行

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;

br=tspan(2)-0.00001*abs(tspan(2));
while t(n)<tspan(2)
    if t(n)+h>br
        h=b-t(n);
    end
    
    [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;
    end
        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
end
t=t'; y=y';