Splitting and ADI. Nonlinear Problems and Problems in Several Space Dimensions

Splitting techniques can be used to efficiently solve certain nonlinear evolution problems and equations in several space dimensions. Consider the following example:

$$\frac{\partial U}{\partial t}=\nabla (a\nabla U)+f\quad 0\le x, y\le 1$$

$$\hfill t>0$$

$$a=a(x,y) \mbox { is bounded in }[0,1]\times [0,1]$$

For simplicity, assume that the grid spacing in $x$ and $y$ is the same: define such grid spacing by $\Delta x$. A naive discretization in space leads to

$$u'_{k,\ell}=\frac{1}{(\Delta x)^2}\left[a_{k-\frac12,\ell}u_{k-1,... ..._{k,\ell-1}+a_{k+\frac12,\ell}u_{k+1,\ell} a_{k,\ell+\frac12}u_{k\ell+1}\right.$$

$$\left. -\left(a_{k-\frac12,\ell}+a_{k,\ell-\frac12}+a_{k+\frac12,\ell}u+ a_{k,\ell+\frac12}\right)u_{k,\ell}\right]+p_{k,\ell}+f_{k,\ell}$$

$$\mbox {with } k, \ell=1\ldots d$$

$$\mbox {and } u'\equiv \frac{\partial }{\partial t} u$$

with $k,\ell=1\ldots d$

$p$ and $f$ are the boundary and forcing term contributions, assumed to be dependent of time. Assume they are 0 for now. In this case the above equation is of the form

$${\bf u}'=\frac{1}{\Delta x}(B_x+B_y){\bf u}\quad t\ge 0, \quad {\bf u}(0) \mbox { given}$$ (153)

$B_x$ and $B_y$ are $d^2\times d^2$ matrices approximating differential operators in $x$ and $y$.

$B_y$ is a block-diagonal matrix and its diagonal is constructed from the tri-diagonal $d\times d$ matrices:

$$\left[\begin{array}{cccccc} -\left(b_{\frac12}+b_{\frac32}\right)... ...0 &b_{d-\frac12} &-\left(b_{d-\frac12}+b_{d+\frac12}\right) \end{array}\right]$$

where $b_{\ell}=a_{k,\ell}$     $k=1,2,\cdots d$

$B_x$ contains all remaining terms.

Its sparsity pattern is block-tridiagonal provided the grid is ordered by rows rather than by columns. Any other ordering will lead to a sparse but large bandwidth system.

The formal solution of (154) is

$$u^{n+1}=e^{\mu(B_x+B_y)}u^n\quad n\ge 0$$

Now we use a Padé approximant to discretize the matrix operators: for example using $\hat r_{1/1}$ to approximate the exponential operator, leads to

$${\bf u}^{n+1}=\Big[I-\frac12\mu(B_x+B_y)\Big]^{-1} \Big[I+\frac12\mu (B_x+B_y)\Big]{\bf u}^n$$ (154)

$$(Crank Nicholson 2D)$$

The following identity is true for matrices:

$$e^{t(Q+S)}=e^{tQ}e^{tS}\quad t\ge 0$$ (155)

ONLY IF $Q$ and $S$ COMMUTE. $Q$ and $S$ are square and of equal dimension. Assume they do, then:

$${\bf u}^{n+1} = \hat r_{1/1}(\mu B_x)\hat r_{1/1}(\mu B_y){\bf u}^n$$ (156)

$$= \underbrace {\left(I-\frac12\mu B_x\right)^{-1}}_I \;(I+\frac12\mu B_x)(I-\frac12\mu B_y) \underbrace{\left(I+\frac12\mu B_y\right)}_{II}{\bf u^n}$$

here $n\ge 0$

The advantage of solving (157) over (155) is that inspite of having to solve 2 linear systems at each time step, with $I-\displaystyle \frac12\mu B_y$ tridiagonal, and $I-\frac12\mu B_x$ tridiagonal is that it can solve (157) in $O(d^2)$ operations!

