# 波动方程的数值求解

## December 4, 2021 • 铃宕 • 阅读设置

$$\partial_{t t} \phi-c^{2} \partial_{x x} \phi=0$$

#### 1.a 数学求解：hyperbolic PDE，解法无聊但适用范围广

$$A \frac{\partial^{2} \phi}{\partial x^{2}}+2 B \frac{\partial^{2} \phi}{\partial x \partial y}+C \frac{\partial^{2} \phi}{\partial y^{2}}+D \frac{\partial \phi}{\partial x}+E \frac{\partial \phi}{\partial y}+F \phi=G$$

• $\Delta > 0$ 双曲型方程

$$\frac{\partial^{2} \phi}{\partial x^{2}} - \frac{\partial^{2} \phi}{\partial y^{2}} = [...] \quad 或 \quad \frac{\partial^{2} \phi}{\partial x \partial y} = [...]$$

• $\Delta = 0$ 抛物线型方程

$$\frac{\partial^{2} \phi}{\partial x^{2}} = [...] \quad 或 \quad \frac{\partial^{2} \phi}{\partial y ^2} = [...]$$

• $\Delta < 0$ 椭圆形方程

$$\frac{\partial^{2} \phi}{\partial x^{2}} + \frac{\partial^{2} \phi}{\partial y^{2}} = [...]$$

$$x- \frac{B- \sqrt{B^{2}-AC}}{A}t= \gamma_1 \ \ \to x+ct = \gamma_1= \xi$$

$$x- \frac{B+ \sqrt{B^{2}-AC}}{A}t= \gamma_{2} \ \ \to x-cx = \gamma_2 = \eta$$

$$\phi _{t}=w_{\xi} \xi_{t}+w_{\eta} \eta_{t}=c\left(w_{\xi}-w_{\eta}\right) \\ \phi _{x}=w_{\xi} \xi_{x}+w_{\eta} \eta_{x}=w_{\xi}+w_{\eta}$$

$$\phi_{t t}-c^{2} \phi_{x x}=-4 c^{2} w_{\xi \eta}=0$$

$$w(\xi, \eta)=F(\xi)+G(\eta)$$

$$\phi (x,t)=F(x+ct)+G(x-ct)$$

$G(x-ct):$ forward wave $\qquad \quad F(x+ct):$ backward wave

#### 1.b 物理求解（by 朗道）

$$\left(\frac{\partial}{\partial t}-c \frac{\partial}{\partial x}\right)\left(\frac{\partial}{\partial t}+c \frac{\partial}{\partial x}\right) \phi=0$$

$$\frac{\partial}{\partial \xi}=\frac{1}{2}\left(\frac{\partial}{\partial t}-c \frac{\partial}{\partial x}\right), \quad \frac{\partial}{\partial \eta}=\frac{1}{2}\left(\frac{\partial}{\partial t}+c \frac{\partial}{\partial x}\right)$$

$$\frac{\partial^{2} \phi}{\partial \xi \partial \eta}=0$$

$$\phi = \phi_1(\xi) + \phi_2(\eta)$$

$$\phi (x,t)=F(x+ct)+G(x-ct)$$

$$x = const + ct$$

$$\phi(x,t=0)=\phi(x-ct,t)$$

#### 2. D'Alembert's formula

\begin{aligned} \partial_{t} \phi &=\Pi \\ \partial_{t} \Pi &=c^{2} \partial_{x x} \phi \end{aligned}

$$\phi(x,t=0)=\theta(x) \\ \Pi(x,t=0)=\psi(x)$$

$$\phi (x,0) = F(x) + G(x) = \theta (x)$$

$$\phi_{t}(x,0)=c F^{\prime}(x)-c G^{\prime}(x)=\psi (x)$$

$$F(x)-G(x)=\frac{1}{c} \int_{0}^{x} \psi(s) \mathrm{d} s+C$$

\begin{aligned} &F(x)=\frac{1}{2} \theta (x)+\frac{1}{2 c} \int_{0}^{x} \psi (s) \mathrm{d} s+\frac{C}{2} \\ &G(x)=\frac{1}{2} \theta (x)-\frac{1}{2 c} \int_{0}^{x} \psi (s) \mathrm{d} s-\frac{C}{2} \end{aligned}

$$\phi(x, t)=\frac{\theta (x+c t)+\theta (x-c t)}{2}+\frac{1}{2 c} \int_{x-c t}^{x+c t} \psi (s) \mathrm{d} s \quad \text{d’Alembert’s formula}$$

$$\phi (x,0)= \theta (x) \\ \phi_{t} (x,0)= \psi(x)=0$$

