[PPT] - GSoC 2016: Exponential Integrators Chiara Segala Mentor: Prof. Marco PowerPoint Presentation

SLIDE 1

GSoC 2016: Exponential Integrators

Octave Conference 2017 March 21, 2017

GSoC 2016: Exponential Integrators

Chiara Segala Mentor: Prof. Marco Caliari

SLIDE 2

GSoC 2016: Exponential Integrators

Outline

Exponential Integrators Matrix functions Implementation of solvers Numerical experiments

SLIDE 3

GSoC 2016: Exponential Integrators Exponential Integrators

1 Exponential Integrators

2 Matrix functions

3 Implementation of solvers

4 Numerical experiments

SLIDE 4

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Integrators

Parabolic problem u′(t) + Au(t) = f (t), u(t0) = u0

SLIDE 5

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Integrators

Parabolic problem u′(t) + Au(t) = f (t), u(t0) = u0 Solution (variation of constants formula) u(tn+1) = e−hnAu(tn) + hn e−(hn−τ)Af (tn + τ)dτ

SLIDE 6

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Integrators

Parabolic problem u′(t) + Au(t) = f (t), u(t0) = u0 Solution (variation of constants formula) u(tn+1) = e−hnAu(tn) + hn e−(hn−τ)Af (tn + τ)dτ First numerical scheme (exponential quadrature rule) un+1 = e−hnAun + hn

s

i=1

bi(−hnA)f (tn + cihn) weights bi(z) can be expressed as linear combinations of the functions ϕk(z) = 1 e(1−θ)z θk−1 (k − 1)!dθ, k ≥ 1. ϕk+1(z) = ϕk(z) − ϕk(0) z , ϕk(0) = 1 k!, ϕ0(z) = ez

SLIDE 7

GSoC 2016: Exponential Integrators Exponential Integrators

Example

s = 1 un+1 = e−hnAun + hnϕ1(−hnA)f (tn + c1hn) = un + hnϕ1(−hnA)(f (tn + c1hn) − Aun)

SLIDE 8

GSoC 2016: Exponential Integrators Exponential Integrators

Example

s = 1 un+1 = e−hnAun + hnϕ1(−hnA)f (tn + c1hn) = un + hnϕ1(−hnA)(f (tn + c1hn) − Aun) Exponential Euler method for c1 = 0 Butcher tableau ϕ1 Exponential midpoint rule for c1 = 1

2

SLIDE 9

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Runge-Kutta methods

Semilinear problem u′(t) = F(t, u) = Au(t) + g(t, u(t)), u(t0) = u0

SLIDE 10

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Runge-Kutta methods

Semilinear problem u′(t) = F(t, u) = Au(t) + g(t, u(t)), u(t0) = u0 Numerical Exponential Runge-Kutta scheme un+1 = un + hn

s

i=1

bi(hnA)(Gni + Aun) Uni = un + hn

s

j=1

aij(hnA)(Gnj + Aun) Gnj = g(tn + cjhn, Unj)

SLIDE 11

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Runge-Kutta methods

Semilinear problem u′(t) = F(t, u) = Au(t) + g(t, u(t)), u(t0) = u0 Numerical Exponential Runge-Kutta scheme un+1 = un + hn

s

i=1

bi(hnA)(Gni + Aun) Uni = un + hn

s

j=1

aij(hnA)(Gnj + Aun) Gnj = g(tn + cjhn, Unj) Reformulation un+1 = un + hnϕ1(hnA)F(tn, un) + hn

s

i=2

bi(hnA)Dni Uni = un + hnciϕ1(cihnA)F(tn, un) + hn

i−1

j=2

aij(hnA)Dnj Dnj = g(tn + cjhn, Unj) − g(tn, un)

SLIDE 12

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Rosenbrock methods

u′(t) = F(t, u), u(t0) = u0

SLIDE 13

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Rosenbrock methods

u′(t) = F(t, u), u(t0) = u0 Numerical Exponential Rosenbrock scheme un+1 = un + hnϕ1(hnJn)F(tn, un) + h2

