Solving partial differential equations with exponential integrators

Jitse Niesen in collaboration with Will Wright (La Trobe University, Melbourne)

Applied maths seminar, Leeds, 21 January 2008

Outline I

The problem I I

I

Exponential integrators I I I I

I

Exponential Euler method Construction of higher-order methods A predictor-corrector method Example using Fourier transform

Computing the matrix exponential I I I

I

Solving PDEs via semi-discretization Stiffness and its effects

Pad´e with scaling and squaring Krylov methods Example using finite differences

Conclusions

Solving ordinary differential equations Consider the ordinary differential equation u 0 = f (u),

u(t0 ) = u0 .

The derivative can be approximated by u 0 (t) ≈

u(t + h) − u(t) h

(h small).

This leads to the Euler method: un+1 = un + hf (un ). The solution of the ODE at time t = nh is approximated by un .

Solving partial differential equations The Allen–Cahn equation:

∂2u ∂u = + u − u3. ∂t ∂x 2

t2 t1 h t0

x0

∆x

x1

x2

· · · xM−1 xM xM+1

 ∂2 1  u(x, t) ≈ u(x − ∆x, t) − 2u(x, t) + u(x + ∆x, t) ∂x 2 ∆x 2

Solving partial differential equations II This approximation converts the PDE into a system of ODEs:      0   u1 − u13 u1 u1 −2 1    u2    u0   1 −2 1 u2 − u23      2  1        ..   .. .. .. .. .. +       . =  . . . . . 2     ∆x    0 3      uM−1   uM−1 − uM−1  uM−1 1 −2 1 3 0 uM 1 −2 uM − uM uM This ODE can be solved by the Euler method. However, the ODE is stiff: Df has very negative eigenvalues. The most negative one is ≈ −4/∆x 2 . 2

∂ Compare σ( ∂x 2 ) = (−∞, 0].

Stiffness and its effects Euler method: un+1 = un + hf (un ). If you solve a stiff ODE with the Euler method then you have to pick a sufficiently small time step h, otherwise the numerical result becomes unstable. In the example, h < 12 ∆x 2 ; this is the CFL condition. Backward Euler method: un+1 = un + hf (un+1 ). The backward Euler method does not suffer from stiffness, but requires solving a system of nonlinear equations. This system can get very large, especially in multiple dimensions. Better methods than (backward) Euler exist, but the basic problem remains the same: How to solve a large system of stiff ODEs?

Outline I

The problem I I

I

Exponential integrators I I I I

I

Exponential Euler method Construction of higher-order methods A predictor-corrector method Example using Fourier transform

Computing the matrix exponential I I I

I

Solving PDEs via semi-discretization Stiffness and its effects

Pad´e with scaling and squaring Krylov methods Example using finite differences

Conclusions

Exponential Euler method u 0 = Lu + N(u)

To solve: Variation of constants

tL

Z

t

u(t) = e u(0) + ≈ etL u(0) +

Z0 t

e(t−τ )L N(u(τ )) dτ e(t−τ )L N(u(0)) dτ

0

etL − 1 N(u(0)). = etL u(0) + tL This leads to the exponential Euler method un+1 = ehL un +

ehL − 1 N(un ). hL

This method is not affected by stiffness in L.

(Certaine 1960)

How to achieve higher order The nonautonomous differential equation u 0 (t) = Lu(t) + v0 +

1 1 tv1 + t 2 v2 + · · · 1! 2!

has exact solution u(t) = etL u(0) + hϕ1 (tL)v0 + h2 ϕ2 (tL)v1 + h3 ϕ3 (tL)v2 + · · · The ϕ functions are defined by ex − 1 x ex − 1 − x ϕ2 (x) = x2 x e − 1 − x − 21 x 2 ϕ3 (x) = x3

ϕ1 (x) =

A predictor-corrector method We are using predictor-corrector methods like I

Predictor is exponential Adams–Bashforth:
∗ un+1 = ehL un + hϕ1 (hL)N(un ) + h2 ϕ2 (hL) N(un ) − N(un−1 )

I

Corrector based on exponential Adams–Moulton: ∗ ∗ + h2 ϕ2 (hL) N(un+1 ) − 2N(un ) + N(un−1 ) un+1 = un+1

I

A posteriori error estimate, used to change step size h

This method has order 2, i.e., error = O(h2 ). We have constructed methods up to order 6. Other exponential methods also exist (without corrector step, based on Runge–Kutta, using exact Jacobian).




Evaluating ϕ functions for scalars Using the definition, e.g., ϕ2 (x) = (ex − 1 − x)/x 2 , leads to loss of precision for small x due to subtraction of two nearly equal numbers. Therefore, use Taylor or Pad´e approximations: ϕ2 (x) ≈

c0 + c1 x + c2 x 2 + c3 x 3 d0 + d1 x + d2 x 2 + d3 x 3

when |x| small.

1 . Also combine ϕ functions using ϕk (z) = zϕk+1 (z) + k! hL For example, e v0 + ϕ1 (hL)v1 + ϕ2 (hL)v2 = ϕ2 (hL)w0 + w1 .

But what to do for matrices? Sometimes, L can be diagonalized. If L approximates Laplacian ∇2 with periodic boundary conditions, L = F −1 DF =⇒ ϕ(L) = F −1 ϕ(D)F.

Gray–Scott equation ut = du ∇2 u − uv 2 + f (1 − u), 2

2

vt = dv ∇ v + uv − (f + k)v with periodic boundary conditions.

x, y ∈ [0, 1]

Outline I

The problem I I

I

Exponential integrators I I I I

I

Exponential Euler method Construction of higher-order methods A predictor-corrector method Example using Fourier transform

