Taming explosive computational instability: Com- pensated explicit - PowerPoint PPT Presentation

Taming explosive computational instability: Com- pensated explicit time-marching schemes in multi- dimensional, nonlinear, well-posed or ill-posed ini- tial value problems for partial differential equations Alfred S. Carasso Applied and Computational Mathematics Division NIST, Gaithersburg, Maryland. OCTOBER 29 2015. NIST Internal Reports # 8027, 8079. Int J Geomath 5 (2014), pp. 1-16. Additional papers to appear in Inverse Problems in Sci- ence and Engineering

Consider well-posed forward problem w t = Lw, t > 0 , with L a 2nd order elliptic differential operator. w n +1 = w n + ∆ tLw n , Pure Explicit ideal. w n +1 = w n − 1 + 2∆ tLw n , Leapfrog better. 2 w n +1 + ∆ tLw n +1 = 2 w n + ∆ tLw n . Crank-Nicholson Pure Explicit ⇒ impractically small ∆ t (Courant). Leapfrog always unstable no matter how small ∆ t . Crank-Nicholson ⇒ large algebraic problem at each n . Stabilize explicit schemes. Apply to difficult multidi- mensional nonlinear well-posed forward problems. Avoid implicit scheme algebraic computations at each n . Stabilized schemes can also be run backward in time. Solve previously intractable ill-posed problems.

EXAMPLES OF IRREVERSIBLE PDE PROBLEMS Advection dispersion equation w t = σ ∆ w − ∇ . ( β w ) Wave propagation in viscous fluid w tt = c 2 ∆ w + γ ∆ w t Coupled sound and heat flow w tt = c 2 ∆ w − c 2 ∆ u, u t = σ ∆ u − λw t Nonlinear variations. Other coupled systems.

ENVIRONMENTAL FORENSICS Groundwater Pollution Identify source of groundwater contamination by solving advection dispersion eqn backward in time .

IMAGES MAKE GOOD TEST PROBLEMS Shirley Temple image Plot of intensity values

Consider PDE problem: w t = Lu, 0 ≤ t ≤ T, w (0) = f. Assume L linear selfadjoint second order elliptic op- erator, with variable coefficients independent of t . Eigenvalues { λ m } ∞ m =1 ; 0 < | λ 1 | ≤ | λ 2 | ≤ · · · ≤ | λ m | ≤↑ ∞ PDE well-posed if all λ m < 0 , ill-posed if all λ m > 0. For any q > 0, operator ( L 2 ) q/ 2 has eigenvalues | λ m | q . For small ω > 0 , q ≥ 2, define smoothing operator S = exp {− ω ∆ t ( L 2 ) q/ 2 } . Define β q = ω 1 / (1 − q ) . Operator S can be synthesized in terms of eigenpairs { λ m , φ m } of L, assumed known or precomputed . Consider compensated pure explicit scheme w n +1 = Sw n + ∆ tSLw n , w 0 = f . n ≥ 0 ,

S = exp {− ω ∆ t ( L 2 ) q/ 2 } , β q = ω 1 / (1 − q ) , w n +1 = Sw n + ∆ tSLw n , w 0 = f . n ≥ 0 , Theorem 1: Compensated pure explicit scheme always stable and � w n +1 � 2 ≤ (1 + β q ∆ t ) � w n � 2 , if n ≥ 0. Hence, if t n = n ∆ t, � w ( t n ) � 2 ≤ exp( β q t n ) � f � 2 . w n +1 = w n − 1 + 2∆ tLw n , n ≥ 1 . Leapfrog scheme Put v n = w n − 1 , and consider equivalent system v n +1 = w n , w n +1 = v n + 2∆ tLw n , n ≥ 1 . Define U n = [ v n , w n ] , � U n � 2 2 ≡� v n � 2 2 + � w n � 2 2 With ω, q, β q , S, as above, compensated leapfrog , v n +1 = Sw n , w n +1 = Sv n + 2∆ tSLw n , n ≥ 1 . Compensated leapfrog scheme always Theorem 2: stable , and � U n +1 � 2 ≤ (1 + 2 β q ∆ t ) � U n � 2 , if n ≥ 1. Hence, if t n = n ∆ t, � U ( t n ) � 2 ≤ exp(2 β q t n ) � U ( t 1 ) � 2 .

