双曲方程的有限差分法

一阶双曲方程 $u_t+au_x=0$

特征方程与特征线：特征方程为 $\frac{dx}{dt}=a$ ，其特征线为 $x = at + C$ 。沿特征线 $x(t)=at + C$ ，对 $u(x(t),t)$ 求关于 $t$ 的全导数，由链式法则 $\frac{du}{dt}=u_x\cdot a + u_t = 0$ ，这表明 $u$ 沿着特征线的值是常数 。若已知初始条件 $u(x,0)=u_0(x)$ ，那么在特征线 $x = at + x_0$ 上，$u(x,t)=u(x - at,0)=u_0(x - at)$

二阶双曲方程 $u_{tt}-a^2u_{xx}=0$（波动方程 ）

特征方程与特征线：特征方程为 $\displaystyle (\frac{dx}{dt})^2 = a^2$ ，即 $\displaystyle \frac{dx}{dt}=\pm a$ ，特征线为 $x = at + C_1$ 和 $x=-at + C_2$ 。

变量变换与化简：令 $\xi = x - at$ ，$\eta = x + at$ ，对 $u$ 进行变量替换 $u(x,t)=\tilde{u}(\xi,\eta)=\tilde{u}(x - at,x + at)$ 。根据链式法则求偏导数
$$
\begin{align}
\frac{\partial u}{\partial t} &=-a\frac{\partial\tilde{u}}{\partial\xi}+a\frac{\partial\tilde{u}}{\partial\eta} \quad\quad
\frac{\partial^2 u}{\partial t^2} =a^2(\frac{\partial^2\tilde{u}}{\partial\xi^2}-2\frac{\partial^2\tilde{u}}{\partial\xi\partial\eta}+\frac{\partial^2\tilde{u}}{\partial\eta^2})\
\frac{\partial u}{\partial x} &=\frac{\partial\tilde{u}}{\partial\xi}+\frac{\partial\tilde{u}}{\partial\eta}\quad \quad \quad\quad\frac{\partial^2 u}{\partial x^2}=\frac{\partial^2\tilde{u}}{\partial\xi^2}+2\frac{\partial^2\tilde{u}}{\partial\xi\partial\eta}+\frac{\partial^2\tilde{u}}{\partial\eta^2}

\end{align}
$$
代入原方程
$$
\frac{\partial^2 u}{\partial t^2}=a^2\frac{\partial^2 u}{\partial x^2}
$$
化简可得 $\displaystyle \frac{\partial^2\tilde{u}}{\partial\xi\partial\eta}=0$

求解方程：对 $\displaystyle \frac{\partial^2\tilde{u}}{\partial\xi\partial\eta}=0$ 关于 $\xi$ 积分得 $\displaystyle \frac{\partial\tilde{u}}{\partial\xi}=\overline{f}_1(\xi)$ ，再关于 $\xi$ 积分得
$$
\begin{align}
\tilde{u} &=\int\overline{f}1(\xi)d\xi + f_2(\eta)\
&=f_1(\xi)+f_2(\eta)
\end{align}
$$
即
$$
u(x,t)=f_1(x - at)+f_2(x + at)
$$
结合初始条件 $u(x,0)=\varphi(x)$ 和 $u_t(x,0)=\psi(x)$ ，通过求导和积分运算确定
$$
\begin{align}
f_1(x)=\frac{1}{2}\varphi(x)-\frac{1}{2a}\int_0^x\psi(s)ds\
f_2(x)=\frac{1}{2}\varphi(x)+\frac{1}{2a}\int_0^x\psi(s)ds
\end{align}
$$
最终得到 d′Alembert公式
$$
u(x,t)=\frac{1}{2}[\varphi(x - at)+\varphi(x + at)]+\frac{1}{2a}\int{x - at}^{x + at}\psi(s)ds
$$
d′Alembert 公式的依赖域、影响域、决定域

依赖域：点 $(x_0,t_0)$ 的依赖域是初始线 $t = 0$ 上由特征线 $x = x_0 - at_0$ 和 $x = x_0 + at_0$ 所界定的区间 $[x_0 - at_0,x_0 + at_0]$ ，解在该点的值仅依赖于初始条件在这个区间的部分

影响域：初始点 $(x_0,0)$ 的影响域是由特征线 $x = x_0 + at$ 和 $x = x_0 - at$ 所界定的区域，初始条件在该点的扰动会影响这个区域内解的值