Identity (156) would be true in this case if $a=1$. In general, the identity is only approximately true:

$$e^{tQ}e^{tS}-e^{t(Q+S)}=(I+tQ+\frac12t^2Q^2+\cdots)(I+tS+\frac12t^2S^2+ \cdots)$$

$$-\left\{I+t(Q+S)+\frac{t^2}{2}(Q+S)^2+\cdots\right\}=\frac12t^2[Q,S]+O (t^3)$$

$$\mbox {where } [Q,S]\equiv QS-SQ \quad \mbox{the cummutator}$$

$$\therefore \mbox {if }[B_x, B_y]\ne 0$$

$$\mbox {then } e^{\mu B_x}e^{\mu B_y}-e^{u(Bx+By)}=0(\mu^2)$$

if $\mu$ is sufficiently small, the approximation is fruitful.

Exercise

Show that ``Strang Splitting''

$$e^{\mu(B_x+B_y)}\approx e^{\frac12\mu B_x}e^{\mu B_y}e^{\frac12 \mu Bx} +O(\mu^3)$$

$\Box$

What to do if boundary conditions and/or forcing non-zero?

Take

$${\bf u}'=\frac{1}{(\Delta x)^2}(B_x+B_y){\bf u}+{\bf h}(t),\; t\ge 0$$

$${\bf u}(0) \mbox { given}$$

then

${\bf u}^{n+1}=e^{\mu(B_x+B_y)}{\bf u}^n+\Delta t\int^{1}_0e^{(1-\tau)\mu(B_x+B_y)} {\bf h}\left((n+\tau)\Delta t\right)d\tau$

$n\ge 0$, and replace integral by trapezoidal rule:

${\bf u}^{n+1}=e^{\mu(Bx+By)}\left[{\bf u}^{n}+\frac12\Delta t {\bf h}(n\Delta t)\right]+\frac12\Delta t {\bf h} \left((n+1)\Delta t\right)$

Example Can use the split-step to efficiently solve the Nonlinear Schrodinger Equation (NLS). First, we use $\hat r_{1/1}$ since it is a unitary approximation and commutes with the Schrodinger operator. This is important in quantum mechanics.

Take

$$u_t=i\nabla^2u+i\sigma\vert u\vert^2u \qquad u(x,y,t) \in {\mathbb{C}}$$

Let ${\cal L}_1=i\nabla^2$      ${\cal L}_2u=i\sigma \vert u\vert^2u$. Here, $\nabla^2 = \partial_{xx} + \partial_{yy}$, in two space dimensions.

1. Advance linear part of NLS
   $$u_t=i\nabla^2u,$$
   i.e. using ${\cal L}_1$ and Fourier Transform, so that the $k^{th}$ spectral component advances as
   $$\hat u_k\left(t+\frac{\Delta t}{2}\right)=e^{-ik^2\frac{\Delta t}{2}}\hat u_k(t)$$
   using an FFT. Then take inverse FFT to obtain a quantity called ${\bf\bar u}\left(t+\Delta t/2\right)$.

2. To propagate under ${\cal L}_2$, solve
   $$\bar{\bar {\bf u}}_t=i\sigma\vert\bar {\bf u}\vert^2 \bar {\bf u}$$
   which has an exact solution since $\vert u\vert^2$ is conserved. The solution is
   $$\bar{\bar {\bf u}}(t+\Delta t)=e^{i\sigma}\Big\vert{\bf\bar{u}}\l... ...t}{2}\right)\Big\vert^2\Delta t \bar{\bf u} \left(t+\frac{\Delta t}{2}\right).$$

3. The final stage is another half-step propagation under ${\cal L}_1$.
   $$\hat u_k(t+\Delta t)=e^{-ik^2\frac{\Delta t}{2}} \hat{ \bar{ \bar {\bf u}}}_k(t+\Delta t).$$
   then inverse FFT of ${\bf\hat u}_k$ gives final ``solution'' after a single time step. Method requires 4 FFT's and 1 experimentation/time step and is $0(\Delta t^3)$ accurate. It is expensive if you do not use FFT's. See Tappert for more details.

