MPDE方法以及Mathematica实现
MPDE方法是一类用来分析数值格式的耗散色散性质的方法，它的思路为： 对于求解问题PDE(a)，我们通过离散设计得到了数值格式， 现在反过来，找出一个与数值格式最契合的PDE(b)，通常是原本的PDE(a)加上一些高阶余项，分析PDE(b)的耗散和色散性质，尤其关注它相比于PDE(a)多出来的高阶项，就可以得到数值格式的耗散和色散性质。
准备
MPDE的分析依赖于下面的结论：对于PDE \[ 0 = u_t - P u = u_t - \sum_{j=1}^n a_j \frac{\partial^j u}{\partial x^j} \] 其中 \(a_i \in \mathbb{R}\) 。一般性的结论如下：
- \(u_x\)：不会对耗散和色散做出任何贡献
- \(u_{xx}\) 或者更高的偶数阶导数项：对耗散性有贡献，对色散性无贡献
- \(u_{xxx}\) 或者更高的奇数阶导数项：对色散性有贡献，对耗散性无贡献
低阶项的贡献是主要的，因此通常可以忽略高阶项，只关注最低阶的系数非零的一个偶数阶导数项和一个奇数阶导数项（不含 \(u_x\) ）。
这里要求系数都是实值，如果系数是复值，相应的结论会更加复杂，需要对实部虚部分开讨论。
MPDE推导（Mathematica实现）
如何从差分格式找到与之最契合的PDE(b)是MPDE方法的核心，也是最繁琐的部分，基本做法就是将精确解在差分格式上展开，得到一个表达式(A)，对(A)不断求导与自身做运算，消去自身除了 \(u_t\) 之外的低价时间导，将时间导转换为空间导，在分析时舍弃高阶项。
MPDE方法的计算量实在是太大了，利用Mathematica的符号计算无疑是更好的选择，编程实现的思路如下：
在程序中按照如下位置将对应系数存储为下三角矩阵 \[ \begin{bmatrix} u \\ u_t & u_x \\ u_{tt} & u_{tx} & u_{xx} \\ u_{ttt} & u_{ttx} & u_{txx} & u_{xxx} \\ u_{tttt} & u_{tttx} & u_{ttxx} & u_{txxx} & u_{xxxx} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix} \] 在实践中只使用 \(n \times n\) 矩阵来存储，对于更高阶系数直接忽略。
例如初始状态的系数矩阵是通过对时间项部分和空间项部分进行泰勒展开得到的 \[ 0 = \left( T_1 u_t + T_2 u_{tt} + T_3 u_{ttt} + \cdots \right) + \left( X_1 u_x + X_2 u_{xx} + X_3 u_{xxx} + \cdots \right) \] 将系数对应存储得到下三角矩阵 \(S_0\) \[ S_0 = \begin{bmatrix} 0 \\ T_1 & X_1 \\ T_2 & & X_2 \\ T_3 & & & X_3 \\ T_4 & & & & X_4 \\ \vdots & & & & & \ddots \\ \end{bmatrix} \] 其中通常 \(T_1=1\)，MPDE的做法是对原始PDE关于时间和空间求导，对应于系数矩阵的变换为：
- 关于时间求导表现为系数矩阵整体向下下移
- 关于空间求导表现为系数矩阵整体向右下角平移
通过逐行从前到后打洞，来消去时间高阶项和时空混合项，最终目标是得到如下形式的系数矩阵（我们只需要最终矩阵的对角线） \[ S = \begin{bmatrix} 0 \\ 1 & X_1' \\ & & X_2' \\ & & & X_3' \\ & & & & X_4' \\ & & & & & \ddots \\ \end{bmatrix} \] MPDE所求的PDE(b)即为 \[ 0 = u_t + X_1' u_x + X_2' u_{xx} + X_3' u_{xxx} + \cdots \]
Mathematica源代码如下 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20Init[n_, tlist_, xlist_] :=
Module[{S, i, j}, S = Table[0, {i, 1, n}, {j, 1, n}];
For[i = 2, i <= n, i++, S[[i, 1]] = tlist[[i - 1]];
S[[i, i]] = xlist[[i - 1]]]; S];
Trans[S_, tn_, xn_, n_] :=
Module[{R, i, j}, R = Table[0, {i, 1, n}, {j, 1, n}];
For[i = 1, i + tn + xn <= n, i++,
For[j = 1, j + xn <= n, j++, R[[i + tn + xn, j + xn]] = S[[i, j]]]];
Return[R]]
Main[S_, n_] :=
Module[{Cur = S, Tmp},
For[i = 3, i <= n, i++,
For[j = 1, j <= i - 1, j++, Tmp = Trans[S, i - j - 1, j - 1, n];
Cur = Cur - Cur[[i, j]]*Tmp]]; Return[Cur]]
MPDE[n_, tlist_, xlist_] :=
Module[{S, R}, S = Init[n, tlist, xlist]; R = Main[S, n];
Return[Diagonal[R] // Factor // Simplify]]
使用时需要设置计算规模 \(n\)：仅保留 \(n\times n\) 的矩阵，还要指定 \(T_i\) 和 \(X_i\) 的具体计算规则，并生成两个足够长的列表tlist和xlist作为参数传入。
对流方程
这一部分考虑对流方程 \[ u_t + \alpha u_x = 0,\quad (\alpha >0) \] 以及对应的数值格式。
例：迎风格式
我们考虑迎风格式（FTBS格式） \[
v_j^{n+1} = v_j^n - \alpha \frac{\Delta t}{\Delta x}(v_j^n - v_j^{n-1})
\] 推导PDE(b)的出发点为如下的泰勒展开式 \[
0 =
\left(
v_t + v_{tt} \frac{\Delta t}2 + v_{ttt} \frac{(\Delta t)^2}6 +
\cdots
\right)
+ \alpha
\left(
v_x - v_{xx} \frac{\Delta x}2 + v_{xxx} \frac{(\Delta x)^2}6 +
\cdots
\right)
\] 因此两组初始系数为 \[
\begin{aligned}
T_m &= \frac{(\Delta t)^{m-1}}{m!}, m = 1,2,\dots\\
X_m &= \alpha \frac{(-\Delta x)^{m-1}}{m!}, m = 1,2,\dots
\end{aligned}
\] 取 \(5 \times 5\)
的系数矩阵计算范围，求解代码为 1
2
3
4
5(*u_t + a u_x = 0 FTBS*)
Tlist1[n_]:= Table[(Deltat) ^(j-1)/j!,{j,1,n}]
Xlist1[n_]:= Alpha * Table[(-Deltax)^(j-1)/j!,{j,1,n}]
n=5;MPDE[n,Tlist1[n],Xlist1[n]]
计算结果的对角线为（记 \(r = \alpha \frac{\Delta t}{\Delta x}\) ） \[ 0,\,\, \alpha ,\,\, \frac{\alpha \Delta x}2 (r-1), \frac{\alpha (\Delta x)^2}6(r-1)(2r-1),\,\, \frac{\alpha (\Delta x)^3}{24}(r-1)(6r^2-6r+1) \]
与它最契合的PDE(b)为（记 \(r = \alpha \frac{\Delta t}{\Delta x}\) ） \[ \begin{aligned} 0 ={}& u_t + \alpha u_x + \frac{\alpha \Delta x}2 (r-1) u_{xx} + \frac{\alpha \Delta x^2}6 (r-1)(2r-1) u_{xxx} \\ &+ \frac{\alpha (\Delta x)^3}{24}(r-1)(6r^2-6r+1) u_{xxxx} + \cdots \end{aligned} \] 我们只需要关注两个主项：\(u_{xx}\) 和 \(u_{xxx}\) 的系数。
对于差分格式自身的分析： 由于 \(u_{xx}\) 的系数为负，差分格式存在耗散性。 关注 \(u_{xxx}\) 的系数，由于PDE(b)的波速为 \[ c(k) = \alpha - \frac{\alpha \Delta x^2}6 (r-1)(2r-1) k^2 \] 因此差分格式自身具有色散性。
差分格式与PDE的性质比较： 易知PDE(a)自身无色散无耗散，因此始终具有数值耗散。 PDE(a)的波速为 \(c(k) = \alpha\)，与PDE(b)的波速进行比较可得：
- \(r \in (\frac12,1)\)，差分格式具有数值正色散
- \(r \in (0,\frac12)\)，差分格式具有数值负色散
扩散方程
这一部分考虑扩撒方程 \[ u_t = \beta u_{xx},\quad (\beta >0) \] 以及对应的数值格式。
例：FTCS格式
我们考虑FTCS格式：（记 \(r = \beta
\frac{\Delta t}{(\Delta x)^2}\) ） \[
v_j^{n+1} = (1-2r) v_j^n + r(v_{j+1}^n + v_{j-1}^n)
\] 推导PDE(b)的出发点为如下的泰勒展开式 \[
0 =
\left(
v_t + v_{tt} \frac{\Delta t}2 + v_{ttt} \frac{(\Delta t)^2}6 +
\cdots
\right)
- \beta
\left(
v_{xx}
+ v_{xxxx} \frac{(\Delta x)^2}{12}
+ v_{xxxxxx} \frac{(\Delta x)^4}{360} + \cdots
\right)
\] 因此两组初始系数为 \[
\begin{aligned}
T_m &= \frac{(\Delta t)^{m-1}}{m!}, m = 1,2,\dots\\
X_m &= - \beta \left[1+(-1)^m\right] \frac{(\Delta x)^{m-2}}{m!}, m
= 1,2,\dots
\end{aligned}
\] 取 \(6 \times 6\)
的系数矩阵计算范围，求解代码为 1
2
3
4
5(*u_t = b u_xx FTCS*)
Tlist2[n_]:= Table[(Deltat)^(j-1)/j!,{j,1,n}]
Xlist2[n_]:= -Beta * Table[(1+(-1)^j) (Deltax)^(j-2)/j!,{j,1,n}]
n=6; MPDE[n,Tlist2[n],Xlist2[n]]
计算结果的对角线为（记 \(r = \beta \frac{\Delta t}{\Delta x^2}\) ） \[ 0,\,\, 0,\,\, -\beta,\,\, 0,\,\, \frac{\beta (\Delta x)^2}{2}(r - \frac{1}6),0 \]
与它最契合的PDE(b)为（记 \(r = \beta \frac{\Delta t}{\Delta x^2}\) ） \[ 0 = u_t - \beta u_{xx} + \frac{\beta (\Delta x)^2}{2}(r - \frac{1}6) u_{xxxx} + \cdots \]
对于差分格式自身的分析： 由于 \(u_{xx}\) 的系数为负，差分格式存在耗散性。 由于不存在奇数阶项，因此差分格式自身没有色散性。
差分格式与PDE的性质比较： 易知PDE(a)自身有耗散无色散，PDE(b)的耗散和PDE(a)相比，还有 \(u_{xxxx}\) 的额外贡献：
- \(r \in (\frac16,\frac12)\)，差分格式具有数值正耗散
- \(r \in (0,\frac16)\)，差分格式具有数值逆耗散
PDE(a)和PDE(b)均没有色散性，因此差分格式没有数值色散。