In this case, the d’Alembert’s solution takes the form

$$\phi(x, t)=\frac{1}{2}[\theta (x+c t)+\theta (x-c t)]$$

$$\phi (x,0)= \theta (x)=0 \\ \phi_{t} (x,0)= \psi(x)$$

For this case the d’Alembert’s solution assumes the form

$$\phi (x,t)=\frac{1}{2 c} \int_{x-c t}^{x+c t} \psi (s) \mathrm{d} s= \frac{1}{2}[g(x+ct) − g(x−ct)]$$

———————————————下面进入数值求解部分—————————————————

#### 3. Finite differencing approximation 有限差分近似

$$u_{i+1}= u_{i} + \left(\frac{\partial u}{\partial x}\right)_{i}(\Delta x)+\frac{1}{2 !}\left(\frac{\partial^{2} u}{\partial x^{2}}\right)_{i}(\Delta x)^{2}+\frac{1}{3 !}\left(\frac{\partial^{3} u}{\partial x^{3}}\right)_{i}(\Delta x)^{3}+\frac{1}{4 !}\left(\frac{\partial^{4} u}{\partial x^{4}}\right)_{i}(\Delta x)^{4} +O\left(\Delta x^{5}\right)$$

$$u_{i-1}=u_{i}+\left(\frac{\partial u}{\partial x}\right)_{i}(-\Delta x)+\frac{1}{2 !}\left(\frac{\partial^{2} u}{\partial x^{2}}\right)_{i}(-\Delta x)^{2}+\frac{1}{3 !}\left(\frac{\partial^{3} u}{\partial x^{3}}\right)_{i}(-\Delta x)^{3}+\frac{1}{4 !}\left(\frac{\partial^{4} u}{\partial x^{4}}\right)_{i}(-\Delta x)^{4} +O\left(\Delta x^{5}\right)$$

$$\frac{u_{i+1}-u_{i-1}}{2 h} = \frac{u'_{i}(2h) +\frac{1}{6} u_{i}^{'''}(2h^3)+O(h^{4})}{2h} \qquad \to \ \frac{du}{dx}=u'_{i}\approx\frac{u_{i+1}-u_{i-1}}{2 h} - \frac{1}{6}u_{i}^{'''}h^{2}$$