Computing the matrix exponential I I I

I

Solving PDEs via semi-discretization Stiffness and its effects

Pad´e with scaling and squaring Krylov methods Example using finite differences

Conclusions

Computing the matrix exponential ϕ functions for matrices are like the matrix exponential eA . “Least dubious” way uses Pad´e approximation (valid when kAk small) with scaling and squaring: (Moler & Van Loan 1978) 1. Scaling: Find least k such that kBk < 5.4 with B = 2−k A. c0 + c1 B + c2 B 2 + · · · + c13 B 13 2. Pad´ e: eB ≈ . d0 + d1 B + d2 B 2 + · · · + d13 B 13 3. Squaring: Compute eA by squaring eB k times. [e2B = (eB )2 ] Constants 5.4 and 13 are optimal.

(Higham 2005)

Use same method for ϕ functions with modified squaring step:   k X 1 1 ϕj (B) . ϕk (2B) = k ϕ0 (B)ϕk (B) + (k − j)! 2 j=1

Krylov method Pad´e approximation with squaring and scaling is too costly for large matrices (N ' 103 ). Idea of Krylov-subspace methods is to project the action of A on the Krylov subspace with dimension m  N: span{v , Av , A2 , . . . , Am−1 v }. The Lanczos (if A symmetric) or Arnoldi algorithm produces an orthogonal N-by-m matrix Vm and an m-by-m matrix Hm such that A ≈ Vm Hm VmT =⇒ ϕ(A) ≈ Vm ϕ(Hm )VmT . Hm is small so we can compute ϕ(Hm ) via Pad´e. (Friesner, Tuckerman & Dornblaser 1989)

Adaptivity in the Krylov process Recall: The nonautonomous differential equation u 0 (t) = Lu(t) + v0 +

1 1 tv1 + t 2 v2 + · · · 1! 2!

has exact solution u(t) = etL u(0) + hϕ1 (tL)v0 + h2 ϕ2 (tL)v1 + h3 ϕ3 (tL)v2 + · · · This can be used to split the computation of ϕ(A) in substeps. We have two parameters to play with: I

The dimension m of the Krylov subspace.

I

The size of the substeps.

Adapt them to achieve the required precision with the least effort, using an error estimate due to Saad (1992) and Sidje (1998).

Overview of the algorithm To solve: u 0 = Lu + N(u) over time interval [0, T ]. 1. Split [0, T ] in smaller intervals. 2. Three stages per step: predictor, corrector, error estimate. 3. Every stage requires computation of a linear combination of ϕ functions, which we reduces to one ϕ function. 4. Split the computation of the ϕ function in substeps. 5. Reduce the size of the matrix by going to the Krylov subspace. 6. Use Pad´e with scaling and squaring to compute the ϕ function of the smaller matrix. We use adaptivity at 1, 4 and 5.

Example: Allen–Cahn in two dimensions The equation we considered at the beginning: ut = α∇2 u + u − u 3 ,

x, y ∈ [0, 1],

t ∈ [0, 1]

Neumann boundary conditions, α = 0.1, 100 × 100 grid. What is the effect of adapting the Krylov subspace dimension m?

steps substeps L·v eA time error

m = {15, 15, 15} m = {32.4, 9.7, 3.1} m = {32, 10, 3} (adaptive) 34 35 34 601 114 1371 9015 1872 10871 679 490 2101 39.41 12.98 58.94 1.5 · 10−4 1.5 · 10−4 1.6 · 10−4

Example: Allen–Cahn in two dimensions II 2D Allen−Cahn α = 0.1

2

10

odeexp5 oderos4 ode15s radau5 exp4

0

10

−2

global error

10

−4

10

−6

10

−8

10

−10

10

−12

10

0

10

1

10 cpu time

2

10

Implicit: ode15s (Matlab), radau5 (Hairer & Wanner) Exponential: odeexp5 (us), oderos4 (Caliari & Ostermann), exp4 (Hochbruck, Lubich & Selhofer)

Discussion I

All methods except ours are given the exact Jacobian.

I

Exponential methods require rhs = Lu + N(u) with stiffness concentrated in L. Performance depends on stiffness:

I

I I

I

I

Explicit methods best on non-stiff problems. Exponential integrators promising on mildly stiff problems, but not enough data to be sure. Implicit methods better on very stiff problems ? (Krylov approximation converges slowly)

Alternative methods for evaluating the ϕ function: I I I I

Leja point interpolation Using contour integrals RD-rational approximations Rational Krylov

(Caliari & Ostermann) (Schmelzer & Trefethen) (Moret & Novati 2004) (Grimm & Hochbruck)

Implicit-explicit (IMEX) method Our comparison did not include IMEX methods. Recall backward Euler method: un+1 = un + hf (un+1 ). Use the splitting f (u) = Lu + N(u) to get
un+1 = un + h Lun+1 + N(un+1 ) . This is a nonlinear equation. However, if stiffness is in linear part, then why not use the linear equation
un+1 = un + h Lun+1 + N(un ) . Can be solved with Krylov techniques. What is faster: ehL v or (I − hL)−1 v ? Theory suggests ehL v . However, in many situations preconditioning helps a lot for (I − hL)−1 v .

Conclusion

Many exponential integrators have been constructed. Some can compete with state-of-the-art implicit methods. However, some important question are still open: I

How do exponential integrators compare to implicit and IMEX methods based on subspace techniques?

I

How do exponential integrators fare on real-world problems?

Conclusion

Many exponential integrators have been constructed. Some can compete with state-of-the-art implicit methods. However, some important question are still open: I

How do exponential integrators compare to implicit and IMEX methods based on subspace techniques?

I

How do exponential integrators fare on real-world problems?

HELP!