nϕ2(hnJn)vn + hn s

i=2

bi(hnJn)Dni Uni = un + hnciϕ1(cihnJn)F(tn, un) + h2

nc2 i ϕ2(cihnJn)vn + hn i−1

j=2

aij(hnJn)Dnj

SLIDE 14

GSoC 2016: Exponential Integrators Exponential Integrators

Exponential Rosenbrock methods

u′(t) = F(t, u), u(t0) = u0 Numerical Exponential Rosenbrock scheme un+1 = un + hnϕ1(hnJn)F(tn, un) + h2

nϕ2(hnJn)vn + hn s

i=2

bi(hnJn)Dni Uni = un + hnciϕ1(cihnJn)F(tn, un) + h2

nc2 i ϕ2(cihnJn)vn + hn i−1

j=2

aij(hnJn)Dnj where Jn = ∂F ∂u (tn, un) vn = ∂F ∂t (tn, un) gn(t, u) = F(t, u) − Jnu − vnt Dnj = gn(tn + cjhn, Unj) − gn(tn, un)

SLIDE 15

GSoC 2016: Exponential Integrators Matrix functions

1 Exponential Integrators

2 Matrix functions

3 Implementation of solvers

4 Numerical experiments

SLIDE 16

GSoC 2016: Exponential Integrators Matrix functions

ϕ-functions

ϕk(z) = 1 e(1−θ)z θk−1 (k − 1)!dθ, k ≥ 1. ϕk+1(z) = ϕk(z) − ϕk(0) z , ϕk(0) = 1 k!, ϕ0(z) = ez

SLIDE 17

GSoC 2016: Exponential Integrators Matrix functions

ϕ-functions

ϕk(z) = 1 e(1−θ)z θk−1 (k − 1)!dθ, k ≥ 1. ϕk+1(z) = ϕk(z) − ϕk(0) z , ϕk(0) = 1 k!, ϕ0(z) = ez ϕ-functions code: scaling and squaring + Padé approximation

SLIDE 18

GSoC 2016: Exponential Integrators Matrix functions

Accuracy

ϕk(±10) phikm(±10) phipade(±10,k) ϕ1(10) 2202.54657948068 2202.54657962449 ϕ2(10) 220.154657948068 220.154657962442 ϕ3(10) 21.9654657948068 21.9654657961828 ϕ4(10) 2.17987991281401 2.17987991294347 ϕ1(−10) 0.0999954600070237 0.0999954600070233 ϕ2(−10) 0.0900004539992977 0.0900004539992175 ϕ3(−10) 0.0409999546000702 0.0409999546000398 ϕ4(−10) 0.0125666712066596 0.0125666712066530

SLIDE 19

GSoC 2016: Exponential Integrators Matrix functions

Augmented matrix

Theorem (see Al-Mohy and Higham, 2011, Th. 2.1) Let A ∈ Cn×n, W = [w1, w2, . . . , wp] ∈ Cn×p, τ ∈ C, and

A =

A W J

∈ C(n+p)×(n+p),

J = Ip−1

∈ Cp×p,

Then for X = ϕl(τ A) with l ≥ 0 we have X(1 : n, n + j) =

j

k=1

τ kϕl+k(τA)wj−k+1, j = 1 : p.

SLIDE 20

GSoC 2016: Exponential Integrators Matrix functions

Augmented matrix for exponential integrators

W (:, p − k + 1) = uk, k = 1 : p, l = 0, τ = t − t0 X = ϕ0((t − t0) A) = e(t−t0)

A =

e(t−t0)A

X12 e(t−t0)J

X(1 : n, n + p) =

p

k=1

ϕk((t − t0)A)(t − t0)kuk

SLIDE 21

GSoC 2016: Exponential Integrators Matrix functions

Augmented matrix for exponential integrators

W (:, p − k + 1) = uk, k = 1 : p, l = 0, τ = t − t0 X = ϕ0((t − t0) A) = e(t−t0)

A =

