|
通用龙格库塔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_BUTTON_TEST_ODE，增加消息相应函数void
CODEDlg::OnButtonTestOde()。在OnButtonTestOde()中给出了四种不同情况的使用示例，分别是：单变量常微分方程（右端项使用函数的形式传递），单变量常微分方程（右端项使用函数对象的形式传递），多变量微分方程组（普通版本），多变量微分方程组（重载版本，包含额外参数）。
a. 龙格库塔算法的完整源代码
// RungeKutta.h: interface for the
RungeKutta method.
//
//////////////////////////////////////////////////////////////////////
#ifndef RUNGE_KUTTA_H_H
#define RUNGE_KUTTA_H_H
#include <valarray>
//单步四级四阶经典龙格库塔算法
template <typename Func>
void ODERungeKuttaOneStep(double
dxn,
//x初值
const
std::valarray<double>&
dyn, //初值y(n)
double
dh,
//步长
Func
func,
//微分方程组右端项
std::valarray<double>&
dynext
//下一步的值y(n+1)
);
//单步四级四阶经典龙格库塔算法，重载版本, 含有额外参数
template <typename Func>
void ODERungeKuttaOneStep(double
dxn,
//x初值
const
std::valarray<double>&
dyn, //初值y(n)
double
dh,
//步长
Func
func,
//微分方程组右端项
std::valarray<double>&
dynext,
//下一步的值y(n+1)
const
std::valarray<double>&
para //额外参数
);
#include "RungeKutta.inl"
//template function implementation
#endif
// RungeKutta.inl: implementation of the
RungeKutta method.
//
//////////////////////////////////////////////////////////////////////
//单步四级四阶经典龙格库塔算法
// 功能：求解常微分方程组初值问题的四级四阶经典龙格库塔算法，向前计算一个步长
// 输入：输入参数有：dxn, x的初值; dyn, 初值向量y(n); dh, 步长;
//
func, 微分方程组右端项
// 输入输出参数有: dynext,
最好初始化为与dyn长度相等的序列
// 输出：无
// 作者：张坤
// 时间：2013.06.13
template <typename Func>
void ODERungeKuttaOneStep(double
dxn,
//x初值
const
std::valarray<double>&
dyn, //初值y(n)
double
dh,
//步长
Func
func,
//微分方程组右端项
std::valarray<double>&
dynext
//下一步的值y(n+1)
)
{
size_t n = dyn.size();
//微分方程组中方程的个数，同时是初值y(n)和下一步值y(n+1)的长度
if (n != dynext.size())
{
dynext.resize(n,
0.0); //下一步的值y(n+1)与y(n)长度相等
}
std::valarray<double>
K1(0.0, n), K2(0.0, n), K3(0.0, n), K4(0.0, n);
func(dxn, dyn,
K1);
//求解K1
func(dxn+dh/2, dyn+dh/2*K1,
K2); //求解K2
func(dxn+dh/2, dyn+dh/2*K2,
K3); //求解K3
func(dxn+dh, dyn+dh*K3,
K4);
//求解K4
dynext = dyn + (K1 + K2 + K3 +
K4)*dh/6.0; //求解下一步的值y(n+1)
}
//单步四级四阶经典龙格库塔算法，重载版本, 含有额外参数
// 功能：求解常微分方程组初值问题的四级四阶经典龙格库塔算法，向前计算一个步长
// 输入：输入参数有：dxn, x的初值; dyn, 初值向量y(n); dh, 步长;
//
func, 微分方程组右端项; para, 微分方程组可能需要的额外参数
// 输入输出参数有: dynext,
最好初始化为与dyn长度相等的序列
// 输出：无
// 作者：张坤
// 时间：2013.06.13
template <typename Func>
void ODERungeKuttaOneStep(double
dxn,
//x初值
const
std::valarray<double>&
dyn, //初值y(n)
double
dh,
//步长
Func
func,
//微分方程组右端项
std::valarray<double>&
dynext,
//下一步的值y(n+1)
const
std::valarray<double>&
para //额外参数
)
{
size_t n = dyn.size();
//微分方程组中方程的个数，同时是初值y(n)和下一步值y(n+1)的长度
if (n != dynext.size())
{
dynext.resize(n,
0.0); //下一步的值y(n+1)与y(n)长度相等
}
std::valarray<double>
K1(0.0, n), K2(0.0, n), K3(0.0, n), K4(0.0, n);
func(dxn, dyn, K1,
para);
//求解K1
func(dxn+dh/2, dyn+dh/2*K1, K2,
para); //求解K2
func(dxn+dh/2, dyn+dh/2*K2, K3,
para); //求解K3
func(dxn+dh, dyn+dh*K3, K4,
para);
//求解K4
dynext = dyn + (K1 + K2 + K3 + K4)*dh/6.0;
//求解下一步的值y(n+1)
}
b.龙格库塔算法的使用示例
void odefunc1(double dx, const
std::valarray<double>&
dyn,
std::valarray<double>&
fai)
{
fai[0] = dyn[0] - 2*dx/dyn[0];
// dy/dx = y - 2*x/y; 单变量微分方程
}
struct odefunc2
{
void operator()(double dx, const
std::valarray<double>&
dyn,
std::valarray<double>&
fai)
{
fai[0] = dyn[0] -
2*dx/dyn[0]; // dy/dx = y - 2*x/y; 单变量微分方程
}
}; //struct 后面可不要丢掉分号“;”
struct odefunc3
{
void operator()(double dx, const
std::valarray<double>&
dyn,
std::valarray<double>&
fai)
{
fai[0] = dyn[1] *
dyn[2]; // dy1/dx = y2 * y3; 3变量微分方程组
fai[1] = -dyn[0] * dyn[2]; //
dy2/dx = -y1 * y3;
fai[2] = -0.51 * dyn[0] *
dyn[1]; // dy3/dx = -0.51 * y1 * y2;
}
}; //struct 后面可不要丢掉分号“;”
struct odefunc4
{
void operator()(double dx, const
std::valarray<double>&
dyn,
std::valarray<double>&
fai, const
std::valarray<double>&
para)
{
fai[0] = para[0] * (para[3] -
dyn[0]); // dy1/dx = c1 * (d1 - y1);
3变量微分方程组，带额外参数
fai[1] = para[1] * (para[4] -
dyn[1]); // dy2/dx = c2 * (d2 - y2);
fai[2] = para[2] * (para[5] -
dyn[2]); // dy3/dx = c3 * (d3 - y3);
}
}; //struct 后面可不要丢掉分号“;”
void odeTest1() //单变量微分方程,使用函数指针的方式传递函数
{
double dh = 0.1; //步长
double dx0 = 0.0;
double dy0 = 1.0; //初值
std::valarray<double>
darrayn(dy0, 1); //初值，向量长度是1
std::valarray<double>
darraynext(0.0, 1); //下一步的值,最好初始化
double dx = dx0; //当前的x
TRACE("\nodeTest1\n"); //调试输出
TRACE("xn
yn\n");
TRACE("%.1f %.7f\n", dx,
darrayn[0]); //调试输出初值
for (int i = 0; i < 5;
++i) //向前计算5步
{
ODERungeKuttaOneStep(dx,
darrayn, dh, odefunc1, darraynext);
dx += dh; //更新x
darrayn =
darraynext;
//更新y(n)，为循环的下一步做准备
TRACE("%.1f
%.7f\n", dx, darrayn[0]); //调试输出计算一步后的值
}
}
void odeTest2() //单变量微分方程,使用函数算子的方式传递函数
{
double dh = 0.1; //步长
double dx0 = 0.0;
double dy0 = 1.0; //初值
std::valarray<double>
darrayn(dy0, 1); //初值，向量长度是1
std::valarray<double>
darraynext(0.0, 1); //下一步的值,最好初始化
double dx = dx0; //当前的x
TRACE("\nodeTest2\n"); //调试输出
TRACE("xn
yn\n");
TRACE("%.1f %.7f\n", dx,
darrayn[0]); //调试输出初值
for (int i = 0; i < 5;
++i) //向前计算5步
{
ODERungeKuttaOneStep(dx,
darrayn, dh, odefunc2(), darraynext);
dx += dh; //更新x
darrayn =
darraynext;
//更新y(n)，为循环的下一步做准备
TRACE("%.1f
%.7f\n", dx, darrayn[0]); //调试输出计算一步后的值
}
}
void odeTest3() //3变量微分方程组, 使用函数算子的方式传递函数
{
double dh = 0.01; //步长
double dx0 = 0.0;
double dy0[] = {0.0, 1.0, 1.0};
//初值
std::valarray<double>
darrayn(dy0, 3); //初值，向量长度是3
std::valarray<double>
darraynext(0.0, 3); //下一步的值,最好初始化
double dx = dx0; //当前的x
TRACE("\nodeTest3\n"); //调试输出
TRACE("xn
y1n
y2n
y3n\n");
//调试输出初值
TRACE("%.2f
%.7f %.7f %.7f\n", dx,
darrayn[0], darrayn[1], darrayn[2]);
for (int i = 0; i < 5;
++i) //向前计算5步
{
ODERungeKuttaOneStep(dx,
darrayn, dh, odefunc3(), darraynext);
dx += dh; //更新x
darrayn =
darraynext;
//更新y(n)，为循环的下一步做准备
//调试输出计算一步后的值
TRACE("%.2f
%.7f %.7f %.7f\n", dx,
darrayn[0], darrayn[1], darrayn[2]);
}
}
void odeTest4() //3变量微分方程组,
带额外参数，使用函数算子的方式传递函数，其意义是求3个通道的阶跃响应
{
double dh = 0.1; //步长
double dx0 = 0.0;
double dy0[] = {0.0, 0.0, 0.0};
//初值
std::valarray<double>
darrayn(dy0, 3); //初值，向量长度是3
std::valarray<double>
darraynext(0.0, 3); //下一步的值,最好初始化
double dcd[] = {1.0, 1.0, 1.0, 1.0, 1.0,
1.0}; //前3个量是惯性环节的开环增益, 后三个量是三个通道的单位阶跃函数
std::valarray<double>
para(dcd, 6);
double dx = dx0; //当前的x
TRACE("\nodeTest3\n"); //调试输出
TRACE("xn
y1n
y2n
y3n\n");
//调试输出初值
TRACE("%.1f
%.7f %.7f %.7f\n", dx,
darrayn[0], darrayn[1], darrayn[2]);
for (int i = 0; i < 5;
++i) //向前计算5步
{
ODERungeKuttaOneStep(dx,
darrayn, dh, odefunc4(), darraynext, para);
dx += dh; //更新x
darrayn =
darraynext;
//更新y(n)，为循环的下一步做准备
//调试输出计算一步后的值
TRACE("%.1f
%.7f %.7f %.7f\n", dx,
darrayn[0], darrayn[1], darrayn[2]);
}
}
void CODEDlg::OnButtonTestOde()
{
// TODO: Add your control notification handler
code here
odeTest1();
odeTest2();
odeTest3();
odeTest4();
}
c.示例调试输出结果
| odeTest1();结果 |
odeTest2();结果 |
odeTest1
xn yn
0.0 1.0000000
0.1 1.0636613
0.2 1.1194055
0.3 1.1677203
0.4 1.2087883
0.5 1.2425414 |
odeTest2
xn yn
0.0 1.0000000
0.1 1.0636613
0.2 1.1194055
0.3 1.1677203
0.4 1.2087883
0.5 1.2425414 |
| odeTest3();结果 |
odeTest4();结果 |
odeTest3
xn y1n y2n y3n
0.00 0.0000000 1.0000000 1.0000000
0.01 0.0066665 0.9999667 0.9999830
0.02 0.0133323 0.9998889 0.9999433
0.03 0.0199970 0.9997667 0.9998810
0.04 0.0266601 0.9996001 0.9997961
0.05 0.0333212 0.9993891 0.9996885 |
odeTest3
xn
y1n
y2n
y3n
0.0 0.0000000
0.0000000 0.0000000
0.1 0.0634542
0.0634542 0.0634542
0.2 0.1228819
0.1228819 0.1228819
0.3 0.1785387
0.1785387 0.1785387
0.4 0.2306638
0.2306638 0.2306638
0.5 0.2794814
0.2794814 0.2794814 |
本文是作者原创，转载必须保证文章的完整性并标明出处（blog.sina.com.cn/zhangkunhn），请尊重作者，支持原创。
|