Remark. Theorems valid even if w t = Lw ill-posed !!! Replace inconvenient smoothing operator S with S ∆ = exp {− ǫ ∆ t ( − ∆) p } . Laplacian-based Known characteristic pairs of ∆ in several geometries. Rectangular domain ⇒ fast FFT synthesis of S ∆ . � S ∆ g � 2 ≤� Sg � 2 . More pre- Major Assumption: cisely, given ω > 0 , and q > 2 , ∃ ǫ > 0 and real ∋ ∀ g ∈ L 2 and sufficiently small ∆ t > 0 , p ≥ q, � exp {− ǫ ∆ t ( − ∆) p } g � 2 ≤� exp {− ω ∆ t ( L 2 ) q/ 2 } g � 2 . Above inequality not proved . Appears verified in many interesting computational examples. Related to Gaus- sian lower bounds for heat kernels . Deep theory. Theorems 1 and 2 remain valid with S ∆ replacing S .

Error bounds for well-posed w t = Lw, 0 ≤ t ≤ T . Given exact initial data f . Schemes stable , but do not converge as ∆ t ↓ 0. Smoothing ⇒ residual error . Let w ex ( t ) be exact solution of w t = Lw, 0 ≤ t ≤ T . Let e ( t n ) = w ex ( t ) − w ( t n ) be the error at t n = n ∆ t . Let B = sup 0 ≤ t ≤ T {� ( − ∆) p w ex ( t ) � 2 } . For compensated pure explicit well-posed case � e ( t n ) � 2 ≤ ( Bǫ/ω )( e β q t n − 1) / ( β q ) q + O (∆ t ) ( Bǫ/ω )( e β q t n − 1) / ( β q ) q ≡ Stabilization Penalty . Involves parameters ǫ, p, ω, q For compensated leapfrog well-posed case � e ( t n ) � 2 ≤ ( Bǫ/ω )( e 2 β q t n − 1) / 2( β q ) q + O (∆ t 2 ) (larger stabilization penalty).

( Bǫ/ω )( e β q t n − 1) / ( β q ) q ≡ Explicit stabilization penalty . If ω = 1 . 0 E − 8 , q = 2 . 75 , T = 1 . 0 E − 4 , Example 1 then β q = 37276 , ( e β q T − 1) = 41 , ( β q ) q = 3 . 73 E 12 , ⇒ Explicit stabilization penalty =( Bǫ/ω ) × 1 . 1 E − 11. Error bounds for ill-posed w t = Lw, 0 ≤ t ≤ T . Noisy initial data f δ , not exact f , with � f δ − f � 2 ≤ δ . Let e ( t n ) = w ex ( t n ) − w n , denote the error at t n = n ∆ t . With B = sup t {� ( − ∆) p w ex ( t ) � 2 } , Pure Explicit error � e ( t n ) � 2 ≤ δe β q t n + ( Bǫ/ω )( e β q t n − 1) / ( β q ) q + O (∆ t ) With prescribed L 2 bound � w ex ( T ) � 2 ≤ M , choosing β q = (1 /T ) log( M/δ ) , ⇒ quasi optimal error bound: � e ( t n ) � 2 ≤ M t n /T δ ( T − t n ) /T + O (∆ t ) + stabilization penalty. Quasi optimal error also in compensated leapfrog .

Linear autonomous selfadjoint analysis Compensated pure explicit and leapfrog can perform well on limited but significant class of problems with small T , and small stabilization penalty. Stablity in time-reversed problem ⇒ Quasi-optimal results, within fundamental uncertainty M ( T − t ) /T δ t/T . Stabilizing pair ( ǫ, p ) for S ∆ located interactively . Nonlinear problems Explicit schemes extremely advantageous in multi- dimensional nonlinear problems on fine meshes. Laplacian smoother S ∆ effective with nonlinear L ??? Small stabilization penalty in nonlinear case ??? Stability in nonlinear time-reversed problems ???

FORWARD LEAPFROG COMPUTATION w t = exp( σw ) ∇ . { q ( x, y, t ) ∇ w } + c ( w ) w x + d ( w ) w y STABILIZED LEAPFROG FORWARD TIME MARCHING IN LARGE ARRAY Sharp 1024x1024 USAF chart image Nonlinear parabolic leapfrog blur Painless Leapfrog O (∆ t 2 ) computation. Crank-Nicholson ⇒ solve nonlinear system of order 10 6 at each n .