e(t−t0)A

X12 e(t−t0)J

X(1 : n, n + p) =

p

k=1

ϕk((t − t0)A)(t − t0)kuk ˆ u(t) = e(t−t0)Au0 +

p

k=1

ϕk((t − t0)A)(t − t0)kuk = e(t−t0)Au0 + X(1 : n, n + p) = [In 0] e(t−t0)

A u0

ep

SLIDE 22

GSoC 2016: Exponential Integrators Matrix functions

Augmented matrix for exponential integrators

W (:, p − k + 1) = uk, k = 1 : p, l = 0, τ = t − t0 X = ϕ0((t − t0) A) = e(t−t0)

A =

e(t−t0)A

X12 e(t−t0)J

X(1 : n, n + p) =

p

k=1

ϕk((t − t0)A)(t − t0)kuk ˆ u(t) = e(t−t0)Au0 +

p

k=1

ϕk((t − t0)A)(t − t0)kuk = e(t−t0)Au0 + X(1 : n, n + p) = [In 0] e(t−t0)

A u0

ep

ˆ

u(t) = [In 0] exp

(t − t0)

A ηW J u0 η−1ep

,

η = 2−⌈log2(W 1)⌉

SLIDE 23

GSoC 2016: Exponential Integrators Matrix functions

Double augmented matrix

Theorem (Double augmented matrix) Let A ∈ Cn×n, W = [w1, w2, . . . , wp] ∈ Cn×p, V = [v1, v2, . . . , vq] ∈ Cn×q, τ ∈ C, and

A =

A W V J K

∈ C(n+p+q)×(n+p+q),

J ∈ Cp×p, K ∈ Cq×q. Then for X = ϕl(τ

A ) with l ≥ 0 we have

X(1 : n, n + j) =

j

k=1

τ kϕl+k(τA)wj−k+1, j = 1 : p. X(1 : n, n + i) =

i−p

k=1

τ kϕl+k(τA)vi−k+1−p, i = (p + 1) : (p + q).

SLIDE 24

GSoC 2016: Exponential Integrators Matrix functions

expmv

eAB ≈ (Tm(s−1A))sB m = ⇒ select_taylor_degree

SLIDE 25

GSoC 2016: Exponential Integrators Matrix functions

expmv

eAB ≈ (Tm(s−1A))sB m = ⇒ select_taylor_degree function [f] = expmv (A, b) . . . M = select_taylor_degree (A, b); . . .

SLIDE 26

GSoC 2016: Exponential Integrators Matrix functions

expmv

eAB ≈ (Tm(s−1A))sB m = ⇒ select_taylor_degree function [f] = expmv (A, b) . . . M = select_taylor_degree (A, b); . . . M = select_taylor_degree (A, b); f = expmv (A, b, M);

SLIDE 27

GSoC 2016: Exponential Integrators Matrix functions

expmv

eAB ≈ (Tm(s−1A))sB m = ⇒ select_taylor_degree function [f] = expmv (A, b) . . . M = select_taylor_degree (A, b); . . . M = select_taylor_degree (A, b); f = expmv (A, b, M); select_taylor_degree = ⇒ Ak1

SLIDE 28

GSoC 2016: Exponential Integrators Matrix functions

Estimate

Theorem (Estimate) Let A = A W J

∈ C(n+p)×(n+p) be the matrix of Higham Theorem. Suppose

that W 1 = γ ≤ 1. Then we have the following estimate

Ak1 ≤ max
Ak1, γ · k−1

j=k−min{k,p} Aj1 + max{min{p − k, 1}, 0}

.

SLIDE 29

GSoC 2016: Exponential Integrators Matrix functions

Estimate

Theorem (Estimate) Let A = A W J

∈ C(n+p)×(n+p) be the matrix of Higham Theorem. Suppose

that W 1 = γ ≤ 1. Then we have the following estimate

Ak1 ≤ max
Ak1, γ · k−1

j=k−min{k,p} Aj1 + max{min{p − k, 1}, 0}

.

