MENU

波动方程的数值求解

December 4, 2021 • 铃宕阅读设置

前记:这其实是我不久前做的人生第一个数值求解偏微分方程组的练习,感觉老师布置的内容实在是太好了,所以整理成笔记记录下来。本篇可以视为数值求解PDE的入门概略。

本文以1维的平面波为例,但均可以直接拓展到3维上

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

1.a 数学求解:hyperbolic PDE,解法无聊但适用范围广

所有二阶线性双变量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 = B^{2}-AC $的取值,可以分成三类, 并均可以化成以下标准形式

  • $\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}} = [...] $$

其中,$[...]$代表所有不含二阶偏导的项。所以波动方程属于双曲型PDE

下面直接套用数学上的结论,可以用特征方程的特征线进行分解,从而求解波动方程。

$$ 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$用两个特征线改写,$w(\xi, \eta)=\phi(x(\xi, \eta), t(\xi, \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 $$

因为 $\left(w_{\xi}\right)_{\eta}=0$, 可以进而假设 $w_{\xi }=f(\xi )$, 然后因为$w=\int f(\xi) \mathrm{d} \xi+G(\eta)$,最终得到

$$ 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 $$

同样引入新变量, $\xi=t-\frac{x}{c}, \quad \eta=t+\frac{x}{c}$

显然,$t=\frac{1}{2}(\eta+\xi), \quad x=\frac{c}{2}(\eta-\xi)$

进一步,不难证明,

$$ \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) $$

现在我们对解的物理意义进行理解,里面随便抽一个部分都是一个单向波,

以G为例子,也就是说对于$x- ct = const$的坐标$x$以及时间$t$

$$ x = const + ct $$

这样的时空坐标下,场就有相同的值,即

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

也就是说$t=0$时刻的$x$点的场$\phi$,会在$t$时间后传播到$x-ct$的位置上。

也就是说无论前向波G还是后向波F,这些电磁波都以光速在空间中沿某一方向传播。

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,t)=F(x+ct)+G(x-ct)$, 我们现在需要利用初始的$\theta\ ,\ \psi$函数求解$F,$

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

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

将上式(22)在 $[0, x]$ 区间积分

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

其中 $C = F(0)-G(0)$.

我们便得到了线性方程组的通解

$$ \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} $$

简单情况一: $\psi = 0$

$$ \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)] $$

简单情况二:$\theta = 0$

$$ \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)] $$

其中 $g(x)=\frac{1}{c}\int^{x}_{x_0} \psi(s) \mathrm{d}s$

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

3. Finite differencing approximation 有限差分近似

所有的数值求解都需要首先离散化,也就是说我们需要将时空坐标网格化,

假设离散后的最小尺度是$h$,那么通常误差项也在$h$量级

但利用泰勒展开,我们可以将任一函数$u$的导数和二阶导数逼近为更低阶的误差项,比如$h^2,h^4$

$$ 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} $$

5. Method of line

利用这种方式,我们可以将偏微分方程直接转换为常微分方程进行求解。比如用RK4进行求解

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}} $$

注意这里差分方程是一个代数方程,它与原方程并不是一个东西。

我们希望网格点趋于无穷多时,差分方程能够还原成原微分方程

也就是说,如果离散误差也在这种情况下趋近为0,我们就认为差分方程和原方程是相容的

  • 显式方法:

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

  • 隐式方法:

    整体联立矩阵求解,没有稳定性限制

7.CFL condition & Von-Neumann stability

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

离散误差:$=A-D$

截断误差:$\text{Round-off}=\varepsilon =N-D$

$$ \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 $$

也就是说在t之后的推进中,误差不会增加,数值解呈现稳定态势,不然误差会持续增大直至爆掉。

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. 边界条件

  • 周期性是最常见的,线性插值的结果与周期性也类似。
  • Advection condition:Impose advection equation on the boundary

$$ \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 $$

但这里我们如何处理Ghost点呢?答案是联立差分方程消掉

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 $$

简化起见,set $\Delta t = \Delta x, c=1$

全反射的情况下,

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

或者更简单的情况是直接把边界设为0,其解是一样的。

Last Modified: February 18, 2024