决定域：初始区间 $[x_1,x_2]$ 的决定域是由特征线 $x = x_1 + at$ 、$x = x_2 - at$ 所围的区域，解在这个区域内的值由初始区间 $[x_1,x_2]$ 上的初始条件决定

二阶波动方程显式差分格式

对于波动方程 $u_{tt}=a^2u_{xx}$ ，初值条件 $u(x,0)=\varphi(x)$ ，$u_t(x,0)=\psi(x)$ ，空间步长 $h$ ，时间步长 $\tau$ ，网格点 $x_j = jh$ ，$t_n = n\tau$

显式差分格式为
$$
\frac{u_j^{n + 1}-2u_j^n + u_j^{n - 1}}{\tau^2}=a^2\frac{u_{j + 1}^n - 2u_j^n + u_{j - 1}^n}{h^2}
$$
截断误差为 $O(\tau^2 + h^2)$ ，令 $\displaystyle r=\frac{a\tau}{h}$ （网格比 ），格式可整理为
$$
u_j^{n + 1}=r^2u_{j + 1}^n+(2 - 2r^2)u_j^n + r^2u_{j - 1}^n - u_j^{n - 1}
$$
初始条件为 $u_j^0=\varphi(x_j)$

$u_j^1$ 的计算

一阶近似：由
$$
\frac{u_j^1 - u_j^0}{\tau}=\psi(x_j)+O(\tau)
$$
可得
$$
u_j^1 = u_j^0+\tau\psi(x_j)+O(\tau^2)
$$
这是利用一阶向前差分近似时间导数 $u_t$

二阶近似：由
$$
\frac{u_j^1 - u_j^{-1}}{2\tau}=\psi(x_j)+O(\tau^2)
$$
结合波动方程的显式差分格式
$$
\frac{u_j^1 - 2u_j^0 + u_j^{-1}}{\tau^2}=a^2\frac{u_{j + 1}^0 - 2u_j^0 + u_{j - 1}^0}{h^2}
$$
消去 $u_j^{-1}$ 后得到
$$
u_j^1=\frac{r^2}{2}u_{j + 1}^0+(1 - r^2)u_j^0+\frac{r^2}{2}u_{j - 1}^0+\tau\psi(x_j)
$$

其中 $\displaystyle r = \frac{a\tau}{h}$ ，此为二阶精度的 $u_j^1$ 计算式，通过联立方程消元，利用差分格式的二阶精度要求推导得出

显式格式的稳定性分析（变量替换法 ）

一阶导数替换：令 $\displaystyle v=\frac{\partial u}{\partial t}$ ，则波动方程 $u_{tt}=a^2u_{xx}$ 可转化为一阶方程组
$$
\begin{cases}\frac{\partial u}{\partial t}=v\\frac{\partial v}{\partial t}=a^2\frac{\partial^2 u}{\partial x^2}\end{cases}
$$
其矩阵形式为
$$
\frac{\partial}{\partial t}\begin{bmatrix}u\v\end{bmatrix}=\begin{bmatrix}0&1\a^2\frac{\partial^2}{\partial x^2}&0\end{bmatrix}\begin{bmatrix}u\v\end{bmatrix}
$$
二阶导数替换（另一种形式 ）：令 $v=\frac{\partial u}{\partial t}$ ，$w = a\frac{\partial u}{\partial x}$ ，则方程组变为
$$
\begin{cases}\frac{\partial v}{\partial t}=a\frac{\partial w}{\partial x}\\frac{\partial w}{\partial t}=a\frac{\partial v}{\partial x}\end{cases}
$$
矩阵形式为
$$
\frac{\partial}{\partial t}\begin{bmatrix}v\w\end{bmatrix}=\begin{bmatrix}0&a\frac{\partial}{\partial x}\a\frac{\partial}{\partial x}&0\end{bmatrix}\begin{bmatrix}v\w\end{bmatrix}
$$
通过将波动方程的二阶显式差分格式
$$
\frac{u_j^{n + 1}-2u_j^n + u_j^{n - 1}}{\tau^2}=a^2\frac{u_{j + 1}^n - 2u_j^n + u_{j - 1}^n}{h^2}
$$
转化为关于时间一阶差分的形式，如
$$
\frac{\frac{u_j^{n + 1}-u_j^n}{\tau}-\frac{u_j^n - u_j^{n - 1}}{\tau}}{\tau}=a^2\frac{\frac{u_{j + 1}^n - u_j^n}{h}-\frac{u_j^n - u_{j - 1}^n}{h}}{h}
$$
令
$$
v_j^n=\frac{u_j^n - u_j^{n - 1}}{\tau} \quad \quad w_{j-\frac{1}{2}}^n=a\frac{u_j^n - u_{j - 1}^n}{h}
$$
将波动方程转化为一阶方程组
$$
\begin{cases}v_j^{n + 1}=r(w_{j+\frac{1}{2}}^n - w_{j-\frac{1}{2}}^n)+v_j^n\w_{j-\frac{1}{2}}^{n + 1}=r(v_j^{n + 1}-v_{j - 1}^{n + 1})+w_{j-\frac{1}{2}}^n\end{cases}
$$
其中 $\displaystyle r = \frac{a\tau}{h}$ ，这是基于对波动方程时间和空间导数的差分近似，把二阶方程转化为一阶方程组形式，便于分析稳定性

