数值分析第七章常微分方程的数值解法

最新推荐文章于 2025-06-28 00:47:47 发布

码农康康

最新推荐文章于 2025-06-28 00:47:47 发布

阅读量9k

点赞数 8

CC 4.0 BY-SA版权

分类专栏：计算机数学文章标签：计算机数学 Euler方程差分方程 Taylor展开式

本文链接：https://blog.youkuaiyun.com/ChiXueZhiHun/article/details/54755022

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

14 篇文章

订阅专栏

1 数值解法相关公式

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

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

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

{y' = f (x, y), a ⩽ x ⩽ b y (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) $f(x,y)$ 是已知函数,

α $\alpha$ 为给定的初值.
若函数

f(x,y) $f(x,y)$ 在区域

{a⩽x⩽b,−∞<y<+∞} $\left \{a \leqslant x \leqslant b, - \infty < y < +\infty \right \}$ 上连续且关于

y $y$ 满足

Lipschitz $Lipschitz$ 条件

| f (x, y) - f (x, y ¯) | ⩽ L | y - y ¯ |, \forall 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>0$ 为

Lipschitz $Lipschitz$ 常数，则初值问题（7.1）有唯一解。

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

假设初值问题(7.1)的解 $y=y(x)$ 唯一存在且足够光滑.对求解区域 $[a,b]$ 做剖分

a = x 0 < x 1 < x 2 < . . . < x n < . . . < x N = b

$a=x_0<x_1<x_2<...<x_n<...<x_N=b$ 其中剖分点

xn=a+nh,n=0,1,...,N $x_n=a+nh,n=0,1,...,N$ ,

h $h$ 称为剖分步长，数值解法就是求精确解

y(x) $y(x)$ 在剖分节点

xn $x_n$ 上的近似值

yn≈y(xn),n=1,2,...,N $y_n\approx y(x_n),n=1,2,...,N$ .
我们采用数值积分方法来建立差分公式.
在区间

[xn,xn+1] $[x_n,x_{n+1}]$ 上对方程(7.1)做积分，则有

y (x n + 1) - y (x n) = \int x n + 1 x n f (x, y (x)) d x

$y(x_{n+1})-y(x_n)=\int_{x_n} ^{x_{n+1}}f(x,y(x))dx$
对右侧分别应用左矩形公式、梯形公式和中矩形公式，可分别得到

Euler $Euler$ 公式、梯形差分公式和

Euler $Euler$ 中点公式。

1.4 $Euler$ 公式

{y n + 1 = y n + h f (x n, y n) y 0 = α, 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 梯形差分公式

⎧ ⎩ ⎨ y n + 1 = y n + h 2 [f (x n, y n) + f (x n + 1, y n + 1)] y 0 = α, 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 $Euler$ 中点公式(双步 $Euler$ 公式)

{y n + 1 = y n - 1 + 2 h f (x n, y n) y 0 = α, 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}$
在

Euler $Euler$ 公式和梯形公式中,为求得

yn+1 $y_{n+1}$ ,只需用到前一步的值

yn $y_n$ ,这种差分方法称为单步法,这是一种自开始方法.

Euler $Euler$ 中点公式则不然, 计算

yn+1 $y_{n+1}$ 时需用到前两步的值

yn $y_n$ ,

yn−1 $y_{n-1}$ ,称其为两步方法,两步以上的方法统称为多步法.在

Euler $Euler$ 公式和

Euler $Euler$ 中点公式中,需要计算的

yn+1 $y_{n+1}$ 已被显式表示出来,称这类差分公式为显式公式,而梯形公式中,需要计算的

yn+1 $y_{n+1}$ 隐含在等式两侧,称其为隐式公式.隐式公式中,每次计算

yn+1 $y_{n+1}$ 都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好.

1.7 改进的 $Euler$ 方法

从数值积分的角度来看,梯形差分公式计算数值解的精度要比 $Euler$ 公式好,但它属于隐式公式,不便于计算.实际上,常将 $Euler$ 公式与梯形差分公式结合使用,通常采用只迭代一次的算法:

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ y n + 1 = y n + h 2 (K 1 + K 2) K 1 = f (x n, y n) K 2 = f (x n + h, y n + h K 1) y 0 = α, 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$ 为非负整数.显然,阶数越高,方法的精度越高.
研究差分公式阶的重要手段是 $Taylor$ 展开式,一元函数和二元函数的 $Taylor$ 展开式为:

2.1 $Taylor$ 一元展开式

$y (x n + 1) = y (x n + h) = y (x n) + y' (x n) h + y '' ( x n ) 2 ! h 2 + y ''' ( x n ) 3 ! h 3 + O (h 4)$ $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' (x n) = f (x n, y (x n)) = f (x n, y n) = f n y'' (x n) = f' (x n, y (x n)) = \partial f n \partial x + \partial f n \partial y f n y''' (x n) = f'' (x n, y (x n)) = \partial 2 f n \partial x 2 + 2 \partial 2 f n \partial x \partial y f n + \partial 2 f n \partial y 2 f 2 n + \partial f n \partial x \partial f n \partial y + (\partial f n \partial y) 2 f n$ $\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 $Taylor$ 二元展开式

$f (x n + h, y n + k) = f (x n + y n) + \partial f ( x n , y n ) \partial x h + \partial f ( x n , y n ) \partial y k + 1 2 ! [\partial 2 f ( x n , y n ) \partial x 2 h 2 + 2 \partial 2 f ( x n , y n ) \partial x \partial y h k + \partial 2 f ( x n , y n ) \partial y 2 k 2] + O (h 3)$ $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)$
对 $Euler$ 方法，有
$y n + 1 y (x n + 1) = y n + h f (x n + y n) = y (x n + h) = y (x n) + y' (x n) h + y '' ( x n ) 2 ! h 2 + y ''' ( x n ) 3 ! h 3 + O (h 4) = y n + f (x n, y n) h + O (h 2)$ $\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)$ ,所以 $Euler$ 是一阶方法。
3 单步方法的收敛性与稳定性

3.1 单步方法的收敛性

对于给定的初值问题

${y' = f (x, y), a ⩽ x ⩽ b y (a) = α$ $\begin{equation} \left \{ \begin{aligned} & y'=f(x,y),a\leqslant x \leqslant b\\ &y(a)=\alpha \end{aligned} \right. \end{equation}$
的单步显示方法可以统一写成如下形式：
$y n + 1 = y n + h Φ (x n, y n, h) (7.1)$ $\begin{equation} y_{n+1}=y_n+h\Phi(x_n,y_n,h) \end{equation} \tag{7.1}$
其中， $\Phi(x,y,h)$ 称为增量函数，对于 $Euler$ 方法，有 $\Phi(x,y,h)=f(x,y)$ ,对于改进 $Euler$ 方法，有 $\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 \to 0 y n = y (x n)$ $\lim_{h \rightarrow0}y_n=y(x_n)$ ，则称单步法(7.1)是收敛的。
可见,若方法(8.5)是收敛的,则当 $h \rightarrow 0$ 时,整体截断误差 $e_n=y(x_n)-y_n$ 将趋于零.
定理7.1

设单步方法(7.1)是 $p\geqslant 1$ 阶方法, 增量函数 $\Phi(x,y,h)$ 在区域 $\left \{ a \leqslant x \leqslant b, -\infty < y < +\infty, 0 \leqslant h \leqslant h_0 \right \}$ 上连续,且关于y满足 $Lipschitz$ 条件,初始近似 $y_0=y(a)=\alpha$ ,则方法(7.1)是收敛的,且存在与h无关的常数 $C$ ,使 $\left | y(x_n)-y_n \right | \leqslant Ch^p$

3.2 单步方法的稳定性

定义7.2 单步方法的稳定性

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

说明：单步显式方法的稳定性与步长密切相关, 在一种步长下是稳定的差分公式,取大一点步长就可能是不稳定的.收敛性是反映差分公式本身的截断误差对数值解的影响;稳定性是反映计算过程中舍入误差对数值解的影响.只有即收敛又稳定的差分公式才有实用价值.