TIME-REVERSED LEAPFROG COMPUTATION w t = exp( σw ) ∇ . { q ( x, y, t ) ∇ w } + c ( w ) w x + d ( w ) w y NONLINEAR LEAPFROG EXPERIMENT DONE ON SEPT 3 2015 SUCCESSFUL STABILIZED LEAPFROG FORWARD AND BACKWARD TIME MARCHING IN NONLINEAR PARABOLIC EQUATION WITH VARIABLE TIME DEPENDENT COEFFICIENTS LEAPFROG NONLINEAR BLUR LEAPFROG NONLINEAR DEBLUR SHARP JOAN CRAWFORD IMAGE Moderate nonlinearity allows full leapfrog backward re- covery, from t = T to t = 0.

ONLY PARTIAL LEAPFROG RECOVERY w t = exp( σw ) ∇ . { q ( x, y, t ) ∇ w } + c ( w ) w x + d ( w ) w y Leapfrog nonlinear parabolic blurring of sharp 512x512 Joan Crawford image, with partial Leapfrog deblurring. (Experiments performed on Sept 18 2015) 30% Partial Leapfrog Deblur Nonlinear Parabolic Leapfrog Blur Stronger nonlinearity only allows partial leapfrog back- ward recovery, from t = T to t = 0 . 7 T .

LEAPFROG RECOVERY MAY FAIL AS t ↓ 0 . w t = exp( σw ) ∇ . { q ( x, y, t ) ∇ w } + c ( w ) w x + d ( w ) w y STABILIZED LEAPFROG SCHEME MARCHED BACKWARD IN TIME Nonlinear parabolic blurring and deblurring of 512x512 image INPUT BLURRED IMAGE t=T PARTIAL DEBLURRING t= 0.7 T PARTIAL DEBLURRING t= 0.6 T PARTIAL DEBLURRING t= 0.5 T PARTIAL DEBLURRING t= 0.4 T PARTIAL DEBLURRING t=0.2 T

Larger nonlinear uncertainty = M 1 − µ ( t ) δ µ ( t ) , where µ ( t ) ↓ 0 exponentially as t decreases from t = T . MAY NOT PERMIT FULL RECONSTRUCTION. Behavior of Holder exponent in backward problems Autonomous, linear self adjoint problem Nonlinear problem

PURE EXPLICIT SCHEME SURPRISE! RECOVERS FROM BLURRED LEAPFROG DATA w t = exp( σw ) ∇ . { q ( x, y, t ) ∇ w } + c ( w ) w x + d ( w ) w y HEAVIER NONLINEAR LEAPFROG/EXPLICIT EXPERIMENT DONE ON SEPT 1 2015 LEAPFROG FORWARD TIME MARCHING FOLLOWED BY PURE EXPLICIT BACKWARD TIME MARCHING IN HEAVIER NONLINEAR PARABOLIC EQUATION AFTER EXPLICIT NONLINEAR DEBLUR SHARP GENE TIERNEY IMAGE AFTER LEAPFROG NONLINEAR BLUR

REMARKABLE DATA SURFACE RECOVERY w t = exp( σw ) ∇ . { q ( x, y, t ) ∇ w } + c ( w ) w x + d ( w ) w y HEAVIER NONLINEAR LEAPFROG/EXPLICIT EXPERIMENT DONE ON SEPT 1 2015 LEAPFROG FORWARD TIME MARCHING FOLLOWED BY PURE EXPLICIT BACKWARD TIME MARCHING IN HEAVIER NONLINEAR PARABOLIC EQUATION SHARP GENE TIERNEY DATA SURFACE AFTER LEAPFROG NONLINEAR BLUR AFTER EXPLICIT NONLINEAR DEBLUR

EXPLICIT RECOVERY FROM LEAPFROG . w t = exp( σw ) ∇ . { q ( x, y, t ) ∇ w } + c ( w ) w x + d ( w ) w y HEAVY NONLINEAR PDE BLURRING NONLINEAR PARTIAL DEBLURRING USING STABILIZED LEAPFROG SCHEME USING STABILIZED EXPLICIT SCHEME