Solving Ill-Posed Cauchy Problems using Rational Krylov Methods - - PowerPoint PPT Presentation

Solving Ill-Posed Cauchy Problems using Rational Krylov Methods

Lars Eldén and Valeria Simoncini GAMM 2008

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 1 / 60

Outline

1

Cauchy Problem Regularization Error Estimate

2

Rational Krylov method Krylov/Regularization

3

Numerical Examples Example 1 Animation Example 2

4

Conclusions Future work

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 2 / 60

Motivating example: Ilmenite iron melting furnace

Thermocouple

Electrode Level K to level D thermocouples highest with level K Under electrodes Between electrodes Center

The furnace material properties are temperature dependent. Problem: find the inner shape of the furnace. Nonlinear, and (rather) complex geometry

PhD thesis: I M Skaar, Monitoring the Lining of a Melting Furnace, NTNU, Trondheim, 2001

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 3 / 60

Inverse Heat Conduction Problem

Steady state heat conduction problem: The upper boundary is unavailable for measurements

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 4 / 60

Ill-Posed Cauchy Problem

Ω: connected in R2 with smooth boundary ∂Ω L: linear, self-adjoint, positive definite elliptic in Ω. uzz − Lu = 0, (x, y) ∈ Ω, z ∈ [0, z1], u(x, y, z) = 0, (x, y) ∈ ∂Ω, z ∈ [0, z1], u(x, y, 0) = g(x, y), (x, y) ∈ Ω, uz(x, y, 0) = 0, (x, y) ∈ Ω. Sought: f(x, y) = u(x, y, z1), (x, y) ∈ Ω.

Formal solution

u(x, y, z) = cosh(z √ L)g BUT: L is a pos. def. unbounded operator ⇒ ILL-POSED!

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 5 / 60

Standard Iterative Procedure

Guess f (1)

1

for k = 1, 2, . . . until convergence

1

Solve uzz − Lu = 0, (x, y) ∈ Ω, z ∈ [0, z1], u(x, y, z) = 0, (x, y) ∈ ∂Ω, z ∈ [0, z1], u(x, y, z1) = f (k), (x, y) ∈ Ω, uz(x, y, 0) = 0, (x, y) ∈ Ω. giving u(k)

2

Evaluate g(·, ·) − u(k)(·, ·, 0) and adjust f (k) − → f (k+1)

In every iteration: Solve a 3D well-posed problem Often slow convergence.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 6 / 60

Other possible methods?

Tikhonov regularization? Impossible, because we do not know the integral operator for equations with variable coefficients and/or complicated geometry. Replace unbounded L by a bounded approximation? Possible in connection with finite difference approximation, but more difficult with finite elements. BUT: Krylov method!

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 7 / 60

Regularization

Formal solution

u(x, y, z) = cosh(z √ L)g BUT: L is a pos. def. unbounded operator ⇒ ILL-POSED! High frequency perturbations in g are blown up

Regularization

Replace unbounded operator by a bounded one!

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 8 / 60

Cut Off High Frequencies

Eigenvalues of L: λ2

j , j = 1, 2, . . . and λj → +∞ as j → ∞

General approach: Compute the k eigenvalues of smallest modulus: LXk = XkDk, where Xk holds orthonormal eigenvectors

Approximate by projection

cosh(z √ L)g ≈ cosh(z √ L)XkX ⊤

k g = Xk cosh(z

• Dk)X ⊤

k g

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 9 / 60

Eigenvalues of L

L is large and sparse (of the order 104 − 105, say) Compute the smallest eigenvalues ⇓ Operate with L−1 ⇓ Solve many standard 2D elliptic problems −Lw = v

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 10 / 60

Error Estimate

L2(Ω) setting, u is an “exact solution” v is an approximate solution with perturbed data gm

Theorem

Assume that u(·, ·, 1) ≤ M and that the data perturbation satisfies g − gm ≤ ǫ. Then if v is computed by projection using the eigenvalues satisfying λj ≤ λc, where λc = (1/z1) log(M/ǫ) then u(·, ·, z) − v(·, ·, z) ≤ 3ǫ1−z/z1Mz/z1, 0 ≤ z ≤ z1. Optimal error bound

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 11 / 60

