飞行仿真--5.微分方程的计算机数值求解方法

chunchun

以下文章来源：微信公众号    无风夜雪
飞行仿真数学建模必然涉及大量的微积分方程，以直升机为例，就涉及有机体六自由度微分方程、桨叶的挥舞运动方程等，若采用叶素法和基于动态入流理论则还需要解算动态入流方程。
《高等数学》、《数值分析》这一类的书尤其是《数值分析》对常见的微分方程的数学求解方法有一套较为成熟的理论。毫无疑问，这些理论归功于那些伟大的数学家们，但是实际更需要使用这套理论的往往是具体行业的工程师们，比如飞行仿真的建模者，他们需要熟练应用这些理论并采用计算机程序实现之，而不像数学家那样take a paper写出演算过程这么简单。
我们大多数人都是关心微分方程的解。其实只有少数简单的微分方程可以求得解析解（即有明确的数学解析式的结果）。不过在无法求得解析解时，可以利用数值分析的方式，利用电脑来找到其数值解。《数值分析》就是这么一门课程，介绍常见微分方程的数值解法。

同样的算法用计算机编程出来结果可能是一致的，但效率却大相径庭。最简单的例子，比如做1到100的算数和，有人喜欢写个for循环循环一百次累加得到结果，有的却能发现规律：1+99=2+98……=100，第二种循环次数明显减半。数学课本永远不可能告速我们怎么用C语言或C++亦或是java语言写具体的算法，它只能告诉我们所谓的"算法"，计算机程序员就没那么简单了，它必须充分理解其数学含义并将"算法"翻译成对应的程序代码，但从这一点讲我觉得我们就应该好好感谢这么多辛苦的程序员们。
如前文《飞行仿真--3.刚体六自由度方程、变换矩阵与四元素》所述的机体六自由度微分方程，就是典型的一阶常微分方程，工程界普遍认为在现在的计算机硬件水平下采用龙格库塔Runge-Kutta方法是最完美的。如四级四阶Runge-Kutta公式的局部截断误差是O(h5),即步长的5次方，已经能够满足一般的飞行仿真需要。
关于龙格库塔Runge-Kutta的具体计算机编程方法，我发现"坤玲2008"的一篇博文《通用龙格库塔Runge-Kutta方法求解常微分方程组初值问题的C++优雅实现》将复杂的龙格库塔Runge-Kutta方法解释的淋漓尽致。
该文不仅给出了龙格库塔Runge-Kutta的具体代码，更重要的是采用了C++标准库中泛型函数的写法，设计了一个比较通用的龙格库塔算法。有带参数和不带参数两种，根据自己需要很容易使用，不再废话，在此极力推荐：
http://blog.sina.com.cn/s/blog_6d4af96701019j3j.html
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)
}

yjy289439

你好 有飞行器仿真的源码吗   最近在研究这个  谢谢啦

FAAC

您好，请问采用c++而非c语言求解微分方程的原因是什么？能否深入解释一下？谢谢

风鸣天下

这个在程序届向来争议已久，应具体而论，这个例子里看似简单，但需要考虑实现几种不同的输入，但但算法基本一直。在C++只要用STL就可以轻松实现这重载，在C++中这叫多态，是面相对象编程的核心区特点之一。如果是用C那就是定义几个版本的函数了。另外数值计算离不开数学函数库，C++的数学函数库比较丰富，如果您有丰富的C数学函数库当然不妨一试。