在决策理论与方法最后一课中，老师讲了连续系统建模与仿真，事实上在本科就已经接触过许多，在这里我把数值仿真方法中的龙格库塔法给详细介绍一下，这是个比较有趣，同时也比较实用的方法。在求解常微分方程初值问题的数值解经常会用到。本文主要参考书为马东升等主编的 《数值计算方法》，机械工业出版社出版。

初值问题

等价于

龙格库塔法的基本思路是用 $f(x,y)$ 在几个不同点的数值加权平均代替 $f(x_n+theta h,y(x_n+ theta h))$ 的值，而使截断误差的阶数尽可能高。也就是说，取不同点的斜率加权平均作为平均斜率，从而提高方法的阶数。这样龙格库塔法保留了泰勒展开法所具有的高阶局部截断误差，同时避免了计算函数 $f(x,y)$ 的高阶函数。

其中，$phi (x_n,y_n,h)$ 是 $f(x,y)$ 在一些点的线性组合，解释为 $f(x,y)$ 的平均斜率。上式中计算了 $r$ 个 $f$ 的函数值，称该式为 $r$ 级的龙格库塔法。这里系数 $c_i,lambda_i,mu_{ij}$ 均为常数，其选取原则是使 $y_{n+1}$ 得展开表达式

与微分方程的解 $y(x_{n+1})$ 在 $(x_n,y_n)$ 处的泰勒展开式

有尽可能多的项相重合，以减少局部截断误差。

经典龙格库塔法

二阶龙格库塔法

已知

适当选取 $alpha,beta,w_1,w_2​$ 的值，使局部截断误差 $y(x_{n+1})-y_{n+1}​$ 的阶数尽可能高。这里仍假定 $y_n=y(x_n)​$，显然 $k_1=y’(x_n)​$ 。

由二元函数的泰勒展开

将 $k_1,k_2$ 的表达式代入

为使 $y(x_{n+1})-y_{n+1}=O(h^3)$

比较 $y(x_{n+1})$ 和 $y_{n+1}$ ，可得

这是四个未知数三个方程的不定方程组。任一未知数可设为自由变量，求出其余三个未知数，这样每一组未知数的组合就确定了一种二阶龙格库塔式。

当取 $w_1=frac{1}{4}$ 时，$w_2=frac{3}{4},alpha=beta=frac{2}{3}$ ，则有

这种方法称为休恩公式。

当取 $w_1=frac{1}{2}$ 时，$w_2=frac{1}{2},alpha=beta=1$ ，此即为改进欧拉公式。

当取 $w_1=0$ 时，$w_2=1,alpha=beta=frac{1}{2}$ ，此时称为变形欧拉公式，

三阶和四阶龙格库塔法

推导过程如下，我们同理可得常用的三阶龙格库塔式为：

常用的四阶龙格库塔公式：

经典的龙格库塔法每一步需要4次计算函数值 $f(x,y)$ ，它具有四阶精度，即局部截断误差是 $O(h^5)$ .

许多计算实例表明，要达到相同的精度，经典公式的步长可以比二阶方法的步长大十倍，而经典公式每步的计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小。正是由于上述原因，工程上常用四阶经典龙格库塔法，高于四阶的方法由于每步计算量增加较多，而精度提高不快，因此使用得比较少。

举个例子

某地区某病菌传染的系统动力学模型为

式中，$x_1$ 表示可能受到传染的人数，$x_2$ 表示已经被传染的病人数，$x_3$ 表示已治愈的人数。使用matlab进行编程，对其进行仿真研究，并绘制出对应的时间相应曲线。

解：直接给出matlab代码1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18clear

x0=[620,10,70]; # 设置状态变量初值

tspan=[0,30]; # 设置仿真时间区间

[t,x]=ode45('funXX',tspan,x0); #调用ode45求仿真解

# 用不同的线型绘制仿真结果曲线

plot(t,x(:,1),t,x(:,2),‘--’,t,x(:,3),‘:’);

xlabel('time(天) t0=0, tf=30'); # 对t-x轴进行标注

ylabel('x(人):x1(0)=620,x2(0)=10,x3(0)=70');

legend('x1','x2','x3');

grid;

# funXX 保存为函数m文件

function xdot=funXX(t,x)

xdot1(1)=-0.001*x(1)*x(2); # 第一个微分方程

xdot1(2)=0.001*x(1)*x(2)-0.072*x(2); # 第二个微分方程

xdot1(3)=0.072*x(2);# 第三个微分方程

xdot=xdot1';

结果作图如下：

参考文献