数值分析(8) 常微分方程初值问题

常微分方程基本概念

常微分方程（ODE） ：未知函数为一元函数的微分方程

常微分方程的阶数表示方程中未知函数的最高阶导数次数

高阶常微分方程可以通过变量代换转化为一阶常微分方程组，例如：

$$y''+y'+y=0$$

可以转化为

$$\begin{cases}y'=u\\u'=-y-u\end{cases}$$

因此，只需考虑一阶常微分方程组，未知函数的变量为$t$时，可记作

$$\boldsymbol y'=\boldsymbol f(t,\boldsymbol y)$$

ODE的解是一个函数簇，一般需要添加一个条件来得到唯一确定的函数，此条件成为边界条件或初始条件，最常见的形式为：

$$y(t_0)=y_0$$

一些特殊的定义：

• 线性常微分方程：$y'=f(t,y)=a(t)y+b(t)$
• 线性齐次常微分方程：$y'=a(t)y$
• 线性齐次常系数常微分方程：$y'=\lambda y$

---

定义（ODE初值问题的稳定性）

对于一个ODE初值问题，若

• $t\to+\infty$时，初值扰动造成的$y(t)$的偏差有界，则称此初值问题稳定
• $t\to+\infty$时，初值扰动造成的$y(t)$偏差随扰动减小趋于零，则称此初值问题渐进稳定
  • 显然，渐进稳定是比稳定更强的条件
• $t\to+\infty$时，初值扰动造成的$y(t)$偏差发散，则称此初值问题不稳定

注 这里要求$t\to+\infty$是因为这是一般实际应用中考虑的情况，这里$t$的含义往往是时间

例 考虑$\begin{cases}y'=\lambda y\\y(t_0)=y_0\end{cases}$的稳定性

准确解为

$$y=y_0e^{\lambda(t-t_0)}$$

初值扰动$\delta$后，得到的解为

$$y=(y_0+\delta)e^{\lambda(t-t_0)}$$

于是偏差为

$$\Delta y=\delta e^{\lambda(t-t_0)}$$

显然，问题稳定等价于

$$\lambda\leq 0$$

---

对一般的非线性ODE方程组

$$\boldsymbol y'=\boldsymbol f(t,\boldsymbol y)$$

此方程组可以线性近似为

$$\boldsymbol y'=\mathbf J\boldsymbol y+\boldsymbol b(t)$$

其中$\mathbf J$为$\boldsymbol f$对$\boldsymbol y$的Jacobi矩阵

稳定性与$\boldsymbol b(t)$项无关，此时问题局部稳定等价于$\mathbf J$的所有特征值实部$\leq 0$

ODE初值问题的数值解法

一般地，初值问题通过对$t$离散化迭代计算来求解数值解，成为离散变量法

每一步迭代计算的$t$记为

$$t_0<t_1<t_2<\cdots$$

相邻$t$的间距$h_n=t_{n+1}-t_n$称为步长

实际近似解求解方程一般形如

$$y_{n+1}=G(y_{n+1}+y_n+y_{n-1}+\cdots+y_{n-k})$$

• 若$k=0$，称为单步法，否则称为多步法
• 若右端不含$y_{n+1}$，称为显格式方法，否则称为隐格式方法
  • 显格式方法可以直接计算，隐格式方法需要求解方程

欧拉法

为显式单步法，公式为：

$$y_{n+1}=y_n+h_nf(t_n,y_n)$$

这实际上直接使用数值微分的向前差分公式得到，即

$$y'=f(t,y)$$

$$\frac{y_{n+1}-y_n}{h_n}=f(t_n,y_n)$$

除了上面的微分理解，也可以使用积分理解：

$$y(t_{n+1})=y(t_n)+\int_{t_n}^{t_{n+1}}f(t,y){\rm d}t\approx y_n+h_nf(t_n,y_n)$$

这里相当于使用左矩形近似积分

数值解法的稳定性

上面我们讨论了ODE初值问题的稳定性，这是与具体方法无关的，这里讨论一个具体的数值解法的稳定性

定义 若在节点$t_n$上的近似值$y_n$有扰动$\delta_n$，则对$m>n$，此扰动引起的$y_m$的误差满足$|\delta_m|\leq|\delta_n|$，则称此方法稳定

---

例 考虑欧拉法解$y'=\lambda y$的稳定性（步长为常数）

$$y_{n+1}=y_n+h\lambda y_n=(1+h\lambda)y_n$$

于是

$$\delta_{n+1}=(1+h\lambda)\delta_n$$

欧拉法稳定等价于

$$|1+h\lambda|\leq 1$$