$$\frac{u_{i+1}-2 u_{i}+u_{i-1}}{h^{2}} = \frac{\frac{1}{2}u^{''}_{i}(2h^{2}) +\frac{1}{24} u_{i}^{''''}(2h^{4})+O(h^{5})}{h^{2}}$$

$$\to \frac{d^{2} u}{d x^{2}}=u^{''}_{i} \approx \frac{u_{i+1}-2 u_{i}+u_{i-1}}{h^{2}}-\frac{1}{12} u^{''''}_{i} h^{2}$$

$$u_{i+1, j}=u_{i, j}+\left(\frac{\partial u}{\partial x}\right)_{i, j} \Delta x+\left(\frac{\partial^{2} u}{\partial x^{2}}\right)_{i, j} \frac{(\Delta x)^{2}}{2}+\left(\frac{\partial^{3} u}{\partial x^{3}}\right)_{i, j} \frac{(\Delta x)^{3}}{6}+\cdots$$

$$\left(\frac{\partial u}{\partial y}\right)_{i+1, j}=\left(\frac{\partial u}{\partial y}\right)_{i, j}+\left(\frac{\partial^{2} u}{\partial x \partial y}\right)_{i, j} \Delta x+\left(\frac{\partial^{3} u}{\partial x^{2} \partial y}\right)_{i, j} \frac{(\Delta x)^{2}}{2}+\left(\frac{\partial^{4} u}{\partial x^{3} \partial y}\right)_{i, j} \frac{(\Delta x)^{3}}{6}+\cdots$$

$$u_{i-1, j}=u_{i, j}+\left(\frac{\partial u}{\partial x}\right)_{i, j}(-\Delta x)+\left(\frac{\partial^{2} u}{\partial x^{2}}\right)_{i, j} \frac{(-\Delta x)^{2}}{2}+\left(\frac{\partial^{3} u}{\partial x^{3}}\right)_{i, j} \frac{(-\Delta x)^{3}}{6}+\cdots$$

$$\left(\frac{\partial u}{\partial y}\right)_{i-1, j}=\left(\frac{\partial u}{\partial y}\right)_{i, j}-\left(\frac{\partial^{2} u}{\partial x \partial y}\right)_{i, j} \Delta x+\left(\frac{\partial^{3} u}{\partial x^{2} \partial y}\right)_{i, j} \frac{(\Delta x)^{2}}{2}+\left(\frac{\partial^{4} u}{\partial x^{3} \partial y}\right)_{i, j} \frac{(\Delta x)^{3}}{6}+\cdots$$

$$\left(\frac{\partial u}{\partial y}\right)_{i+1, j}-\left(\frac{\partial u}{\partial y}\right)_{i-1, j}=2\left(\frac{\partial^{2} u}{\partial x \partial y}\right)_{i, j} \Delta x+\left(\frac{\partial^{4} u}{\partial x^{3} \partial y}\right)_{i, j} \frac{(\Delta x)^{3}}{6}+\cdots$$

$$\to \left(\frac{\partial^{2} u}{\partial x \partial y}\right)_{i, j}=\frac{(\partial u / \partial y)_{i+1, j}-(\partial u / \partial y)_{i-1, j}}{2 \Delta x}-\left(\frac{\partial^{4} u}{\partial x^{3} \partial y}\right)_{i, j} \frac{(\Delta x)^{2}}{12}+\cdots$$

rewite with

$$\left(\frac{\partial u}{\partial y}\right)_{i+1, j}=\frac{u_{i+1, j+1}-u_{i+1, j-1}}{2 \Delta y}+O(\Delta y)^{2} \quad,\quad \left(\frac{\partial u}{\partial y}\right)_{i-1, j}=\frac{u_{i-1, j+1}-u_{i-1, j-1}}{2 \Delta y}+O(\Delta y)^{2}$$

Finally get the result

$$\left(\frac{\partial^{2} u}{\partial x \partial y}\right)_{i, j}=\frac{u_{i+1, j+1}-u_{i+1, j-1}-u_{i-1, j+1}+u_{i-1, j-1}}{4 \Delta x \Delta y}+\mathcal{O}\left[(\Delta x)^{2},(\Delta y)^{2}\right]$$

#### 4. Convergence

$$f(x) − f^{(h)}(x) = Ch^2 + \mathcal{O}(h^{3})$$

convergence test，如果有理论值，我们就可以通过这种最直接的方式

$$\frac{\left|f(x)-f^{(h)}(x)\right|}{\left|f(x)-f^{(h / 2)}(x)\right|}=Q=2^{p}=4$$

self convergence test，大部分情况下解析值无法获得，这时我们也可通过如下形式判断收敛性

$$f^{(h)}(x) − f^{(h/2)}(x) = -\frac{3}{4}Ch^2 + \mathcal{O}(h^{3})$$

$$\frac{\left|f^{(h)}(x)-f^{(h/2)}(x)\right|}{\left|f^{(h/2)} (x)-f^{(h / 4)}(x)\right|}= Q = 2^{p}$$

#### 6.误差与稳定性分析

$$\frac{\phi_{i}^{n+1}-2 \phi_{i}^{n}+\phi _{i}^{n-1}}{\Delta t^{2}}=c^{2} \frac{\phi _{i+ 1}^{n}-2 \phi _{i}^{n}+\phi _{i-1}^{n}}{\Delta x^{2}}$$

• 显式方法：

直接推进求解，一旦$\Delta x$取定，$\Delta t$便不再独立，需要满足稳定性条件，一般要取很小

• 隐式方法：

节点整体联立矩阵求解，没有稳定性限制，对不追求时间精确度的常态解友好

#### 7.CFL condition & Von-Neumann stability

let $A=$偏微分方程精确解 $D$ = 差分方程精确解, $N$ = 有限精度的计算机上实际计算的解

$$\frac{\phi_{i}^{n+1}-2 \phi {i}^{n}+\phi _{i}^{n-1}}{\Delta t^{2}}=c^{2} \frac{\phi _{i+ 1}^{n}-2 \phi _{i}^{n}+\phi _{i-1}^{n}}{\Delta x^{2}}$$

$$\frac{\left(D_{i+1}^{n}+\varepsilon_{i+1}^{n}\right)-2\left(D_{i}^{n}+\varepsilon_{i}^{n}\right)+\left(D_{i-1}^{n}+\varepsilon_{i-1}^{n}\right)}{(\Delta t)^{2}}=c^{2} \frac{\left(D_{i+1}^{n}+\varepsilon_{i+1}^{n}\right)-2\left(D_{i}^{n}+\varepsilon_{i}^{n}\right)+\left(D_{i-1}^{n}+\varepsilon_{i-1}^{n}\right)}{(\Delta x)^{2}}$$

$$\frac{\varepsilon_{i}^{n+1}-2 \varepsilon_{i}^{n}+\varepsilon_{i}^{n-1}}{(\Delta t)^{2}}=c^{2}\, \frac{\varepsilon_{i+1}^{n}-2 \varepsilon_{i}^{n}+\varepsilon_{i-1}^{n}}{(\Delta x)^{2}}$$

$$\left|\frac{\varepsilon_{i}^{n+1}}{\varepsilon_{i}^{n}}\right| \leqslant 1$$

Ansatz: $\quad\varepsilon_{m}(x, t)=\mathrm{e}^{a t} \mathrm{e}^{\mathrm{i} k_{m} x}$

$$\frac{\mathrm{e}^{a(t+\Delta t)} \mathrm{e}^{\mathrm{ik}_{m^{x}}}-2\mathrm{e}^{a t} \mathrm{e}^{\mathrm{ik}_{m^{x}}}+\mathrm{e}^{a(t-\Delta t)} \mathrm{e}^{\mathrm{ik}_{m^{x}}} }{(\Delta t)^{2}}=c^{2}\, \frac{\mathrm{e}^{a t} e^{i k_{m}(x+\Delta x)}-2 \mathrm{e}^{a t} \mathrm{e}^{\mathrm{i} k_{m} x}+\mathrm{e}^{a t} \mathrm{e}^{i k_{m}(x-\Delta x)}}{(\Delta x)^{2}}$$

$$\left|\frac{\varepsilon_{i}^{n+1}}{\varepsilon_{i}^{n}}\right|=\left|\mathrm{e}^{\alpha \Delta t}\right|=\left|1-\frac{2 c^{2} (\Delta t)^{2}}{(\Delta x)^{2}} \sin ^{2} \frac{k_{m} \Delta x}{2}\right| \leqslant 1$$

$$\to \ \alpha=\frac{c \Delta t}{\triangle x} \leq 1$$

same with CFL condition.

CFL存在物理解释！

The solution of the hyperbolic IVP at a given point depends on the information in its domain of dependency; numerical stability is guaranteed if the time-step is sufficiently small such that it contains the domain of dependency for that point.

#### 8. 边界条件

• 周期性是最常见的，线性插值的结果与周期性也类似。

$$\partial _t \phi \pm \partial _x \phi =0$$

Analytical Solution: $\phi (x,t)= \phi(x \pm ct)\ \,$ , with $\ c=\alpha=1$

$$\quad \phi^{n+1}_{i\pm 1}= \phi^{n}_{i}$$

• characteristics analysis: 全反射

$$0 = (\partial _t \mp \partial _x)[(\partial_{t} \pm \partial_{x} )\Phi ] = (\partial _t \mp \partial _x)w_{\pm}$$

General solution

$$w_{\pm}(x,t)=w_{\pm}(x \pm ct)\qquad (+:incoming\ -:outgoing)$$

Since there are no incoming/outgoing waves in the exterior region, it follows $w_{\pm}(L/0,t)=0$

no incoming

$$\left.\left(\partial_{t} + \partial_{x}\right) \Phi(x, t)\right|_{x=L}=0$$

$$\left(\frac{ \phi_{i}^{n+1} - \phi_{i}^{n-1}}{2h_{t}} + \frac{ \phi_{i+1}^{n} - \phi_{i-1}^{n}}{2h_{x}}\right) |_{x=L} =0 \ \to \ \phi_{i}^{n+1} - \phi_{i}^{n} = - \alpha( \phi_{i+1}^{n} - \phi_{i-1}^{n})$$

no outgoing

$$\left.\left(\partial_{t} - \partial_{x}\right) \Phi(x, t)\right|_{x=0}=0$$

$$\left(\frac{ \phi_{i}^{n+1} - \phi_{i}^{n-1}}{2h_{t}} - \frac{ \phi_{i+1}^{n} - \phi_{i-1}^{n}}{2h_{x}}\right) |_{x=0} =0 \ \to \ \phi_{i}^{n+1} - \phi_{i}^{n} = \alpha( \phi_{i+1}^{n} - \phi_{i-1}^{n})$$

$$\phi_{i+1}^{n} = \phi_{i-1}^{n}$$

$$\frac{\partial \phi}{\partial n} = \boldsymbol{n} \cdot \nabla \phi=0$$

centered finite difference stencil

$$\partial _{tt}\phi -c^{2}\partial_{rr}\phi =0$$

$$\frac{\phi_{i}^{n+1}-2 \phi_{i}^{n}+\phi_{i}^{n-1}}{(\Delta t)^{2}} - c^{2} \frac{\phi _{i+1}^{n}-2 \phi _{i}^{n}+\phi _{i-1}^{n}} {(\Delta x)^{2}}=0$$

$$\phi_{i}^{n+1} = 2\phi_{i\pm1}^{n-1} - \phi_{i}^{n-1}$$