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