由于算法稳定一般首先要求问题稳定，我们要求$\lambda\leq 0$

于是

$$h\leq-\frac2\lambda$$

即：要想让欧拉法稳定，步长不能太大

---

定义（局部截断误差） 假设$y_{n+1}$计算中依赖的所有前序值 全部精确 ，则数值解法的局部截断误差定义为：

$$l_{n+1}=y(t_{n+1})-y_{n+1}$$

注 局部截断误差反映了计算中 一步 的误差

例 欧拉法解$y'=\lambda y$的局部截断误差

$$l_{n+1}=y(t_{n+1})-y_{n+1}=y(t_n)[e^{h\lambda}-1-h\lambda]=O(h^2)$$

对一般的$y'=f(t,y)$，有相同的结论：

$$l_{n+1}=O(h^2)$$

定义 若某种解法的局部截断误差$l_{n+1}=O(h^{p+1})$，则称该方法具有$p$阶准确度

注 局部截断误差的累计结果，即整体截断误差，比局部截断误差要低一阶。因此准确度定义也比局部截断误差低一阶

注 准确度表示的是迭代点不变的情况下，随$h$减小，误差的减小情况，至少具有1阶准确度表示方法是收敛的

变形欧拉法

欧拉法的积分理解为使用 左矩形 近似积分，但我们还有其他的做法：

• 向后欧拉法，使用 右矩形 近似积分，即$y_{n+1}=y_n+h_nf(t_{n+1},y_{n+1})$
• 梯形法，使用 梯形 近似积分，即$y_{n+1}=y_n+h_n/2[f(t_n,y_n)+f(t_{n+1}+y_{n+1})]$

两者均为单步隐格式方法

---

向后欧拉法分析

对$y'=\lambda y$，计算公式为

$$y_{n+1}=y_n+h\lambda y_{n+1}$$

即

$$y_{n+1}=\frac{1}{1-h\lambda}y_n$$

于是

$$\delta_{n+1}=\frac{1}{1-h\lambda}\delta_n$$

稳定条件为

$$|1-h\lambda|\geq 1$$

只要问题稳定，即$\lambda\leq0$，上述不等式必然成立，故向后欧拉法无条件稳定

其局部截断误差满足

$$l_{n+1}=y(t_{n+1})-y_{n+1}$$

$$=y(t_{n+1})-y_n-hf(t_{n+1},y_{n+1})$$

$$=y(t_{n+1})-y(t_n)-hf(t_{n+1},y_{n+1})$$

$$=hf(t_{n+1},y(t_{n+1}))-hf(t_{n+1},y_{n+1})+O(h^2)$$

$$=hf_y'(t_{n+1},\xi)[y(t_{n+1})-y_{n+1}]+O(h^2)$$

$$=hf_y'(t_{n+1},\xi)l_{n+1}+O(h^2)$$

故

