通用龙格库塔Runge-Kutta方法求解…

通用龙格库塔Runge-Kutta方法求解常微分方程组初值问题的C++优雅实现

1. 算法简介

a. 事情的起因

前一段时间在C++项目过程中，需要求解一个微分方程组，看了相关的数值分析教程（《数值分析》，欧阳洁等编著，北京：高等教育出版社，2009.9），又看了别人设计好的龙格库塔程序，觉得写得比较繁琐，而且不够通用。索性自己编写一个，借鉴了C++标准库中泛型函数的写法，设计了一个比较通用的龙格库塔算法。

b. 初值问题

常微分方程组初值问题是：

    dy/dx = f(x,y), a <= x <= b.

    y(a) = y0.

其中f(x,y)为x,y的已知实值函数，y0为给定的初始值。这里，在控制系统中，x常用时间t表示，x是标量；若y和f(x,y)是标量，则上述初值问题是常微分方程的初值问题；若y和f(x,y)是向量，则上述初值问题是常微分方程组的初值问题.可以将常微分方程的初值问题看作是常微分方程组的初值问题的特例。初值问题更详细的介绍请参考维基百科。

c. 龙格库塔算法（方法）

是龙格库塔算法（方法）求解常微分方程组初值问题的优秀方法，具有很高的求解精度。比如，四级四阶经典Runge-Kutta公式的局部截断误差是O(h5),即步长的5次方。可以根据Runge-Kutta方法的思路推导出更精确的数值算法。这里实现的Runge-Kutta算法特指四级四阶经典Runge-Kutta公式。四级四阶经典Runge-Kutta公式表示如下：

    y(n+1) = y(n) + h * (K1 + 2*K2 + 2*K3 + K4) / 6,

    K1 = f(x(n), y(n)),

    K2 = f(x(n) + h/2, y(n) + h*K1/2),

    K3 = f(x(n) + h/2, y(n) + h*K2/2),

    K4 = f(x(n) + h, y(n) + h*K3)。

式子中的h表示步长。

2. C++实现

a. <valarray>标准库简介

<valarray>是C++ STL（标准模板库）的一个组件，可以理解为用于存储要进行数值运算的数据的数组（value array，数值数组）。该模板类重载了<math>或者<math.h>中对数据进行运算的函数和四则运算等操作符。对valarray<type>对象的某个操作（运算）对施加到对象内的每个元素。因此，valarray<type>对象相当于一个向量。熟悉Matlab的读者会很好理解该对象，因为Matlab中的函数和运算符都可以用在向量上。比如，

std::valarray<double> x, y;

y = sin(x);

这种表达在Matlab中是很常见的，对于x向量的每个元素，求其正弦值，放在y中。使用了<valarray>之后，在C++中，也可以获得这种形式上的简介，省去了书写和阅读循环体的麻烦。

注意：在VC6.0版本中（其他版本VC未测试），引入<valarray>库会导致错误，因为在<windef.h>中定义了min和max宏，与<valarray>中的min和max函数有冲突。最简单的解决方法是：在"stdafx.h"中，增加

#define MOMINMAX  ，将该语句放在<stdafx.h>文件中所有的#include<*.h>之前。这种解决方法的原因可以从<windef.h>中看出来。以下直接拷贝<windef.h>中的相关部分。

#ifndef NOMINMAX

#ifndef max
#define max(a,b)            (((a) > (b)) ? (a) : (b))
#endif

#ifndef min
#define min(a,b)            (((a) < (b)) ? (a) : (b))
#endif

#endif  / * NOMINMAX * /

所以定义了NOMINMAX后就不会再定义min和max宏了。

b. 函数对象

函数对象是对函数的抽象，功能上类似函数，但比函数功能更强大。C++ STL的定义和实现使用了函数对象，使用C++ STL中也会大量遇到函数对象。这里使用函数对象来给Runge-Kutta算法函数传递微分方程组的右端项。读者可以参考C++ Primer了解更多关于函数对象的概念。

c. 单步计算

这里给出的Runge-Kutta算法函数只进行单步计算，即给出初值和步长，计算下一步的函数值。若想计算多步，需要使用循环体，在循环体中调用单步Runge-Kutta算法函数。

d. 额外参数

有时候，微分方程组右端项会含有一些与自变量x和因变量y无关的参数，所以这里在实现Runge-Kutta算法函数时给出了重载版本，可以传递额外的参数。

3.附录

附录分为三部分。在第一部分给出龙格库塔算法的完整源代码，在第二部分给出了该算法函数的使用示例，在第三部分给出了该示例程序调试输出的结果。给出的示例基于MFC的对话框程序。在对话框界面上添加Button，设置其ID为IDC_BU