SLIDE 30

GSoC 2016: Exponential Integrators Implementation of solvers

1 Exponential Integrators

2 Matrix functions

3 Implementation of solvers

4 Numerical experiments

SLIDE 31

GSoC 2016: Exponential Integrators Implementation of solvers

dexprk23

c2 c2ϕ1,2

2 3 2 3ϕ1,3 − 8 9ϕ2,3 8 9ϕ2,3

(1 −

1 2c2 )ϕ1 1 2c2 ϕ1

ϕ1 − 3

2ϕ2 3 2ϕ2

ϕj,k = ϕj(ckhA) and ϕj = ϕj(hA)

SLIDE 32

GSoC 2016: Exponential Integrators Implementation of solvers

dexprk23

c2 c2ϕ1,2

2 3 2 3ϕ1,3 − 8 9ϕ2,3 8 9ϕ2,3

(1 −

1 2c2 )ϕ1 1 2c2 ϕ1

ϕ1 − 3

2ϕ2 3 2ϕ2

ϕj,k = ϕj(ckhA) and ϕj = ϕj(hA) Un1 = un Un2 = un + 1 2hnϕ1(1 2hnA)F(tn, un) Un3 = un + 2 3hnϕ1(2 3hnA)F(tn, un) + 2 3hnϕ2(2 3hnA)2 3(2)Dn2 un+1 = un + hnϕ1(hnA)F(tn, un) + hnϕ2(hnA)3 2Dn3 ˆ un+1 = un + hnϕ1(hnA)[F(tn, un) + 1 2(2)Dn2]

SLIDE 33

GSoC 2016: Exponential Integrators Implementation of solvers

exp_runge_kutta_23

function [t_next, x_next, x_est, k] = exp_runge_kutta_23 (f, t, x, dt,

ptions,

k_vals = [], t_next = t + dt) x = x.'; M = options.GetInfo; myexpmv = @(t,A,v) expmv (t, A, v, M); % myexpmv = @(t,A,v) expmv (t, A, v); myexpmv = @(t,A,v) expm (t*A)*v; A = options.LinOp; g = options.GFun; d = size (A, 1); c2 = 1/2; k = zeros (length (x), 3); Fn = feval (f, t, x(:)); gn = feval (g, t, x(:)); [Atilde, eta] = augmatrix ( A, Fn); X = feval(myexpmv, c2*dt, Atilde, [zeros(size (Atilde)−1, 1); 1/eta])(1:d, :); U2 = x(:) + X; g2 = feval (g, t + c2*dt, U2); D2 = g2 − gn; [Atilde, eta] = augmatrix ( A, [Fn, (2/(3*c2)*D2)./(2/3*dt)]); X = feval(myexpmv, 2/3*dt, Atilde, [zeros(size (Atilde)−1, 1); 1/eta])(1:d, :); U3 = x(:) + X; g3 = feval (g, t + 2/3*dt, U3); D3 = g3 − gn; [Atilde, eta] = augmatrix ( A, [Fn, (3/2*D3)./dt], [Fn + 1/(2*c2)*D2]); ep = [[zeros(1, 1); 1/eta; 0], [zeros(2, 1); 1/eta]]; X = feval(myexpmv, dt, Atilde, [zeros(d, 1), zeros(d, 1); ep])(1:d, :); x_next = x(:) + X(:,1); x_est = x(:) + X(:,2); k(:, 1) = Fn; k(:, 2) = A*U2 + g2; k(:, 3) = A*U3 + g3;

SLIDE 34

GSoC 2016: Exponential Integrators Implementation of solvers

dexprb34

1 2 1 2ϕ1( 1 2·)

1 ϕ1 ϕ1 − 14ϕ3 + 36ϕ4 16ϕ3 − 48ϕ4 −2ϕ3 + 12ϕ4 ϕ1 − 14ϕ3 16ϕ3 −2ϕ3 ϕj = ϕj(hJn)

SLIDE 35

GSoC 2016: Exponential Integrators Implementation of solvers

dexprb34

