CS205A Lecture 17 Time-stepping strategies

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第十七讲，这一讲继续介绍常微分方程（ODE）。

梯形方法

首先回顾一些基本定义：

$$y_i =y(t+ih)$$

利用中点差分公式不难得到

$${y}_{k+1}={y}_{k}+h F\left({y}_{k+1 / 2}\right)+O\left(h^{3}\right)$$

利用梯形法可得

$$\frac{F\left({y}_{k+1}\right)+F\left({y}_{k}\right)}{2}=F\left({y}_{k+1 / 2}\right)+O\left(h^{2}\right)$$

所以得到梯形方法计算公式

$${y}_{k+1}={y}_{k}+h \frac{F\left({y}_{k+1}\right)+F\left({y}_{k}\right)}{2}$$

类似之前讨论，继续对$y’=ay,a<0$分析该算法稳定性，带入可得

$$\begin{aligned} y_{k+1}&=y_{k}+\frac{1}{2} h a\left(y_{k+1}+y_{k}\right)\\ (1-\frac 1 2 ha) y_{k+1}&= (1+\frac 1 2 ha)y_k \\ y_{k+1}&= \frac{1+\frac 1 2 ha}{1-\frac 1 2 ha} y_k \end{aligned}$$

因此

$$y_{k}=\left(\frac{1+\frac{1}{2} h a}{1-\frac{1}{2} h a}\right)^{k} y_{0}$$

当满足下式的时候该算法稳定：

$$\left|\frac{1+\frac{1}{2} h a}{1-\frac{1}{2} h a}\right|<1$$

解得

$$ha <0$$

因为$a<0,h>0$，所以该算法总是稳定的。

龙格库塔法

首先由积分定义可得：

$$\begin{aligned} {y}_{k+1} &={y}_{k}+\int_{t_{k}}^{t_{k}+\Delta t} {y}^{\prime}(t) d t \\ &={y}_{k}+\int_{t_{k}}^{t_{k}+\Delta t} F[{y}(t)] d t \end{aligned}$$

如果使用梯形公式计算上述积分，那么得到

$${y}_{k+1}={y}_{k}+\frac{h}{2}(F\left({y}_{k}\right)+F\left({y}_{k+1}\right)) + O(h^3)$$

此即为之前推导的公式。注意我们有

$$y_{k+1}=y_k +hF(y_k)+O(h^2)$$

代入上式得到

$${y}_{k+1}={y}_{k}+\frac{h}{2}(F\left({y}_{k}\right)+F\left(y_k +hF(y_k)\right)) + O(h^3)$$

忽略$O(h^3)$的项得到Heun方法。依旧对$y’=ay,a<0$分析该算法稳定性：

$$\begin{aligned} y_{k+1} &=y_{k}+\frac{h}{2}\left(a y_{k}+a\left(y_{k}+h a y_{k}\right)\right) \\ &=\left(1+\frac{h}{2} a(2+h a)\right) y_{k} \\ &=\left(1+h a+\frac{1}{2} h^{2} a^{2}\right) y_{k} \end{aligned}$$

如果满足如下条件，那么该算法稳定：

$$\begin{aligned}-1 & \leq 1+h a+\frac{1}{2} h^{2} a^{2} \leq 1 \\ \Longleftrightarrow-4 & \leq 2 h a+h^{2} a^{2} \leq 0 \end{aligned}$$

因为$h>0$，所以左边的不等号恒成立，右边的不等号成立当且仅当$h\le \frac 2 {|h|}$。

Heun方法是龙格库塔法的一个例子，龙格库塔法的思路是逼近$F({y}_k+lh)$，最著名的是四次龙格库塔法：

$$\begin{aligned} {y}_{k+1} &={y}_{k}+\frac{h}{6}\left({k}_{1}+2 {k}_{2}+2 {k}_{3}+{k}_{4}\right) \\ \text { where } {k}_{1} &=F\left({y}_{k}\right) \\ {k}_{2} &=F\left({y}_{k}+\frac{1}{2} h {k}_{1}\right) \\ {k}_{3} &=F\left({y}_{k}+\frac{1}{2} h {k}_{2}\right) \\ {k}_{4} &=F\left({y}_{k}+h {k}_{3}\right) \end{aligned}$$

指数法

考虑如下方程

$$\vec{y}^{\prime}=A \vec{y}+G[\vec{y}]$$

移项得到

$$\vec y' -A\vec y =G[\vec y]$$

利用ODE的方法（积分因子），左乘$e^{-At}​$，化简可得

$$e^{-A t}\left(\vec{y}^{\prime}-A \vec{y}\right)=\frac{d}{d t}\left(e^{-A t} \vec{y}(t)\right)$$

所以

$$\frac{d}{d t}\left(e^{-A t} \vec{y}(t)\right)= \int_{0}^{t} e^{-A \tau} G[\vec{y}(\tau)] d \tau$$

利用$Ae^{At}=e^{At}A$，对左边积分得到

$$e^{-A t} \vec{y}(t)-\vec{y}(0)=\int_{0}^{t} e^{-A \tau} G[\vec{y}(\tau)] d \tau$$

即

$$\begin{aligned} \vec{y}(t) &=e^{A t} \vec{y}(0)+e^{A t} \int_{0}^{t} e^{-A \tau} G[\vec{y}(\tau)] d \tau \\ &=e^{A t} \vec{y}(0)+\int_{0}^{t} e^{A(t-\tau)} G[\vec{y}(\tau)] d \tau \end{aligned}$$

推广上述方法不难得到

$$\vec{y}_{k+1}=e^{A h} \vec{y}_{k}+\int_{t_{k}}^{t_{k}+h} e^{A\left(t_{k}+h-t\right)} G[\vec{y}(t)] d t$$

如果近似$G[\vec{y}(t)] \approx G\left[\vec{y}_{k}\right]​$，那么上式可以化为

$$\vec{y}_{k+1} \approx e^{A h} \vec{y}_{k}+\left[\int_{0}^{h} e^{A(h-t)} d t\right] G\left[\vec{y}_{k}\right]$$

所以

$$\vec{y}_{k+1}=e^{A h} \vec{y}_{k}+A^{-1}\left(e^{A h}-I_{n \times n}\right) G\left[\vec{y}_{k}\right]$$

多元情形

Newmark Schemes见讲义，这里忽略。