Operator Splitting Methods in PDEs

求解偏微分方程时，常用的一个方法，今天才开始好好学习一下，做个笔记和大家分享下。
学习的书籍是：Helge Holden et al. (2010):
Splitting Methods for Partial Differential Equations with Rough Solutions
Analysis and MATLAB programs

为什么要使用算符分解方法求解偏微分方程

• 对于许多问题的模拟，最后往往归结为求解偏微分方程。由于现在计算能力的大大提高，求解的偏微分方程越来越复杂
• 求解过程可能涉及到不同的物理过程，这些过程反映在PDE中，表现为不同的数学项。这些不同的项数学上的表现也大不相同，需要做不同的处理，使得不论进行解析求解还是数值模拟都很困难
• 如果求解的方程是一个新类型的PDE，现有的解算方法可能会失效
• 利用算符分解方法可以将一个复杂的PDE分解为一些简单的PDE，这些简单的PDE一般都有现成的算法求解，最后再将这些简单PDE的解综合起来得到最终结果（类似divide and conquer的思想）

基本思路

• 将PDE方程分解成一系列的子方程，这些子方程是容易求解者有现成的解法
• 使用一定的算符分解算法，将这些子方程的解连起来，形成最终的解

以一种抽象的方式举例，比如求解柯西问题：

$$\frac{\mathrm{d}U}{\mathrm{d}t} + \mathcal{A}(U)=0, ~~ U(0)=U_0$$
其中 $\mathcal{A}$代表一些未指明的算符，以后会举例说明（比如在热流传输过程中的对流项、输运项）

使用分离变量法（我自己的理解，有问题大家请指正），可以获得下列形式的解：

$$U(t) = e^{-t\mathcal{A}}U_0$$

假设算符能够分解，即$\mathcal{A}=\mathcal{A}_1 + \mathcal{A}_2$，并且可以很容易的求解下列子问题：

$$\frac{\mathrm{d}U}{\mathrm{d}t} + \mathcal{A}_j(U)=0, ~~ U(0)=U_0, ~~~ j=1,2$$
获得的解形式为：
$$U(t) = e^{-t\mathcal{A}_j}U_0$$

考虑最简单的形式，假设$t_n=n\Delta t, ~~ n=1,2,...$,$\Delta t >0$，且足够小，则希望获得如下近似：

$$U(t_{n+1})=e^{-\Delta t\mathcal{A}_2}e^{-\Delta t\mathcal{A}_1}U(t_n)$$
假设 $\mathcal{A}_1$和 $\mathcal{A}_2$是对易算符，即满足
$$\frac{\partial \mathcal{A}_2}{\partial U}\mathcal{A}_1 - \frac{\partial \mathcal{A}_1}{\partial U}\mathcal{A}_2 = 0$$
则 $e^{-\Delta t\mathcal{A}_2}e^{-\Delta t\mathcal{A}_1} = e^{-\Delta t\mathcal{A}}$, 获得的 $U(t_{n+1})$是精确解。
但实际上， $\mathcal{A}_1$和 $\mathcal{A}_2$往往是不对易的，此时我们希望
$$U(t_{n+1})=e^{-t\mathcal{A}_j}U_0 = \lim\limits_{\Delta t \to 0,t=n\Delta t} e^{-\Delta t\mathcal{A}_2}e^{-\Delta t\mathcal{A}_1}U(t_n)$$
这个就是著名的Lie–Trotter–Kato公式。将解中的算符 $e^{t\mathcal{A}_j}$替换为数值近似，就获得了求解PDE的数值方法。所有更精确的方法，都是在这个基础上发展得到的。采用不同的数值近似，可以获得不同的方法。

方法的优点

这个方法看上去非常粗糙，但是具有下列优点

• 对子方程，可采用不同的数值或分析方法。在某些情况下，这是非常必要的，可以使求解更简单有效
• 更方便利用一些典型PDE的现有算法，具有更大的灵活性。
• 减少计算过程中对内存的需求，提高求解的稳定性
• 对某些问题，可能是唯一可行的求解方法
• 可以很方便的扩充模拟问题的模型，加入不同的物理过程

方法的缺点

• 如果子系统耦合度非常高，比如所代表的物理过程演化的时间尺度非常短，则算符分解方法可能会失效，分解过程中步长的选取会受到严重限制
• 必须考虑分解造成的误差，以防止出现严重的错误

补充：
对方程：$U_t = U_x$，使用前向欧拉公式，对所有的$\lambda$，解是不稳定的；使用leap-frog方法（2阶时间精度，2阶空间精度），对$\lambda < 1$，解是稳定的
然后，不幸的是，对方程：$U_t = U_{xx}$，使用前向欧拉公式，对$\lambda < 1/2$，解是稳定的；使用leap-frog方法（2x2精度），对所有$\lambda$，解都是不稳定的

现在我们想解下列方程：$U_t = U_x + U_{xx}$，这两种方法还可以用吗？！！

使用算符分解方法，可以做如下处理：
（1）求解$U_t = U_x$，初值为$U_0$，获得解$U_1$（可以使用leap-frog）
（2）求解$U_t = U_{xx}$，用上一步的$U_1$作为初值，获得解$U_2$（可以使用欧拉方法）
（3）将$U_2$当作（1）中的初值$U_0$，循环执行（1），（2），获得下个步长的解，
END

First Order Splitting: Lie-Trotter Splitting

问题仍然如前面所述。
算法如下：

(1)

$$\frac{\partial U(t)} { \partial t} = \mathcal{A}_1 U(t), t [t^n,t^{n+1}], u(t^n)=u^n_{sp}$$
(2)
(3)

First Order Splitting: Additive Splitting

————————————-万恶分割线————————————————-
由于时间和水平所限，错误在所难免，大家批判阅读，特此敬告！