1 2 1 2ϕ1( 1 2·)

1 ϕ1 ϕ1 − 14ϕ3 + 36ϕ4 16ϕ3 − 48ϕ4 −2ϕ3 + 12ϕ4 ϕ1 − 14ϕ3 16ϕ3 −2ϕ3 ϕj = ϕj(hJn) Un1 = un Un2 = un + 1 2hnϕ1(1 2hnJn)F(tn, un) + 1 2hnϕ2(1 2hnJn)1 2hnvn Un3 = un + hnϕ1(hnJn)[F(tn, un) + Dn2] + hnϕ2(hnJn)hnvn un+1 = un + hnϕ1(hnJn)F(tn, un) + hnϕ2(hnJn)hnvn + hnϕ3(hnJn)[16Dn2 − 2Dn3] + hnϕ4(hnJn)[−48Dn2 + 12Dn3] ˆ un+1 = un + hnϕ1(hnJn)F(tn, un) + hnϕ2(hnJn)hnvn + hnϕ3(hnJn)[16Dn2 − 2Dn3]

SLIDE 36

GSoC 2016: Exponential Integrators Implementation of solvers

exp_rosenbrock_34

function [t_next, x_next, x_est, k] = exp_rosenbrock_34 (f, t, x, dt,

ptions,

k_vals = [], t_next = t + dt) x = x.'; k = zeros (length (x), 4); if (isempty (options.Jacobian)) fu = @(xx) f(t,xx); Jn = numjac ( x(:), fu); else J = options.Jacobian; Jn = feval (J, t, x(:)); end M = select_taylor_degree_aug (Jn, [], 2, [], [], [], false, false, true); myexpmv = @(t,A,v) expmv (t, A, v, M); d = size (Jn, 1); Fn = feval (f, t, x(:)); if (isempty (options.DFdt)) ft = @(tt) f(tt, x(:)); vn = numjac ( t, ft); else V = options.DFdt; vn = feval (V, t, x(:)); end g = @(t,x) f(t,x) − Jn*x − vn*t; gn = feval (g, t, x(:)); [Atilde, eta] = augmatrix ( Jn, [Fn, 1/2*vn]); X = feval(myexpmv, 1/2*dt, Atilde, [zeros(size (Atilde)−1, 1); 1/eta])(1:d, :); U2 = x(:) + X; g2 = feval (g, t + 1/2*dt, U2); D2 = g2 − gn; [Atilde, eta] = augmatrix ( Jn, [Fn+D2, vn]); X = feval(myexpmv, dt, Atilde, [zeros(size (Atilde)−1, 1); 1/eta])(1:d, :); U3 = x(:) + X; g3 = feval (g, t + dt, U3); D3 = g3 − gn; M = select_taylor_degree_aug (Jn, [], 4, [], [], [], false, false, true); myexpmv = @(t,A,v) expmv (t, A, v, M); [Atilde, eta] = augmatrix ( Jn, [Fn, vn, (16*D2−2*D3)./(dt^2), (−48*D2+12*D3)./(dt^3)], [Fn, vn, (16*D2−2*D3)./(dt^2)]); ep = [[zeros(3, 1); 1/eta; zeros(3, 1)], [zeros(6, 1); 1/eta]]; X = feval(myexpmv, dt, Atilde, [zeros(d, 1), zeros(d, 1); ep])(1:d, :); x_next = x(:) + X(:, 1); x_est = x(:) + X(:, 2); k(:, 1) = Fn; k(:, 2) = feval (f, t + 1/2*dt, U2); k(:, 3) = feval (f, t + dt, U3);

SLIDE 37

GSoC 2016: Exponential Integrators Implementation of solvers

Derivative approximation

u′ = F(u) u(t0) = u0 J = ∂F

∂u (u0)

δ = √eps = ⇒ J(:, j) ≈

F(u0+δej )−F(u0) δ

= ⇒ J(:, j) ≈ Im(F(u0 + iδej)/δ)

SLIDE 38