$$l_{n+1}=\frac{1}{1-hf_y'(t_{n+1},\xi)}O(h^2)=O(h^2)$$

具有1阶准确度

梯形法类似，也具有无条件稳定性，但有2阶准确度

Runge-Kutta方法

引入

这是由欧拉法改进的一种显式方法，其思路为类似数值积分中的中矩形或梯形公式

如果使用中矩形公式，我们需要做：

$$y_{n+1}=y_n+\int_{t_n}^{t_{n+1}}f(t,y)$$

$$\approx y_n+f(t_n+\frac{h_n}2,y(t_n+\frac{h_n}2))$$

然而， 中点处的$y$值无法确定，这里使用欧拉法估算：

$$k_1=f(t_n,y_n)$$

$$k_2=f(t_n+\frac h2,y_n+\frac h2k_1)$$

$$y_{n+1}=y_n+hk_2$$

上面的方法称为中点公式

---

类似地，我们也可以先用欧拉法估算终点处的函数值：

$$k_1=f(t_n,y_n)$$

$$k_2=f(t_n+h,y_n+hk_1)$$

$$y_{n+1}=y_n+\frac h2(k_1+k_2)$$

上面的方法称为Heun方法或改进的欧拉法

---

中点公式和Heun方法的共同点是：

• 都使用了两次函数求值，其中一次使用了欧拉法估计
• 都比欧拉法准确
• 都属于2级R-K公式的特例

Runge-Kutta方法概述

Runge-Kutta是一般意义在的显式、单步数值法，我们从以下的数值积分出发：

$$y_{n+1}=y_n+\int_{t_n}^{t_n+h}f(s,y(s)){\rm d}s$$

使用某个机械求积公式，则

$$y_{n+1}=y_n+h\sum\limits_{i=1}^rc_if(t_n+\lambda_ih,y(t_n+\lambda_ih))$$

然而，式中的$y(t_n+\lambda_i h)$只能估算，设所有这些$f$值的估计值分别为$k_1,k_2,\cdots,k_r$，则得到迭代公式

$$y_{n+1}=y_n+h\sum\limits_{i=1}^rc_ik_i$$

---

一般的R-K方法描述如下：

• 首先在$[t_n,t_{n+1}]$上取$r$个不同的$t$值，此时的公式称为$r$级公式，估算每个$f(t,y)$
• 第一个估算值固定为$k_1=f(t_n,y_n)$
• 第二个估算值使用欧拉法估算，即$k_2=f(t_n+\lambda_2h,y_n+\lambda_2hk_1)$
• 后续估算值要综合前面所有信息，例如第三个值如下法计算：$k_3=f(t_n+\lambda_3h,y_n+\mu_{3,1}hk_1+\mu_{3,2}hk_2)$
  • 其中$\mu_{3,j}$为待定系数，这些系数选择需要时公式准确度阶数尽可能高
• 最终计算公式为$y_{n+1}=y_n+h\sum_{i=1}^rc_ik_i$

于是，上面介绍的中点公式和改进的欧拉法为$\lambda_2$取值不同的2级R-K公式

---

一种经典的4级、4阶的R-K法为：

$$k_1=f(t_n,y_n)$$

$$k_2=f(t_n+\frac h2,y_n+\frac h2k_1)$$

$$k_3=f(t_n+\frac h2,y_n+\frac h2k_2)$$

$$k_4=f(t_n+h,y_n+hk_3)$$

$$y_{n+1}=y_n+\frac h6(k_1+2k_2+2k_3+k_4)$$

---

R-K法参数的推导

以2级R-K公式为例，对以下公式

$$y_{n+1}=y_n+h(c_1k_1+c_2k_2)$$

$$k_1=f(t_n,y_n)$$

$$k_2=f(t_n+\lambda_2h,y_n+\lambda_2hk_1)$$

求$c_1,c_2,\lambda_2$使公式准确度阶数尽可能高

对模型问题$y'=\lambda y$，考虑到

$$y_{n+1}=y_n+c_1hk_1+c_2hk_2$$

$$=y_n+c_1h\lambda y_n+c_2h\lambda(y_n+\lambda_2h\lambda y_n)$$

$$=y(t_n)[1+(c_1+c_2)h\lambda+c_2h^2\lambda^2\lambda_2]$$

而

$$y(t_{n+1})=y(t_n)+hy'(t_n)+\cdots$$

$$=y(t_n)[1+h\lambda+\frac{h^2\lambda^2}{2}+\cdots]$$

此时最高只能达到2阶准确度，即$y(t_{n+1})=y_{n+1}+O(h^3)$，此时的条件为

$$c_1+c_2=1$$

$$c_2\lambda_2=\frac12$$

对中点公式而言

$$c_1=0,c_2=1,\lambda_2=\frac12$$

对改进的欧拉法而言

$$c_1=c_2=\frac12,\lambda_2=1$$

故它们都是有2阶准确度的2级R-K公式

---

事实上，不超过4级的R-K公式的准确度阶数等于其级数，而超过4级的R-K公式的准确度阶数无法达到级数，故一般使用的最高级公式就是4级R-K公式

---

显式单步法的相容性条件

一般的显式单步法形如

$$y_{n+1}=y_n+h\phi(t_n,y_n,h)$$

注 显然$\lim\limits_{h\to0^+}y_{n+1}=y_n$，故递增项应至少是$h$的一阶量

最简单的单步法即欧拉法具有1阶准确度，因而更复杂的单步法如果要有效，必须至少具有1阶准确度

对欧拉法，我们有

$$y(t_{n+1})=y(t_n)+hf(t_n,y(t_n))+O(h^2)$$

要想让

$$y(t_{n+1})=y(t_n)+h\phi(t_n,y(t_n),h)+O(h^2)$$

只需让

$$\phi(t_n,y(t_n),h)=f(t_n,y(t_n))+O(h)$$

一般情况下$\phi$都是$h$的连续函数，于是这等价于

$$\phi(t_n,y(t_n),0)=f(t_n,y(t_n))$$

上述条件是显式单步法有效的基本条件，称为相容性条件

对于常见的单步法形式：

$$y_{n+1}=y_n+h\sum\limits_{i=1}^rc_ik_i$$

$$k_i=f(t_n+\lambda_ih,y_n+hX)$$

相容性条件等价于

$$\sum\limits_{i=1}^rc_i=1$$