$\Box$

ADI (Alternating Implicit Direction) Methods

A splitting method. Take

$$U_t=b_1U_{xx}+b_2U_{yy}$$

on a unit square.

Let $A_1 U=b_1 U_{xx}\quad A_2U=b_2 U_{yy}$

$$\mbox {then }U_t=A_1 U+A_2U$$

as before, supposing we used $\hat r_{1/1}$, with $k\equiv$ time step:

$$\left(I-\frac{k}{2}A_1, -\frac{k}{2}A_2\right){\bf u}^{n+1}=\left(I+\frac{k}{2}A_1+ \frac{k}{2}A_2\right){\bf u}^n+{\cal O}(b^3)$$

Since $(1\pm a_1)(1\pm a_2)=1\pm a_1\pm a_2+a_1 a_2$

we add $$\frac{k^2 A_1 A_2}{4}{\bf u}^{n+1}$$ to both sides

$$\left(I-\frac{k}{2}A_1-\frac{k}{2}A_2+\frac{k^2}{4}A_1 A_2\right){\bf u}^{n+1}=$$

$$\left(I+\frac{k}{2}A_1+\frac{k}{2}A_2+\frac{k^2}{4}A_1 A_2\right){\bf u}^n+\frac{k^2}{4}A_1 A_2({\bf u}^{n+1}{\bf -u}^n)+{\cal O}(k^3)$$

which can be factored

$\left(I-\displaystyle \frac{k}{2}A_1\right)\left(I-\displaystyle \frac{k}{2}A_... ...{({\bf u}^{n+1}-{\bf u}^n)}_ {{\cal O}(k)}\;\underbrace{{\cal O}(k^3)}_{0(k^3)}$

If $A_1$ and $A_2$ are discretized using the 2$^{nd}$ order stencil, we obtain tridiagonal matrices that are sparse and easily solved. Let $A_{1h}$ and $A_{2h}$ be $2^{nd}$ order approximations. Then

$$\left(I-\frac{k}{2}A_{1h}\right)\left(I-\frac{k}{2}A_{2h}\right){\bf u}^{n+1} = \left(I+\frac{k}{2}A_{1h}\right)(I+\frac{k}{2}\Delta_{2h}){\bf u}^n$$

$$+ {\cal O}(k^3)+{\cal O}(kh^2)$$

or

$$\left(I-\frac{k}{2}A_{1h}\right) \left(I-\frac{k}{2}A_{2h}\right)... ...12} =\left(I+\frac{k}{2}A_{1h} \right)\left(I+\frac{k}{2}A_{2h}\right){\bf u}^n$$ (157)

The Peaceman-Ratchford Algorithm to solve (158)

$$\left\{\begin{array}{c}\left(I-\frac{k}{2}A_{1h}\right) \tilde {\... ...lde {\bf u}^{n+\frac12}\\ \mbox {(unconditionally stable)}. \end{array}\right.$$

And we can see why it is called ADI $\cdots$

Another algorithm to solve (158):

Douglas-Rachford Algorithm: $(1,2)$ scheme, unconditionally stable.

$$\left\{\begin{array}{c}\left(I-\frac{k}{2}A_{1h}\right) {\bf\tild... ...rac12}-kA_{2h}{\bf u}^n\\ \mbox {(unconditionally stable)}. \end{array}\right.$$

Implementation Comments:

1. ADI schemes require intermediate values on the boundary. The approximation must be chosen so that no instabilities are introduced. Start with something simple.
2. The code can be made very fast if the row-column switching discussed previously is implemented.
3. Higher-order matrix representations or non-sparse representations may require an iterative technique for solution: for most cases an SOR or Conjugate Gradient (if symmetric) are quite efficient.

$\Box$