Fourier 方法分析稳定性

设周期边值条件下的通项解
$$
\begin{align}
v_j^n &= V_1^n e^{i\lambda x_j}\
w_{j}^n &= V_2^n e^{i\lambda x_{j}}
\end{align}
$$
$\displaystyle \lambda=\frac{2p\pi}{L}$ ，$p$ 为整数

将通项解代入差分方程组，经过化简得到关于 $V_1^{n + 1}$ 和 $V_2^{n + 1}$ 的递推关系
$$
\begin{cases}
\begin{align}
V_1^{n + 1}&=icV_2^n + V_1^n\
V_2^{n + 1}&=icV_1^{n + 1}+V_2^n\
&=icV_1^{n }+(1-c^2)V_2^n
\end{align}
\end{cases}
$$
其中 $\displaystyle c = 2r\sin\frac{\lambda h}{2}$

写成矩阵向量形式
$$
\begin{bmatrix}v^{n + 1}\w^{n + 1}\end{bmatrix}=\begin{bmatrix}1&ic\ic&1 - c^2\end{bmatrix}\begin{bmatrix}v^n\w^n\end{bmatrix}=G\begin{bmatrix}v^n\w^n\end{bmatrix}
$$

其中 $\displaystyle c = 2r\sin\frac{\lambda h}{2}$ ，$\displaystyle r=\frac{a\tau}{h}$

显式格式稳定等价于 ${G^n}$ 一致有界，当 $\displaystyle r=\frac{a\tau}{h}\leq1$ 时满足，此为 CFL（Courant - Friedrichs - Lewy）条件

CFL 条件几何解释：数值解依赖域 $[x_j - n\tau,x_j + n\tau]$ 需包含精确解依赖域 $[x_j - a n\tau,x_j + a n\tau]$ 。若数值依赖域小于精确依赖域（如 $\overline{P’Q’}>\overline{PQ}$ ），初始条件在 $\overline{P’Q’}\setminus\overline{PQ}$ 部分改变时，数值解不变但精确解改变，不符合实际。故需 $\overline{PQ}\geq\overline{P’Q’}$ ，推导出 $\displaystyle r=\frac{a\tau}{h}\leq1$ ，保证数值解能捕捉精确解的依赖关系

隐式格式

格式形式：用 $n - 1,n,n + 1$ 层中心差分加权平均逼近 $u_{xx}$ ，格式为
$$
\frac{u_j^{n + 1}-2u_j^n + u_j^{n - 1}}{\tau^2}=\frac{a^2}{h^2}[\theta\delta_x^2u_j^{n + 1}+(1 - 2\theta)\delta_x^2u_j^n+\theta\delta_x^2u_j^{n - 1}] \quad (0\leq\theta\leq1)
$$
整理得
$$
u_j^{n + 1}-r^2\theta\delta_x^2u_j^{n + 1}=(2 + r^2(1 - 2\theta)\delta_x^2)u_j^n+(r^2\theta\delta_x^2 - 1)u_j^{n - 1}
$$
当 $\theta = 0$ 时退化为显式格式；$\theta>0$ 为隐式格式

计算 $\boldsymbol{\theta=\frac{1}{4}}$ 时隐式格式稳定性（Fourier 方法）

