MATLAB龙格-库塔方法解微分方程

龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例，计算格式如下

()()()11234121324

3226,,22,+22,n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎛⎫=++⎨ ⎪⎝⎭⎪⎪⎛⎫=+⎪ ⎪⎝⎭⎪⎪=++⎩

1龙格-库塔法解一阶ODE 对于形如()()0, dy f x y a x b dx y a y ⎧=<≤⎪⎨⎪=⎩

的一阶ODE 初值问题，可以直接套用公式，如今可以借助计算机方便的进行计算，下面给出一个实例

()2 0101dy x y x dx y y ⎧=-<≤⎪⎨⎪=⎩

取步长h=0.1

，此处由数学知识可得该方程的精确解为y =。在这里利用MATLAB 编程，计算数值解并与精确解相比，代码如下：

(1)写出微分方程，便于调用和修改

function val = odefun( x,y )

val = y-2*x/y;

end

(2)编写runge-kutta 方法的函数代码

function y = runge_kutta( h,x0,y0 )

k1 = odefun(x0,y0);

k2 = odefun(x0+h/2,y0+h/2*k1);

k3 = odefun(x0+h/2,y0+h/2*k2);

k4 = odefun(x0+h,y0+h*k3);

y = y0+h*(k1+2*k2+2*k3+k4)/6;

end

(3)编写主函数解微分方程,并观察数值解与精确解的差异clear all

h = 0.1;

x0 = 0;

y0 = 1;

x = 0.1:h:1;

y(1) = runge_kutta(h,x0,y0);

for k=1:length(x)

x(k) = x0+k*h;

y(k+1) = runge_kutta(h,x(k),y(k));

end

z = sqrt(1+2*x);

plot(x,y,’*’);

hold on

plot(x,z,'r');

结果如下图，数值解与解析解高度一致

2龙格-库塔法解高阶ODE

对于高阶ODE来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处仍以一个实例进行说明。

++=

5002007502000

y y y

这是一个二阶ODE,描述的是一个物体的有阻尼振动情况。

初始条件为000 0y y == ，将方程降阶，引入一个向量型变量Y

y Y y ⎡⎤=⎢⎥⎣⎦故有()2000200750500y y dY Y y y y dt ⎡⎤⎡⎤⎢⎥===--⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦

记()()1 2y Y y Y ==则()()()()2200020027501500Y Y Y Y ⎡⎤⎢⎥=--⎢⎥⎢⎥⎣⎦

至此，二阶方程降阶为一阶

方程组。值得注意的是此时再用龙格-库塔法进行求解时，代入的将是一个Y 向量。同样利用MATLAB 进行计算，步长h=0.05，时间周期为[0,20].

(1) 编写ODE 函数

function Y = odefun1( ~,Y0 )

% 此处Y0为一个列向量，因为时间t 未显含在一阶方程组中 % 所以ode 函数的第一个参数为空，要根据具体情况而定。 Y = [Y0(2);

(2000-200*Y0(2)-750*Y0(1))/500;];

end

(2) 编写runge-kutta 函数

function Y = rkfa( h,t0,Y0 )

k1 = odefun1(t0,Y0);

k2 = odefun1(t0+h/2,Y0+h/2*k1);

k3 = odefun1(t0+h/2,Y0+h/2*k2);

k4 = odefun1(t0+h,Y0+h*k3);

Y = Y0+h*(k1+2*k2+2*k3+k4)/6;

end

(3)编写主函数

clear all

h = 0.05;

t = 0.05:h:20;

t0 = 0;

Y0 = [0;

0];%初值

Y = cell(1,length(t));

Y{1} = rkfa( h,t0,Y0 );

z = zeros(2,length(t));

for k=1:length(t)

Y{k+1} = rkfa( h,t0,Y{k}); z(1,k) = Y{k}(1);

z(2,k) = Y{k}(2);

end

plot(t,z(1,:),'r');%位移y的图像hold on

plot(t,z(2,:));%速度y’的图像

求解结果如下图