Eigenvalue method

Is it necessary to compute the eigenvalues and eigenvectors accurately? Do we need all the information that we get in the eigenvalues? Can we take advantage of the fact that we want to compute an approximation of cosh(z √ L)g for this particular vector?

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 12 / 60

Eigenvalue method

Is it necessary to compute the eigenvalues and eigenvectors accurately? NO! Do we need all the information that we get in the eigenvalues? NO! Can we take advantage of the fact that we want to compute an approximation of cosh(z √ L)g for this particular vector? YES!

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 13 / 60

Lanczos tridiagonalization

Choose q1 and iterate L−1qk = qk−1βk−1 + qkαk + qk+1βk, k = 1, 2, . . . , with αk = q⊤

k L−1qk and βk = q⊤ k+1L−1qk;

One matrix-vector multiply L−1qk per step One standard 2D elliptic solve (black box) per step −Lw = qk

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 14 / 60

Lanczos properties

Initial convergence influenced by starting vector v1. Choose q1 = 1/β gm, β = gm Faster for largest eigenvalues of L−1 ⇒ Fast convergence for some of the smallest eigenvalues of L Optional: To get faster convergence for eigenvalues in [0, λc]

• perate with (L − τI)−1, where τ = λc/2

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 15 / 60

Lanczos reduction

L−1Qk = QkTk + βk+1qk+1e⊤

k+1 ≈ QkTk

Approximation

cosh(z √ L)g ≈ Qk cosh(zT −1/2

k

)Q⊤

k g

Problem: We cannot prevent Tk from approximating large eigenvalues! Solution: Regularize Tk: cut off large eigenvalues1

1Krylov+regularization: O’Leary & Simmons (1981), Björck, Grimme & Van Dooren

(1994)

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 16 / 60

Projected and Truncated Approximation

Let ((θ(k)

j

)2, y(k)

j

), j = 1, . . . , k be the eigenpairs of T −1

k

Define F(z, λ) = cosh(zλ1/2) and Sk = T −1

k

Truncated approximation: vk(z) = QkF(z, Sc

k )Q⊤ k gm

:= Qk

• θ(k)

j

≤λc

y(k)

j

cosh(zθ(k)

j

)(y(k)

j

)⊤e1gm.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 17 / 60

Error Estimate

Recall F(z, λ) = cosh(zλ1/2) and Sk = T −1

k

Theorem

Let u be the “exact solution” and vk(z) = QkF(z, Sc

k )Q⊤ k gm

Under the same hypotheses as earlier, u(z) − vk(z) ≤ 3ǫ1−z/z1Mz/z1 + 2[F(z, Lc) − QkF(z, Sc

k )Q⊤ k ]g.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 18 / 60

Krylov/Regularization, First version

Starting vector q1 = 1/β gm for k = 2, 3, . . . until “stable”

[Qk, Tk] = krylovstep(L−1, Qk−1, Tk−1) Compute vk(z) = QkF(z, Sc

k)Q⊤ k gm

end Check residual Kvk − gm < ǫ Kvk is the solution of the 3D problem with u = vk at the upper boundary and uz = 0 at the lower. Expensive!

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 19 / 60

Residual: Kvk − gm < ǫ

Solve 3D problem (denote solution uk) uzz − Lu = 0, (x, y) ∈ Ω, z ∈ [0, z1], u(x, y, z) = 0, (x, y) ∈ ∂Ω, z ∈ [0, z1], u(x, y, 1) = vk(x, y), (x, y) ∈ Ω, uz(x, y, 0) = 0, (x, y) ∈ Ω. Well-posed but expensive! Kvk = uk(x, y, 0) We only want to compute this when we are sure that Kvk − gm < ǫ

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 20 / 60

Krylov/Regularization: “stable”

Starting vector q1 = 1/β gm for k = 1, 2, . . . until “stable”

[Qk, Tk] = krylovstep(L−1, Qk−1, Tk−1) Compute vk(z) = QkF(z, Sc

k)Q⊤ k gm

end Check residual Kvk − gm < ǫ How can we quantify “stable”?

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 21 / 60

Error Estimate for Krylov Procedure