GSoC 2016: Exponential Integrators Numerical experiments

1 Exponential Integrators

2 Matrix functions

3 Implementation of solvers

4 Numerical experiments

SLIDE 39

GSoC 2016: Exponential Integrators Numerical experiments

Semilinear parabolic problem

Semilinear parabolic problem ∂U ∂t (x, t) − ∂2U ∂x2 (x, t) = 1 1 + U(x, t)2 + Φ(x, t) (1)

SLIDE 40

GSoC 2016: Exponential Integrators Numerical experiments

Semilinear parabolic problem

Semilinear parabolic problem ∂U ∂t (x, t) − ∂2U ∂x2 (x, t) = 1 1 + U(x, t)2 + Φ(x, t) (1) x ∈ [0, 1] t ∈ [0, 1] homogeneous Dirichlet boundary conditions Φ chosen such that the exact solution is U(x, t) = x(1 − x)et space discretization: standard finite differences with 200 grid points we expect to see order three for odexprk23 and order four for odeprb34

SLIDE 41

GSoC 2016: Exponential Integrators Numerical experiments

Example

% RUNGE KUTTA a = 0; b = 1; m = 200; x = linspace (a, b, m+2)'; x = x(2:m+1); h = (b − a)/(m + 1); A = toeplitz (sparse([1, 1], [1, 2], [−2, 1]/h^2, 1, m)); g = @(t,U) [ 1./(1 + U.^2) + exp(t)*(2 + x − x.^2) − 1./(1 + x.^2.*(1 − x).^2 * exp(2*t))]; F = @(t,U) A*U + g(t,U); U0 = (x.*(1 − x)); trange = [0; 1];

pt = odeset ('AbsTol',

10^−2, 'GFun', g, 'LinOp', A, 'RelTol', 10^−2); [t,U] = odexprk23 (F, trange, U0, opt); % ROSENBROCK a = 0; b = 1; m = 200; x = linspace (a, b, m+2)'; x = x(2:m+1); h = (b − a)/(m + 1); A = toeplitz (sparse([1, 1], [1, 2], [−2, 1]/h^2, 1, m)); g = @(t,U) [ 1./(1 + U.^2) + exp(t)*(2 + x − x.^2) − 1./(1 + x.^2.*(1 − x).^2 * exp(2*t))]; F = @(t,U) A*U + g(t,U); J = @(t,U) A − diag ([ 2*U./(1 + U.^2).^2 ]); V = @(t,U) [exp(t)*(2 + x − x.^2) + ( 2*x.^2.*(1 − x).^2 * exp(2*t) )./( (1 + x.^2.*(1 − x).^2 * exp(2*t)).^2 )]; U0 = (x.*(1 − x)); trange = [0; 1];

pt = odeset ('AbsTol',

10^−2, 'DFdt', V, 'Jacobian', J, 'RelTol', 10^−2); [t,U] = odexprb34 (F, trange, U0, opt);

SLIDE 42

GSoC 2016: Exponential Integrators Numerical experiments

Convergence

Test with semilinear parabolic problem (1) Runge-Kutta method Rosenbrock method

SLIDE 43

GSoC 2016: Exponential Integrators Numerical experiments

Schrödinger equation

i ∂ψ ∂t = H(x, t)ψ H(x, t) = −1 2 ∂2 ∂x2 + k x2 2 + µ(sin t)2x, Runge-Kutta method Rosenbrock method

SLIDE 44

GSoC 2016: Exponential Integrators Numerical experiments

expmv.m and Oct-Files

%.cc file DEFUN_DLD (taylor, args, , "Computing the Action of the Matrix Exponential") %.m file f = taylor (s, t, A, b, ...);

SLIDE 45

GSoC 2016: Exponential Integrators Numerical experiments

expmv.m and Oct-Files

%.cc file DEFUN_DLD (taylor, args, , "Computing the Action of the Matrix Exponential") %.m file f = taylor (s, t, A, b, ...);

SLIDE 46

GSoC 2016: Exponential Integrators Numerical experiments