matlab龙格库塔法调用,MATLAB龙格-库塔方法解微分方程.pdf

龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了

高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例，计算格式如下

 h

y y  K  2K  2K K 

n1 n 1 2 3 4

 6

K f x ,y

1  n n 



  h h 

K f x  ,y  K

2  n n 1 

  2 2 

  h h 

K f x  ,y + K

 3  n n 2 

  2 2 

K f x h,y hK

 4  n n 3 

1 - ODE

龙格库塔法解一阶

dy

f x,y ax b

  

对于形如dx 的一阶ODE 初值问题，可以直接套用公式，如今可以



y a y

   0

借助计算机方便的进行计算，下面给出一个实例

dy 2x

 y  0  x  1

dx y



y 0 1

  

取步长h=0.1，此处由数学知识可得该方程的精确解为y 12x 。在这里利用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 来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处

仍以一个实例进行说明。

 

500y 200y 750y 2000

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



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

0 0



 y 

y y

Y   故有Y dY   2000 200y 750y  

    