Recall F(z, λ) = cosh(zλ1/2) and Sk = T −1

k

u(z) − vk(z) ≤ 3ǫ1−z/z1Mz/z1 + 2[F(z, Lc) − QkF(z, Sc

k )Q⊤ k ]g.

Second term: Krylov approximation error for low frequency operator applied to the exact right hand side!

Heuristic

When the approximate solution does not change much between successive Krylov steps, then this error is small.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 22 / 60

Krylov/Regularization

Recall F(z, λ) = cosh(zλ1/2) and Sk = T −1

k

Starting vector q1 = 1/β gm for k = 1, 2, . . . maxit

[Qk, Tk] = krylovstep(L−1, Qk−1, Tk−1) Compute wk(z) = F(z, Sc

k)Q⊤ k gm

if wk − wk−1 < tol then

if Kvk − gm < ǫ then stop iterating

endif

end Compute vk = Qkwk

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 23 / 60

Test example 1: Laplace equation

Ω: unit square uzz + ∆u = 0, (x, y, z) ∈ Ω × [0, 0.1], u(x, y, z) = 0, (x, y, z) ∈ ∂Ω × [0, 0.1], u(x, y, 0) = g(x, y), (x, y) ∈ Ω, uz(x, y, 0) = 0, (x, y) ∈ Ω. Determine the values at the upper boundary, f(x, y) = u(x, y, 0.1), (x, y) ∈ Ω. Data perturbation: g − gm/g ≈ 0.0085 98 eigenvalues are smaller than the tolerance eigs performs approximately 300 2D elliptic solves.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 24 / 60

Solution and Exact Data

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 25 / 60

Convergence History (0.75λc)

5 10 15 20 25 30 35 40 10

−3

10

−2

10

−1

10 step Difference True error Residual Spectral

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 26 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 3

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 27 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 4

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 28 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.1 0.1 0.2 0.3 0.4 0.5 0.6

k= 5

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 29 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 6

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 30 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 7

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 31 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 8

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 32 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.1 0.1 0.2 0.3 0.4 0.5 0.6

k= 9

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 33 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 10

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 34 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 11

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 35 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 12

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 36 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 13

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 37 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 14

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 38 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 15

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 39 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 16

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 40 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 17

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 41 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 18

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 42 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 19

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 43 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 20

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 44 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 21

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 45 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 22

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 46 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 23

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 47 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 24

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 48 / 60

Stopping criterion satisfied here

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 25

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 49 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 26

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 50 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 27

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 51 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 28

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 52 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 29

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 53 / 60

Solutions as function of iteration index

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.05 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

k= 30

exact solution

• approx. sol.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 54 / 60

Example 2

Finite element discretization of L with variable coefficients on the ellipse Ω = {(x, y, z) | x2 + y2/4 ≤ 1, 0 ≤ z ≤ z1 = 0.6}. The stiffness matrix has dimension 8065 Data perturbation: ∼ 1%

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 55 / 60

Solution and Exact Data

−2 −1 1 2 −2 −1 1 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 −2 −1 1 2 −2 −1 1 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 56 / 60

Perturbed Data

−2 −1 1 2 −2 −1 1 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 57 / 60

Convergence history (0.5λc). Solution after 14 steps

5 10 15 20 25 10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10 step Difference True error Residual −2 −1 1 2 −2 −1 1 2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Step 14

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 58 / 60

Conclusions

3D Cauchy problem: complex 2D geometry + cylinder in z Krylov method + black box 2D elliptic solver Stability theory = ⇒ recipe for cut-off level Exponential of small matrix is computed (cheap) Safe-guarded stopping criterion: Solution difference (cheap) + Residual (rather expensive) Much fewer 2D elliptic solves than eigenvalue computation: 98 eigenvalues were smaller than the tolerance. MATLAB’s eigs: 300 2D solves Krylov: 25 solves

Highly accurate eigenvalues are not needed The data influence the basis (projection) vectors

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 59 / 60

Remaining work

Variable coefficients in z? Other Cauchy problems (parabolic, Helmholtz, transient electromagnetics)? Talk and paper are available on my web page.

Lars Eldén and Valeria Simoncini () Rational Krylov Method GAMM 2008 60 / 60