色散关系
考虑如下两个方程：
$$
\begin{cases}
u_t=u_{xx} \\
u_t=u_{xxx}
\end{cases}
$$
取试探解 $u=e^{\sigma t+ikx}=e^{i(\omega t+kx)}$，分别代入可得两个关系：
$$
\begin{cases}
\sigma_1=(ik_1)^2=-k_1^2 \\
\sigma_2=(ik_2)^3=-ik_2^3
\end{cases}
$$
可以看到，$\sigma_1 \in \mathbb{R}, \sigma_2 \in \mathbb{C}$，前者为实数，后者为虚数，因此：
$$
\begin{cases}
u_1=e^{-k_1^2 t}e^{ik_1x} \\
u_2 = e^{i(-k_2^3t+k_2x)}
\end{cases}
$$
所以 $u_1$ 的振幅随时间衰减，$u_2$ 的相位随时间改变。
更一般的，偶数阶方程 $u_t=u_{2nx}$ 会表现为振幅随时间改变，可以看作行波扩散；
奇数阶方程 $u_t=u_{(2n+1)x}$ 会表现为相位随时间改变，可以看作行波色散。
这里的色散是指色散关系：
$$
\begin{cases}
v_{phase} = \frac{\omega}{k}=\frac{\sigma}{ik} \\
v_{group} = \dv{\omega}{k}=\frac{1}{i}\dv{(\sigma)}{k}
\end{cases}
$$
对于奇数阶方程 $u_t=u_{(2n+1)x},n \ge 0$ 来说，有：
$$
\sigma=(ik)^{2n+1}=i(-1)^nk^{2n+1}
$$
所以：
$$
\begin{cases}
v_{phase} = \frac{\sigma}{ik}=(-1)^nk^{2n} \\
v_{group} = \dv{\omega}{k}
= \frac{1}{i}\dv{(i(-1)^nk^{2n+1})}{k}
=(-1)^{n}(2n+1)k^{2n}
\end{cases}
$$
当波包的群速度和相速度不相等时，即 $v_{phase} \ne v_{group}$ 时，会发生色散（不同行波的速度不同）。
特别的，当 $n=0$ 时，波动方程 $\pdv{u}{t}+a\pdv{u}{x}=0$ 是不会发生色散的。
可以看到，数值耗散和数值色散是PDE的固有性质。
如果我们对一个PDE进行有限差分，那么由于截断误差的存在，方程右端会出现关于 $x$ 的导数，
而导数的阶数往往会同时存在偶数阶和奇数阶（甚至是一个无穷级数）。
因此，数值解必然会存在数值耗散和数值色散的问题，
这是截断误差所引发的。
扩散方程
考虑扩散方程：
$$
u_{t}=au_{xx}, a>0
$$
一般而言，我们要求扩散系数 $a>0$，
一方面这是符合物理意义（物质从浓度高的地方向浓度低的地方进行扩散，即符合统计意义下的物质的无规则运动）的，
另一方面如果 $a<0$ 会导致PDE的不稳定。
考虑一个波 $e^{\sigma t+ikx}$，代入方程得到：
$$
\sigma = -ak^2
$$
如果 $\sigma>0$，会导致振幅随时间不断增大，最终导致系统不稳定。
从物理意义上来看，如果扩散系数 $a>0$，说明物质沿着浓度梯度从高浓度向低浓度扩散，这抑制了物质的不断积累，
一个区域内的物质浓度越高，那么它向周围扩散的趋势越大，最终减少了该域的物质；
但反过来，如果扩散系数 $a<0$，说明物质会从浓度低的区域向浓度高的区域进行汇聚，
这导致的结果是整个系统由多个 $\delta$ 函数组成，最终形成若干聚点。
波动方程
考虑波动方程：
$$
u_t+au_{x}=0,a>0
$$
对于空间项，有如下几种离散方法：
$$
\begin{aligned}
& \pdv{u}{x}\cong \frac{\Delta}{h}u, \quad \Delta=E-1 \\
& \pdv{u}{x}\cong \frac{\nabla}{h}u, \quad \nabla=1-E^{-1} \\
& \pdv{u}{x}\cong \frac{\Delta_0}{2h}u, \quad \Delta_0=E^{1}-E^{-1}
\end{aligned}
$$
可以看到：
$$
\begin{aligned}
\Delta=&E-1=e^{hD}-1=hD+h^2D^2/2+O(h^3) \\
\nabla =& 1-E^{-1}=1-e^{-hD}=hD-h^2D^2/2+O(h^2) \\
\Delta_0=&E^{1}-E^{-1}
=e^{hD}-e^{-hD} \\
=&\qty(1+hD/1+h^2D^2/2+h^3D^3/6+O(h^4)) \\
& \quad -\qty(1-hD/1+h^2D^2/2-h^3D^3/6+O(h^4)) \\
=&2hD+h^3D^3/3+O(h^4)
\end{aligned}
$$
所以：
$$
\begin{aligned}
& \frac{\Delta}{h}u = \dv{u}{x}+\frac{h}{2}\dv[2]{u}{x}+O(h^2) \\
& \frac{\nabla}{h}u = \dv{u}{x}-\frac{h }{2}\dv[2]{u}{x}+O(h^2) \\
& \frac{\Delta_0}{2h}u = \dv{u}{x}+\frac{h^2}{6}\dv[3]{u}{x}+O(h^3)
\end{aligned}
$$
对于时间项，一般的离散方法为：
$$
\frac{\Delta_t}{\tau}u=\frac{E_\tau-1}{\tau}u
$$
注意到 $u_t=-au_x, D_t=-aD_x$，所以：
$$
\begin{aligned}
& E_\tau-1=e^{\tau D_t}-1=\tau D_t + \tau^2D_t^2/2 +\tau^3D_t^3/6 + O(\tau^4) \\
\Rightarrow
& \frac{E_\tau-1}{\tau}
=D_t+\frac{\tau a^2}{2}D_x^2-\frac{\tau^2 a^3}{6}D_x^3+O(\tau^3)
\end{aligned}
$$
所以：
$$
u_t+au_{x}
=
\underbrace{\qty(\frac{E_\tau-1}{\tau}+a\widetilde{\partial_x})}_{\text{Finite Difference, =0}}
-\begin{cases}
\qty(\frac{\tau a^2}{2} \pm \frac{ah}{2})D_x^2 - \frac{\tau^2a^3}{6}D_x^3+O(\tau^3,h^2),
\quad& \Delta,\nabla\\
\frac{\tau a^2}{2} D_x^2 + \frac{ah^2-\tau^2 a^3}{6}D_x^3+O(\tau^3,h^3), \quad& \Delta_0
\end{cases}
$$
定义网格比 $\lambda=\frac{a\tau}{h}$。
对于中心差分格式，总有扩散系数：
$$
-\frac{\tau a^2}{2} < 0
$$
因此中心差分格式是无条件不稳定的。
对于迎风格式，扩散系数为：
$$
-\frac{\tau a^2}{2} \mp \frac{ah}{2}
=\frac{ah}{2}\qty(-\lambda \mp 1)
$$
令扩散系数大于 $0$，则：
- $a>0(\lambda>0)$，则使用 $\nabla$，要求 $\lambda<1$；
- 若 $a<0(\lambda<0)$，则使用 $\Delta$，要求 $-\lambda<1$。
实际上，$D_t,D_x$ 的近似中，后面还有更多的余项：
$$
\begin{aligned}
& \frac{E_\tau-1}{\tau}
=D_t+\sum_{n=1}^{\infty} \frac{\tau^{n-1}a^n}{n!}D_x^n \\
& \frac{E-1}{h}=D_x+\sum_{n=1}^{\infty} \frac{h^{n-1}}{n!}D_x^n
\end{aligned}
$$
对于偶数项，要求系数小于 $0$，即：
$$
\begin{aligned}
& ah^{n-1} + \tau^{n-1} a^n < 0 \\
\Rightarrow
& a\qty(1+\qty(\frac{a\tau}{h})^{n-1}) < 0
\end{aligned}
$$
显然必有 $a<0,-\lambda<1$。
考虑使用二阶迎风格式，对波动方程分析：
$$
u_t+vu_x=0
$$
为了得到二阶迎风格式，考虑到 $E=e^{hD}=1+\Delta$，需要将 $D_t$ 展开两次：
$$
\begin{aligned}
hD
&=\ln(1+\Delta) \\
&=\Delta-\frac{\Delta^2}{2} \\
&=E-1-\frac{(E-1)^2}{2} \\
&=\frac{2E-2-E^2-1+2E}{2} \\
&=\frac{4E-E^2-3}{2} \\
\end{aligned}
$$
所以 $D_x=\frac{4E_h-E_h^2-3}{2h}$。
考虑波动方程中，有 $D_t=-vD_x$，对时间项使用显格式，有：
$$
\begin{aligned}
\frac{E_{\tau}-1}{\tau}
&=\frac{e^{\tau D_t}-1}{\tau} \\
&=D_t+\sum_{n=2}^{\infty}\frac{\tau^{n-1}}{n!}D_t^{n} \\
&=D_t+\sum_{n=2}^{\infty}\frac{\tau^{n-1}(-v)^{n}}{n!}D_x^{n} \\
&=D_t+\frac{\tau v^2}{2}D_x^2+\sum_{n=3}^{\infty}\frac{\tau^{n-1}(-v)^{n}}{n!}D_x^{n} \\
\end{aligned}
$$
对于空间项，差分格式对应的截断误差为：
$$
\begin{aligned}
\frac{4E_h-E_h^2-3}{2h}
&=\frac{4e^{hD_x}-e^{2hD_x}-3}{2h} \\
&=\frac{1}{2h}\qty(\sum_{n=0}^{\infty}\frac{(4-2^n)h^n}{n!}D_x^n-3) \\
&=\frac{1}{2h}\qty(2hD+\sum_{n=3}^{\infty}\frac{(4-2^n)h^n}{n!}D_x^n) \\
&=D_x+\sum_{n=3}^{\infty} \frac{2(1-2^{n-2})h^{n-1}}{n!}D_x^n
\end{aligned}
$$
可以看到，该差分格式存在扩散系数为 $-\frac{\tau v^2}{2}$ 的二阶扩散项，
这是一个反常扩散，因此该格式是不稳定的。
可以看到，空间项的二阶精度导致了系统的不稳定，
即空间项的截断误差不存在 $D_x^2$ 与时间项的截断误差产生的 $D_x^2$ 的系数进行整合。
考虑时间项进行二阶差分：
$$
\begin{aligned}
D_t
&=\frac{1}{\tau}\ln(\frac{1}{1-\nabla}) \\
&=\frac{-1}{\tau}\ln(1-\nabla) \\
&=\frac{-1}{\tau}\qty(-\nabla-\frac{\nabla^2}{2}) \\
&=\frac{2-2E^{-1}+\qty(1-E^{-1})^2}{2\tau} \\
&=\frac{E^{-2}+3-4E^{-1}}{2\tau}
\end{aligned}
$$
截断误差为：
$$
\begin{aligned}
\frac{E_\tau^{-2}+3-4E_\tau^{-1}}{2\tau}
&=\frac{1}{2\tau}\qty(e^{-2\tau D_t}+3-4e^{-\tau D_t}) \\
&=D_t+\frac{1}{2\tau}\qty(\sum_{n=2}^{\infty} \frac{(-2\tau)^n-4(-\tau)^n}{n!}D_t^n) \\
&=D_t+\sum_{n=2}^{\infty}\frac{\tau^{n-1}\qty(-(-2)^{n-1}-2(-1)^n)}{n!}(-v)^nD_x^n \\
&=D_t+\sum_{n=2}^{\infty}\frac{2v^n\tau^{n-1}\qty(2^{n-2}-1)}{n!}D_x^n \\
&=D_t+\sum_{n=3}^{\infty}\frac{2v^n\tau^{n-1}\qty(2^{n-2}-1)}{n!}D_x^n \\
\end{aligned}
$$
考虑扩散主项的系数，其为：
$$
\begin{aligned}
\sigma
&=-\frac{6v^4\tau^3}{24}-v\frac{2(1-4)h^3}{24} \\
&=\frac{vh^3-v^4\tau^3}{4} \\
&=\frac{vh^3}{4}\qty(1-\qty(\frac{v\tau}{h})^3)
\end{aligned}
$$
但需要注意的是，我们使用 $\nabla_\tau$ 做前向差分后，由于对流项的时间偏移为 $E^0$，
所以这种差分格式为隐格式。
我们考虑采用 $E^{-1},E^{0},E^{1}$ 来构造格式，这样可以得到二阶精度显格式的构造。
可以采用时间项的中心差分：
$$
\frac{\Delta_{0,\tau}}{\tau}=\frac{E_\tau^{1}-E_\tau^{-1}}{2\tau}
$$
则有：
$$
\begin{aligned}
\frac{\frac{E_\tau-1}{\tau}+\frac{1-E_\tau^{-1}}{\tau}}{2}
&=\frac{E_\tau-E_\tau^{-1}}{2\tau} \\
&=\frac{e^{\tau D_\tau}-e^{-\tau D_\tau}}{2\tau} \\
&=D_\tau+\frac{\sum_{n=2}^{\infty}\frac{\tau^n(1-(-1)^n)}{n!}D_\tau^n}{2\tau} \\
&=D_\tau+\frac{\sum_{n=3}^{\infty}\frac{\tau^n(1-(-1)^n)}{n!}D_\tau^n}{2\tau} \\
&=D_\tau+\sum_{n=3}^{\infty}\frac{\tau^{n-1}(1-(-1)^n)}{2n!}D_\tau^n
\end{aligned}
$$
可以看到，中心差分不具有扩散项，所以中心差分格式往往不具有稳定性。
考虑如下双重中心差分格式：
$$
\begin{aligned}
\frac{E_\tau-1+E_\tau^{-1}-E_\tau^{-2}}{2\tau}
&=\frac{1}{2\tau}\sum_{n=1}^{\infty} \frac{\tau^n+(-\tau)^{n}-(-2\tau)^{n}}{n!}D_t^n \\
&=\frac{1}{2\tau}\sum_{n=1}^{\infty} \frac{\tau^n\qty(1+(-1)^{n}-(-2)^{n})}{n!}(-vD_x)^n \\
&=D_t+\frac{-v^2\tau}{2}D_x^2+\frac{1}{2\tau}\sum_{n=3}^{\infty} \frac{(-\tau v)^n\qty(1+(-1)^{n}-(-2)^{n})}{n!}D_t^n
\end{aligned}
$$
