用于解决ODE系统。

function [tt,xx,k]=rkf45 (ftys,a,Alpha,b,h,delta,M)

%---------------------------------------------------------------------------
%RKF45   % Runge-Kutta-Fehlberg solution for a system Y' = F(T,Y) with Y(a) = Y_a.
% Sample call
%   [tt,xx,k] = rkf45('f',a,Alpha,b,h,delta,M)
% Inputs
%   f   name of the function
%   a      left  endpoint of [a,b]
%   Alpha  initial value
%   b      right endpoint of [a,b]
%   h      initial step size
%   delta  threshhold for stepsize acceptance 
%   M      maximum number of iterations allowed
% Outputs
%   tt     solution: array of points during the computation
%   xx     solution: array of corresponding approximations at the points.
%   k      number of iterations excuted
% Math5466: Numerical Analysis. Second project(c) Hong Zhang 2010
%
% Algorithm Runge-Kutta-Fehlberg Method RKF45
%---------------------------------------------------------------------------

ka=[16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55];
kamb=[1/360, 0, -128/4275, -2197/75240, 1/50, 2/55];
kc=[0, 1/4, 3/8, 12/13, 1, 1/2];
kd=[0, 0, 0, 0 ,0;1/4, 0, 0, 0, 0;3/32, 9/32, 0, 0, 0;1932/2197, -7200/2197, 7296/2197, 0, 0;...
    439/216, -8, 3680/513, -845/4104, 0;-8/27, 2, -3544/2565, 1859/4104, -11/40];
t=a;
X=Alpha;
iflag=1;
k=1;
tt (1)=t;
xx (1,:)=X;
while (k<M)
    d=b-t;
    if (abs (d)<= abs (h))
       iflag=0;
       h=d;
    end
    s=t;
    Y=X;
    Y1=feval (ftys, s+kc (1)*h, Y);
    Y2=feval (ftys, s+kc (2)*h, Y+kd (2,1)*h*Y1);
    Y3=feval (ftys, s+kc (3)*h, Y+(kd (3,1)*Y1+kd (3,2)*Y2)*h);
    Y4=feval (ftys, s+kc (4)*h, Y+(kd (4,1)*Y1+kd (4,2)*Y2+kd (4,3)*Y3)*h);
    Y5=feval (ftys, s+kc (5)*h, Y+(kd (5,1)*Y1+kd (5,2)*Y2+kd (5,3)*Y3+kd (5,4)*Y4)*h);
    Y6=feval (ftys, s+kc (6)*h, Y+(kd (6,1)*Y1+kd (6,2)*Y2+kd (6,3)*Y3+kd (6,4)*Y4+kd (6,5)*Y5)*h);
    X=Y+h*(ka (1)*Y1+ka (2)*Y2+ka (3)*Y3+ka (4)*Y4+ka (5)*Y5+ka (6)*Y6); 
    E=kamb (1)*Y1+kamb (2)*Y2+kamb (3)*Y3+kamb (4)*Y4+kamb (5)*Y5+kamb (6)*Y6;% error
    t=t+h;
    k=k+1;
    tt (k)=t;
    xx (k,:)=X;
    if (iflag==0)
       break;
    end

    if (max (abs (E))>=delta)% decline
       t=s;
       h=h/2;
       X=Y;
       k=k-1;
    elseif (max (abs (E))< delta/128)
       h=2*h;
    end   
end

转载本文请联系原作者获取授权，同时请注明本文来自张弘科学网博客。
链接地址：https://blog.sciencenet.cn/blog-398465-326487.html

上一篇：A Numerical Integration Problem
下一篇：A simple method "perl pie "to replace strings in multiple files on Linux system