龙格库塔法RKF45Matlab实现

龙格库塔法RKF45的Matlab实现

2007-08-16 14:03:32| 分类：MatLab/Maple/Mat|字号订阅

4阶5级龙格库塔法用于解一阶微分方程（组），对于高阶微分方程，可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写（数值方法matlab版），但只能解微分方程，不能解微分方程组。由LiuLiu@uestc修改，使之能够解微分方程组。该程序精度比matlab自带的ode45更高。

rkf45.m:

function [Rt Rx]=rkf45(f,tspan,ya,m,tol)

% Input:

% - f function column vector

% - tspan[a,b] left & right point of [a,b]

% - ya initial value column vector

% -m initial guess for number of steps

% -tol tolerance

% Output:

% - Rt solution: vector of abscissas

% - Rx solution: vector of ordinates

% Program by John.Mathews, improved by liuliu@uestc

if length(tspan)~=2

error('length of vector tspan must be 2.');

end

if ~isnumeric(tspan)

error('TSPAN should be a vector of integration steps.');

end

if ~isnumeric(ya)

error('Ya should be a vector of initial conditions.');

end

h = diff(tspan);

if any(sign(h(1))*h <= 0)

error('Entries of TSPAN are not in order.') ;

end

a=tspan(1);

b=tspan(2);

ya=ya(:);

a2 = 1/4; b2 = 1/4; a3 = 3/8; b3 = 3/32; c3 = 9/32; a4 = 12/13;

b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197; a5 = 1;

b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104; a6 = 1/2;

b6 = -8/27; c6 = 2; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40;

r1 = 1/360; r3 = -128/4275; r4 = -2197/75240; r5 = 1/50;

r6 = 2/55; n1 = 25/216; n3 = 1408/2565; n4 = 2197/4104; n5 = -1/5;

big = 1e15;

h = (b-a)/m;

hmin = h/64;% 步长自适应范围下限

hmax = 64*h;% 步长自适应范围上限

max1 = 200;% 迭代次数上限

Y(1,:) = ya;

T(1) = a;

j = 1;

% tj = T(1);

br = b - 0.00001*abs(b);

while (T(j)

if ((T(j)+h)>br), h = b - T(j); end

%caculate values of k1...k6,y1...y6

tj = T(j);

yj = Y(j,:);

y1 = yj;

k1 = h*feval(f,tj,y1);

y2 = yj+b2*k1;

if big

k2 = h*feval(f,tj+a2*h,y2);

y3 = yj+b3*k1+c3*k2; if big

k3 = h*feval(f,tj+a3*h,y3);

y4 = yj+b4*k1+c4*k2+d4*k3; if big

k4 = h*feval(f,tj+a4*h,y4);

y5 = yj+b5*k1+c5.*k2+d5*k3+e5*k4; if big

y6 = yj+b6*k1+c6.*k2+d6*k3+e6*k4+f6*k5; if big

k6 = h*feval(f,tj+a6*h,y6);

err = abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);

ynew = yj+n1*k1+n3*k3+n4*k4+n5*k5;

% error and step size control

if ( (err

Y(j+1,:) = ynew;

if ((tj+h)>br),

T(j+1) = b;

else

T(j+1) = tj + h;

end

j = j+1;

tj = T(j);

end

if (max(err)==0),

s = 0;

else

s1 = 0.84*(tol.*h./err).^(0.25);% 最佳步长值

s=min(s1);

end

if ((s<0.75)&(h>2*hmin)), h = h/2; end

if ((s>1.50)&(2*h

if ( (big

end

% [Rt Rx]=[T' Y];

Rt=T';

Rx=Y;

使用方法：

首先编写方程（组）文件（注意与ode45不同，这儿方程组为1Xn数组：

function dx= fun(t,x)

dx=zeros(1,2);

dx(1)=x(1)+x(2)*2+0*t;

dx(2)=3*x(1)+x(2)*2+0*t;

然后使用：

[Rt,Rx]=rkf45(@fun,[0,0.2],[6;4],100,1e-7)