Operator Splitting

math

---

Problem

$$\frac{du}{dt} + A(u,t) = 0, \quad u(0) = u_0.$$ The function $u$ can further dependen on the spatial variable $x$ and $A(u)$ could a be linear/nonlinear operator involving partial differential operator in space, i.e. (1) could be a PDE not just ODE. Indeed the operator splitting method is mostly used as PDE not ODE solvers.

Forward vs Backward Euler Methods

Going from $t_n$ to $t_{n+1} = t_n + \delta t$, we can use either

• Forward Euler: explicit method.
  
  • Pro: simple matrix-vector product
  • Con: stability constraint requiring small $\delta t$
• Backward Euler: implicit method
  
  • Pro: unconditional stable
  • Con: solve a non-trivial equation at each step

Basic Idea of Operator Splitting Methods

Write $A = A_1 + A_2$ with ``simpler" operators $A_1$ and $A_2$. Then evolve them using implicit and explicit scheme alternatively.

The Peaceman-Rachford method
From $t_n$ to $t_{n+1} := t_n + \delta t$, we set $t_{n+1/2} = t_n + \delta t/2$.

• $t_n \to t_{n+1/2}$: implicit $A_1$ and explicit $A_2$
• $t_{n+1/2} \to t_{n+1}$: explicit $A_1$ and implicit $A_2$

The Douglas-Rachford method

• $t_n \to t_{n+1}$: implicit $A_1$ and explicit $A_2$
• $t_n \to t_{n+1}$: explicit $A_1$ and implicit $A_2$

Write out detailed algorithm.

Example: C-N and Euler Methods
$A_1=A, A_2=0$ or $A_1=0, A_2=0$
- P-R will give two versions of C-N
- D-R will give two versions of Euler's method

Example: ADI (Alternating Direction Implicit) Method
$u_t = \Delta u$.

$$A = -\partial_{xx} - \partial_{yy}, A_1 = -\partial_{xx}, A_2 = - \partial_{yy}.$$
Each step solve a block tri-diagonal matrix which is much easier than solving $A$.

Stability

Consider the simple linear problem

$$\frac{du}{dt} + Au = 0, \quad u(0) = u_0,$$ with a SPD operator $A$. After a time discretization, we could write generitically $u^{n+1} = R(\delta t,A)u^n$.

Stablity: for some norm $\|\cdot \|$, we have

$$\|u^{n+1}\|\leq \|u^n\|.$$ With such stability, we could obtain the error estimate
$$\|U^n - u_I^n\| \leq \delta t \sum_i \tau_i(\delta t, h),$$ where $\tau_i(\delta t, h)$ is the truncation error.

The stability is equivalent to study the norm $\|R(\delta t,A)\|$ for some approriate norms (e.g. $L^2$, $L^{\infty}$ or $A$-norm). As $A$ is SPD and $R$ is a rational polynomail, it is reduced to study the function $R(\delta t,\lambda) = R(\delta t \lambda) = R(\xi)$, where $\lambda>0$ is an eigenvalue of $A$ and $\xi = \delta t \lambda >0$.

Consider the inequality $|R(\xi)|\leq 1$. If it holds for all $\xi>0$, we call it unconditionally stable. Otherwise the upper bound $\xi < M$ put a restriction on the step size $\delta t$ which is usually in the form $\delta t = \mathcal O(\lambda_{\max}^{-1})$.

For general operator $A$, we could extend the domain of $R(\xi)$ from the real line to the complex plane and consider range of $Re(\xi)$, the real part of $\delta t \lambda$.

Remark In ODE solvers, the eqn is in the form $du/dt = Au$ which is a sign difference with our setting. For example, the unconditional stable means $|R(\xi)|\leq 1$ for all $\xi$ with $Re(\xi)<0$ not $Re(\xi)>0$ in our setting.

Stronger Stability

The exact solution of (*) satisfy $u(t) = e^{-At}u_0 = \sum_{i} c_ie^{-\lambda_i t} \to 0$ no matter what the initial condition is. After a time discretization, we would like to preserve the decay property. To do so, we need a stronger stability (also known as A-stable or L-stable) result

$$|R(\xi)|<1, \text{ and } \lim_{\xi\to \infty}|R(\xi)|< 1.$$

It might happen $\lim_{\xi\to 0}|R(\xi)| = 1$ which is acceptable. As $\xi \ll 1$ corresponds to small eigenvalues for which the component decay with a smaller rate $e^{-\lambda t} \approx 1 - \xi$ (the graph is relative flat) and thus the contraction rate $|R(\xi)|$ close to $1$ is acceptable for $\xi \ll 1$, i.e. for small eigenvalues of $A$ (fast transition vs slow transition). For large eigenvalues of $A$, the exact decay rate $e^{-\lambda t}$ is large. Then $|R(\xi)|$ close $1$ is not accurate.

Consider the case of $|R(\xi)|<1$ but $\lim_{\xi\to \infty}|R(\xi)|= 1$. Then to capture the decay property in the discretization, $\delta t \lambda_{\max} \leq M$ is still needed. When $\lambda_{\max}(A)$ is large, a small step size is needed to have the correct decay property even the solution is smooth after a fast transition. This difficulty is known as stiffness.

Consider

$$\frac{du}{dt} + Au = f, \quad u(0) = u_0.$$ The exact solution $$u(x,t) = \sum c_ie^{-\lambda_i t} + g(x),$$ where $g(x)$ is the so-called steady-state solution (independent of $t$). The solution $u(x,t)$ will converge to $g(x)$ as $t\to \infty$.

Stronger stable scheme will be also more suitable for computing the steady-state solution.

Crank-Nicolson method is unconditional stable and second order time discretization but it is not strong stability. The polynomial $R(\xi) = (1- \xi/2)/(1+\xi/2) \to 1$ as $\xi \to +\infty$.

Backward Euler is strong stable as $R(\xi) = 1/(1+\xi)$ but only first order.

A second order and strong stable scheme is BDF_2: suppose $u^n, u^{n-1}$ are known, we obtain $u^{n+1}$ by solving

$$\left (\frac{3}{2}u^{n+1} - 2 u^n + \frac{1}{2} u^{n-1}\right )/\delta t + A(u^{n+1}, t_{n+1}) = 0.$$

Accuracy

Formally this is done by comparing Taylor series. First

$$e^{-\delta t} = 1 - \delta t + \frac{1}{2}\delta t^2 + \frac{1}{6}\delta t^3 + ...$$
Then we expand $R(\delta t)$ at $\delta t = 0$. If $|R(\delta t) - e^{-\delta t}| = \mathcal O(\delta t^{k+1})$, then we say the scheme of order $k$. The power is reduced by one since we consider the finite difference approximation to partitial derive of $t$ and cancel one $\delta t$ in the denominator.

Example

• Euler
• C-N
• BDF_2

A $\theta$-method

Now we set two inter middle points in one time intervale. Given $\theta\in (0,1/2)$, we set $t_{n+\theta}$ and $t_{n+1-\theta}$. Again we split $A=A_1+A_2$ and iterative:

• $t_n \to t_{n+\theta}$: implicit $A_1$, explicit $A_2$
• $t_{n+\theta}\to t_{n+1-\theta}$: implicit $A_2$, explicit $A_1$
• $t_{n+1-\theta}\to t_{n+1}$: implicit $A_1$, explicit $A_2$

What is the optimal choice of $\theta$? How do we split $A$? We study a simplified problem and splitting scheme: $A(u,t) = Au$ is linear and $A_1=\alpha A, A_2=\beta A$.

Write out $R(\xi)$ and study the stability and accuracy.