微分方程的数值解法之Runge-Kutta方法

简介

当我们在纸面上画曲线时会遇到三种情况: 第一种是曲线可以初等函数通过各种组合获得, 第二种是曲线可以由初等函数及其导函数确定,但是无法求出明确的函数关系,第三种就是没有任何关系的随机函数。而物理规律往往存在一定关系的,也就是各种各样的微分方程,但是有限的数学函数不可能构建出各种各样的物理规律,所以这就是根据基本微积分原理使用列表法罗列出自变量与因变量的对应关系,把这些分立的点在纸上描出,只要足够分离间隔足够小,基本上就能得到精确的数值解。然后根据曲线的变化趋势就可以找到物理规律了。这就包括了前两种情况,所以数值解法可以应用的范围极广,但是若有解析解,我们可以方便的通过一些数学手段来分析其规律,没有解析解,也就是没有初等函数组合出的具体函数,我们只能牺牲时间,来逐点计算,列表以找到映射关系。

Runge-Kutta 法的理论基础是 Taylor 展开,通过增加Taylor的幂级数项来增加精度,但是注意这隐含着那个不可解析的函数也必须是收敛的才有意义!这个需要根据具体情况判断,换言之出现无穷大的情况就要小心了。

理论细节

问题数学化

已经初始值:\(x_0\), \(y_0\) 和微分方程关系式 \(y'=f(x,y)\), 求\(y=y(x)\)的映射关系(列表法)

具体步骤

建立的基础是比卡逐次逼近法,不断增加Taylor级数的项。具体操作为:

  1. \(x_1=x_0+h\), 采用线性项做初步近似,即 \[\begin{equation}\label{eq:jie1} g_1(x_1)=y_0+hf(x_0,y_0) \end{equation}\]

  2. 若我们能求出二阶导数,则在式\(\eqref{eq:jie1}\)的基础上再追加上二阶项后记为\(g_2(x_1)\), 显然比式\(\eqref{eq:jie1}\)更加精确。按导数的定义,只要令\(h\)足够小,\((x_0,y_0)\)处的二阶导数就可以取为 \[\begin{equation}\label{eq:daoshu1} f'(x_0,y_0)=\frac{f(x_1,g_1(x_1))-f(x_0,y_0)}{h} \end{equation}\] 于是式\(\eqref{eq:jie1}\)按Taylor级数追加下一阶项后为 \[\begin{equation}\label{eq:jie2} g_2(x_1)=y_0+hf(x_0,y_0)+\frac{h^2}{2!}f'(x_1,g_1(x_1)) \end{equation}\]

  3. 按照上述步骤一直做下去,假设已经精确到了\(n\)阶项,即 \[\begin{equation}\label{eq:jieN} g_n(x_1)=y_0+hf(x_0,y_0)+\frac{h^2}{2!}f'(x_1,g_1(x_1))+\cdots + \frac{h^n}{n!}f^{(n)}(x_1,g_{n-1}(x_1)) \end{equation}\] 进一步可以求解第\(n+1\)阶导数 \[\begin{equation}\label{eq:jieN1} f^{(n+1)}(x_0,y_0)=\frac{f^{(n)}(x_1,g_n(x_1))-f^{(n)}(x_0,y_0)}{h} \end{equation}\] 将式\(\eqref{eq:jieN1}\)按Taylor展开式代入式\(\eqref{eq:jieN}\)\[\begin{align} g_{n+1}(x_1)&=y_0+hf(x_0,y_0)+\frac{h^2}{2!}f'(x_1,g_1(x_1))+\cdots +\\ &\frac{h^n}{n!}f^{(n)}(x_1,g_{n-1}(x_1))+\frac{h^{n+1}}{(n+1)!}f^{(n+1)}(x_1,g_{n}(x_1)) \label{eq:jieN2} \end{align}\]

  4. 按式\(\eqref{eq:jieN}\)到式\(\eqref{eq:jieN2}\)一直递推下去就可以得到越来越接近\(y(x_1)\)的值,所认也可以说近似精度在提高。在保持\(n\)阶近似时,可以取下一个函数点\((x_1,y_1)\), 其中\(y_1=g_n(x_1)\). 在此基础上,以\((x_1,y_1)\)为起点,将递推关系中的\(x_0\)换成\(x_1\),\(x_1\)换成\(x_2\), 重复递推过程则获得\((x_2,y_2)\), 一直执行下去就可以获得列表法的结果\((x_i,y_i)\), 其中\(i=0\),\(1\),\(2\),\(\cdots\)

  5. 上面分析了Runge-Kutta法的理论,并未直接从程序入手分析,也没有列出程序编码,其目的是理解数值计算的理论,在今后的学习和工作中更好的处理物理问题。

参考文章