双曲方程的有限差分法

一阶双曲方程 $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}$