格式推导与变量替换

针对波动方程隐式格式，当 $\theta=\frac{1}{4}$
$$
\frac{\displaystyle \frac{u_j^{n+1} - u_j^n}{\tau} - \frac{u_j^n - u_j^{n-1}}{\tau}}{\tau} = a^2 \cdot \frac{1}{4} \cdot \left[
\frac{\displaystyle \frac{u_{j+1}^{n+1} - u_j^{n+1}}{h} - \frac{u_j^{n+1} - u_{j-1}^{n+1}}{h}}{h} + \frac{1}{2} \cdot \frac{\displaystyle \frac{u_{j+1}^{n+1} - u_j^{n+1}}{h} - \frac{u_j^{n+1} - u_{j-1}^{n+1}}{h}}{h} + \frac{1}{4} \cdot \frac{\displaystyle \frac{u_{j+1}^{n+1} - u_j^{n+1}}{h} - \frac{u_j^{n+1} - u_{j-1}^{n+1}}{h}}{h}
\right]
$$

$$
= a \cdot \frac{1}{2} \cdot \frac{
a \cdot \left( \frac{1}{2} \cdot \frac{u_{j+1}^{n+1} - u_j^{n+1}}{h} + \frac{1}{2} \cdot \frac{u_{j+1}^n - u_j^n}{h} \right) - a \cdot \left( \frac{1}{2} \cdot \frac{u_j^{n+1} - u_{j-1}^{n+1}}{h} + \frac{1}{2} \cdot \frac{u_j^n - u_{j-1}^n}{h} \right)
}{h} \

  • a \cdot \frac{1}{2} \cdot \frac{
    a \cdot \left( \frac{1}{2} \cdot \frac{u_{j+1}^{n+1} - u_j^{n+1}}{h} + \frac{1}{2} \cdot \frac{u_{j+1}^n - u_j^n}{h} \right) - a \cdot \left( \frac{1}{2} \cdot \frac{u_j^{n+1} - u_{j-1}^{n+1}}{h} + \frac{1}{2} \cdot \frac{u_j^n - u_{j-1}^n}{h} \right)
    }{h}
    $$

令
$$
V_j^n=\frac{u_j^n - u_j^{n - 1}}{\tau}\quad \quad w_{j+\frac{1}{2}}^n=a(\frac{1}{2}\frac{u_{j + 1}^n - u_j^n}{h}+\frac{1}{2}\frac{u_{j + 1}^{n - 1}-u_j^{n - 1}}{h})
$$
则差分方程组化为
$$
\frac{V_j^{n + 1} - V_j^n}{\tau}
= a \left[
\frac{1}{2} \cdot \frac{W_{j + \frac{1}{2}}^{n + 1} - W_{j - \frac{1}{2}}^{n + 1}}{h}

  • \frac{1}{2} \cdot \frac{W_{j + \frac{1}{2}}^n - W_{j - \frac{1}{2}}^n}{h}
    \right]
    $$

$$
\frac{W_{j - \frac{1}{2}}^{n + 1} - W_{j - \frac{1}{2}}^n}{\tau}
= a \left[
\frac{1}{2} \cdot \frac{V_j^{n + 1} - V_{j - 1}^{n + 1}}{h}

  • \frac{1}{2} \cdot \frac{V_j^n - V_{j - 1}^n}{h}
    \right]
    $$

Fourier 方法分析

设周期解形式
$$
V_j^n = v^n e^{i\lambda x_j}\quad w_{j}^n = w^n e^{i\lambda x_{j}}
$$
代入差分方程组，化简得到增长矩阵
$$
G_n=\frac{1}{1 + \frac{c^2}{4}}\begin{bmatrix}1 - \frac{c^2}{4}&ic\ic&1 - \frac{c^2}{4}\end{bmatrix}
$$
其中 $\displaystyle c = 2r\sin\frac{\lambda h}{2}$ ，$\displaystyle r=\frac{a\tau}{h}$

计算增长矩阵的特征值，可得 $G^HG=1$ ，且矩阵的 $L^2$ 范数 $|G_n|_2 = 1$ 。根据稳定性判定，增长矩阵的幂次 $G_n^n$ 一致有界，故当 $\theta=\frac{1}{4}$ 时，该隐式差分格式恒稳定