【数值分析干货】第七章常微分方程的数值解法

原创已于 2025-08-17 13:31:13 修改 · 9.2k 阅读

34 ·

CC 4.0 BY-SA版权

文章标签：

#计算机数学 #Euler方程 #差分方程 #Taylor展开式

于 2017-01-27 12:08:25 首次发布

计算机数学专栏收录该内容

12 篇文章

订阅专栏

1 数值解法相关公式

1.1 为什么要研究数值解法？

所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.

1.2 问题 7.1 一阶常微分方程初值问题的一般形式

${y′=f(x,y),a⩽x⩽by(a)=α\begin{equation} \left \{ \begin{aligned} & y'=f(x,y),a\leqslant x \leqslant b\\ &y(a)=\alpha \end{aligned} \right. \end{equation}$
其中 $f (x, y)$ 是已知函数, $α\alpha$ 为给定的初值.
若函数 $f (x, y)$ 在区域
${a⩽x⩽b,−∞<y<+∞}\left \{a \leqslant x \leqslant b, - \infty < y < +\infty \right \}$ 上连续且关于 $y$ 满足 $L i p sc hi t z$ 条件 $∣f(x,y)−f(x,yˉ)∣⩽L∣y−yˉ∣,∀y,yˉ\left | f(x,y)-f(x, \bar y) \right | \leqslant L \left | y-\bar y \right |, \forall y, \bar y$
其中 $L > 0$ 为 $L i p sc hi t z$ 常数，则初值问题（7.1）有唯一解。

1.3 构造数值解法的基本思想

假设初值问题(7.1)的解 $y = y (x)$ 唯一存在且足够光滑.对求解区域 $[a, b]$ 做剖分 $a=x_0<x_1<x_2<...<x_n<...<x_N=b$ 其中剖分点 $x_n=a+nh,n=0,1,...,N$ , $h$ 称为剖分步长，数值解法就是求精确解 $y (x)$ 在剖分节点 $x_n$ 上的近似值 $yn≈y(xn),n=1,2,...,Ny_n\approx y(x_n),n=1,2,...,N$ .
我们采用数值积分方法来建立差分公式.
在区间 $x_n,x_{n+1}]$ 上对方程(7.1)做积分，则有
$y(xn+1)−y(xn)=∫xnxn+1f(x,y(x))dxy(x_{n+1})-y(x_n)=\int_{x_n} ^{x_{n+1}}f(x,y(x))dx$
对右侧分别应用左矩形公式、梯形公式和中矩形公式，可分别得到 $E u l er$ 公式、梯形差分公式和 $E u l er$ 中点公式。

1.4 $E u l er$ 公式

${yn+1=yn+hf(xn,yn)y0=α,n=0,1,2,...,N−1\begin{equation} \left \{ \begin{aligned} &y_{n+1}=y_n+hf(x_n,y_n) \\ &y_0=\alpha,n=0,1,2,...,N-1 \end{aligned} \right. \end{equation}$

1.5 梯形差分公式

${yn+1=yn+h2[f(xn,yn)+f(xn+1,yn+1)]y0=α,n=0,1,2,...,N−1\begin{equation} \left\{ \begin{aligned} &y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1})] \\ &y_0=\alpha,n=0,1,2,...,N-1 \\ \end{aligned} \right. \end{equation}$

1.6 $E u l er$ 中点公式(双步 $E u l er$ 公式)

${yn+1=yn−1+2hf(xn,yn)y0=α,n=0,1,2,...,N−1\begin{equation} \left\{ \begin{aligned} &y_{n+1}=y_{n-1}+2hf(x_n,y_n) \\ &y_0=\alpha,n=0,1,2,...,N-1 \end{aligned} \right. \end{equation}$
在 $E u l er$ 公式和梯形公式中,为求得 $y_{n+1}$ ,只需用到前一步的值 $y_n$ ,这种差分方法称为单步法,这是一种自开始方法. $E u l er$ 中点公式则不然, 计算 $y_{n+1}$ 时需用到前两步的值 $y_n$ , $y_{n-1}$ ,称其为两步方法,两步以上的方法统称为多步法.在 $E u l er$ 公式和 $E u l er$ 中点公式中,需要计算的 $y_{n+1}$ 已被显式表示出来,称这类差分公式为显式公式,而梯形公式中,需要计算的 $y_{n+1}$ 隐含在等式两侧,称其为隐式公式.隐式公式中,每次计算 $y_{n+1}$ 都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好.
###1.7 改进的 $E u l er$ 方法
从数值积分的角度来看,梯形差分公式计算数值解的精度要比 $E u l er$ 公式好,但它属于隐式公式,不便于计算.实际上,常将 $E u l er$ 公式与梯形差分公式结合使用,通常采用只迭代一次的算法:
${yn+1=yn+h2(K1+K2)K1=f(xn,yn)K2=f(xn+h,yn+hK1)y0=α,n=0,1,2,...,N−1\begin{equation} \left \{ \begin{aligned} &y_{n+1}=y_n + \frac{h}{2}(K_1+K_2) \\ &K_1=f(x_n,y_n) \\ &K_2=f(x_n+h,y_n+hK_1)\\ &y_0=\alpha,n=0,1,2,...,N-1 \end{aligned} \right. \end{equation}$

2 差分公式的误差分析

在节点 $x_{n+1}$ 的误差 $y(x_{n+1}-y_{n+1})$ ，不仅与 $y_{n+1}$ 这一步计算有关，而且与前 $n$ 步计算值 $y_n,y_{n-1},...,y_1$ 都有关。
为了简化误差分析，着重研究进行一步计算时产生的误差，即假设 $y_n=y(x_n)$ ，求误差 $y(x_{n+1})-y_{n+1}$ ,这时的误差称为局部截断误差，它可以反映出差分公式的精度。
如果单步差分公式的局部截断误差为 $O(h^{p+1})$ ,则称该公式为** $p$ 阶方法**.这里 $p$ 为非负整数.显然,阶数越高,方法的精度越高.
研究差分公式阶的重要手段是 $T a y l or$ 展开式,一元函数和二元函数的 $T a y l or$ 展开式为:

2.1 $T a y l or$ 一元展开式

$y(xn+1)=y(xn+h)=y(xn)+y′(xn)h+y′′(xn)2!h2+y′′′(xn)3!h3+O(h4)y(x_{n+1})=y(x_n+h)=y(x_n)+y'(x_n)h+\frac{y''(x_n)}{2!}h^2+\frac{y'''(x_n)}{3!}h^3+O(h^4)$
其中，
$y′(xn)=f(xn,y(xn))=f(xn,yn)=fny′′(xn)=f′(xn,y(xn))=∂fn∂x+∂fn∂yfny′′′(xn)=f′′(xn,y(xn))=∂2fn∂x2+2∂2fn∂x∂yfn+∂2fn∂y2fn2+∂fn∂x∂fn∂y+(∂fn∂y)2fn\begin{equation} \begin{aligned} &y'(x_n)=f(x_n,y(x_n))=f(x_n,y_n)=f_n \\ &y''(x_n)=f'(x_n,y(x_n))=\frac{\partial f_n}{\partial x} + \frac{\partial f_n}{\partial y}f_n \\ &y'''(x_n)=f''(x_n,y(x_n))=\frac{\partial ^2f_n}{\partial x^2} + \frac{2\partial ^2 f_n}{\partial x \partial y}f_n + \frac{\partial ^2f_n}{\partial y^2}f_n^2 + \frac{\partial f_n}{\partial x}\frac{\partial f_n}{\partial y} +\left ( \frac{\partial f_n}{\partial y} \right )^2 f_n \end{aligned} \end{equation}$

2.2 $T a y l or$ 二元展开式

$f(xn+h,yn+k)=f(xn+yn)+∂f(xn,yn)∂xh+∂f(xn,yn)∂yk+12![∂2f(xn,yn)∂x2h2+2∂2f(xn,yn)∂x∂yhk+∂2f(xn,yn)∂y2k2]+O(h3)f(x_n+h,y_n+k)=f(x_n+y_n)+\frac{\partial f(x_n,y_n)}{\partial x}h+\frac{\partial f(x_n,y_n)}{\partial y}k+\frac{1}{2!}\left [ \frac{\partial ^2f(x_n,y_n)}{\partial x^2}h^2+2\frac{\partial ^2f(x_n,y_n)}{\partial x \partial y}hk+\frac{\partial ^2 f(x_n,y_n)}{\partial y^2}k^2 \right ]+O(h^3)$
对 $E u l er$ 方法，有
$yn+1=yn+hf(xn+yn)y(xn+1)=y(xn+h)=y(xn)+y′(xn)h+y′′(xn)2!h2+y′′′(xn)3!h3+O(h4)=yn+f(xn,yn)h+O(h2)\begin{equation} \begin{aligned} y_{n+1}&=y_n+hf(x_n+y_n) \\ y(x_{n+1})&=y(x_n+h)=y(x_n)+y'(x_n)h+\frac{y''(x_n)}{2!}h^2+\frac{y'''(x_n)}{3!}h^3+O(h^4) \\ &=y_n+f(x_n,y_n)h+O(h^2) \end{aligned} \end{equation}$
从而 $y(x_{n+1})-y_{n+1}=O(h^2)$ ,所以 $E u l er$ 是一阶方法。

3 单步方法的收敛性与稳定性

3.1 单步方法的收敛性

对于给定的初值问题
${y′=f(x,y),a⩽x⩽by(a)=α\begin{equation} \left \{ \begin{aligned} & y'=f(x,y),a\leqslant x \leqslant b\\ &y(a)=\alpha \end{aligned} \right. \end{equation}$
的单步显示方法可以统一写成如下形式：
$yn+1=yn+hΦ(xn,yn,h)(7.1)\begin{equation} y_{n+1}=y_n+h\Phi(x_n,y_n,h) \end{equation} \tag{7.1}$
其中， $Φ(x,y,h)\Phi(x,y,h)$ 称为增量函数，对于 $E u l er$ 方法，有 $Φ(x,y,h)=f(x,y)\Phi(x,y,h)=f(x,y)$ ,对于改进 $E u l er$ 方法，有 $Φ(x,y,h)=12[f(x,y)+f(x+h,y+hf(x,y))]\Phi(x,y,h)=\frac{1}{2}[f(x,y)+f(x+h,y+hf(x,y))]$

定义7.1 单步法的收敛性

设 $y (x)$ 是初值问题(7.1)的解 , $y_n$ 是单步法 (7.1)产生的近似解.如果对任意固定的点 $x_n$ ,均有 $lim⁡h→0yn=y(xn)\lim_{h \rightarrow0}y_n=y(x_n)$ ，则称单步法(7.1)是收敛的。
可见,若方法(8.5)是收敛的,则当 $\rightarrow 0$ 时,整体截断误差 $e_n=y(x_n)-y_n$ 将趋于零.

定理7.1

设单步方法(7.1)是 $p⩾1p\geqslant 1$ 阶方法, 增量函数 $Φ(x,y,h)\Phi(x,y,h)$ 在区域 ${a⩽x⩽b,−∞<y<+∞,0⩽h⩽h0}\left \{ a \leqslant x \leqslant b, -\infty < y < +\infty, 0 \leqslant h \leqslant h_0 \right \}$ 上连续,且关于y满足 $L i p sc hi t z$ 条件,初始近似 $y0=y(a)=αy_0=y(a)=\alpha$ ,则方法(7.1)是收敛的,且存在与h无关的常数 $C$ ,使 $∣y(xn)−yn∣⩽Chp\left | y(x_n)-y_n \right | \leqslant Ch^p$

3.2 单步方法的稳定性

定义7.2 单步方法的稳定性

对于初值问题(7.1),取定步长 $h$ ,用某个差分方法进行计算时,假设只在一个节点值 $y_n$ 上产生计算误差 $δ\delta$ ,即计算值 $ynˉ=yn+δ\bar {y_n}=y_n+\delta$ ,如果这个误差引起以后各节点值 $y_m(m>n)$ 的变化均不超过 $δ\delta$ ,则称此差分方法是绝对稳定的.讨论数值方法的稳定性,通常仅限于典型的试验方程 $y′=λyy'=\lambda y$ ,其中 $λ\lambda$ 是复数且 $Re(λ)<0Re(\lambda)<0$ .在复平面上，当方法稳定时要求变量 $λh\lambda h$ 的取值范围称为方法的绝对稳定域，它与实轴的交集称为绝对稳定区间。