Jan Verwer

1. On the Numerical Integration of Maxwell’s Equations Centrum voor Wiskunde en Informatica Jan Verwer Research in progress with Mike Botchev, TU-Twente Workshop ‘Time Integration of Evolution Equations’ Innsbruck, Sept. 12 – 15, 2007

2. Contents • The lecture is about numerical integration of • Maxwell’s equations with a damping (conduction) term • Aim: methods of higher order (beyond standard two) • I’ll discuss • The Maxwell equations and the semi-discrete system • Second-order methods (rehearsal of reference methods) • How to obtain higher order? • Recently started joint work • Related results for another damped wave equation, viz the coupled sound and heat flow problem

3. The Maxwell equationsand the semi-discrete system

4. Maxwell’s equations – σE is the conduction term 3D, so six time-dependent PDEs

5. General semi-discrete system - Curl matrix K may not be square - Mass matrices are spd - Conduction matrix S is also spd

6. Stability test model Gives the (homoge- geneous) sd-system Which can be reduced to the 2x2 stability test model Numerical stability for this test model is “equivalent” to L2 – stability for the semi-discrete system

7. Rehearsal of three second-orderintegration methods combining theleapfrog - and trapezoidal rule

8. Implicit trapezoidal rule ITR Mimics stability However, at the cost of dealing with this huge saddle-type matrix Under investigation by Mike B.

9. To avoid the saddle-type system we will treat the wave terms explicitly and the conduction term implicitly (IMEX approach) (i) Two-step leapfrog-trapezoidal rule In terms of b and e, this gives the following TW rule S fully explicit like K is no option!

10. (ii) Staggered leapfrog-trapezoidal rule Using time-staggering gives the ST rule S = 0, Yee-scheme Implicit in S -- explicit in K, and gives the much simpler block-diagonal matrix instead of Nice matrix for Krylov solvers

11. (iii) Composite leapfrog-trapezoidal rule Eliminating through gives

12. (iii) Composite leapfrog-trapezoidal rule I call this the CO rule as it is a (symmetric) COmposition scheme:

13. (iii) Composite leapfrog-trapezoidal rule • CO rule • is symmetric by composition and is a common one-step • scheme in the sense that it steps from • costs the same as the staggered ST rule since the third-stage • derivative evaluation can be reused for the next step • NB. ST is “identical” to CO only with the right initialization. • Otherwise the two methods differ.

14. Stability for the ST rule The test model gives Characteristic eq. The root condition is satisfied Hence unconditional stability for the conduction term

15. Summary of twostep TW, staggered ST, and composite CO rules All three • are of second order • cost the same per time step • treat the curl terms explicitly and the conduction term implicitly • are unconditionally stable for the conduction term and • conditionally stable for the curl term for the test model • and ST identical to CO with right initialization Recall that

16. Numerical illustration (1D) Numerical results for t = 0.1 (zero bndry) and t = 0.5 (nonzero bndry)

17. Fourth-order compact scheme in space Dirichlet boundary values prescribed from exact solution. This discretization fits in the general semi-discrete Maxwell form For stability analysis the test model applies which gives the step size restrictions

18. Efficiency plot: accuracy versus work, t = 0.5 Here α = 1 Note that τ and h decrease simultaneously We see here the temporal, second-order convergence TW CO ST

19. How to design higher-order methods,which, like the LF-TR methods, treat the wave terms explicitly?

20. First for moderate (non-stiff) conduction High-order, standard explicit RK-methods make a good choice, since we treat the wave terms explicitly anyhow ! Perform stability check for our test model through the stability region

21. Stability region of RK4

22. First for moderate (non-stiff) conduction A second possibility are high-order composition methods For the Maxwell system the resulting scheme is again explicit in K and implicit in S

23. Stability region of CO4 CO4: 4th-order scheme, s = 5, coefs. from McLachlan ‘95

24. Comparison: RK4 - CO2 - CO4 Fourth-order compact in space gives the step size restrictions

25. Accuracy versus work at t = 0.1 α = 1 and τ,h decrease simultaneously CO2 RK4 CO4

26. Accuracy versus work at t = 0.5 α = 1 and τ,h decrease simultaneously Order reduction for RK4 and CO4 due to time- dependent BC’s CO2 RK4 CO4

27. Next, for large (stiff) conduction Stiffness requires an implicit treatment of the conduction term, while we wish to keep handling the wave terms explicitly • High-order composition is now impossible due to the negative substeps →instability

28. Stability region of CO4 CO4: 4th-order scheme, s = 5, coefs. from McLachlan ‘95 hole Hole is due to negative substep (for any CO of order > 2)

29. Next, for large (stiff) conduction Stiffness requires an implicit treatment of the conduction term, while we wish to keep handling the wave terms explicitly • High-order composition is now impossible due to the negative substeps →instability • Alternative: global R-extrapolation (output only) of CO2 • -- CO2 determines stability • -- Exploits even τ-expansion • Results for the 4th-order extrapolation (3 times more expensive than CO2)

30. Accuracy versus work at t = 0.5 α= 10^8 and τ , h decrease simultaneously Same order behaviour for non-stiff case (small α) CO2 CO4E

31. Accuracy versus work at t = 0.5 α= 10^8 and τ , h decrease simultaneously Same order behaviour for non-stiff case (small α) CO2 CO4E Extrapolation seems free of order reduction?

32. Error analysis results on order reduction Thm. For the linear system Covers Maxwell 2nd-order CO2 is free of h-dependent order reduction and thus show common ODE convergence. This also holds under extrapolation while the order extrapolates to 4, 6, 8, …. NB. CO2 enjoys error cancellation. With nonlinearity this might not happen and reduction may occur under extrapolation ?

33. Is extrapolation in general free of order reduction? Counter example for CO4E: the Sine-Gordon equation formulated in the system form Numerical results for L = 10π (zero bndry) and L = π (nonzero bndry) (using 4th-order in space)

34. Accuracy versus work at t =πfor L = 10π τ , h decrease simultaneously CO2 no order reduction CO4E 4

35. Accuracy versus work at t =πfor L = π τ , h decrease simultaneously CO2 tiny order reduction 3 CO4E 4

36. Another damped wave equation system:coupled sound and heat flow

37. Coupled sound and heat flow Scaled linearized equations expressing conservation of energy, mass and momentum, cf. Richtmyer & Morton, Sect. 10.4: Write and define the 1st-order Euler-type scheme

38. Coupled sound and heat flow The composition then defines the symmetric, 2nd-order, one-step integration scheme • effectively 3 stages per step • explicit in velocity and volume, implicit in energy • no stability restrictions from energy eq.; stab. is determined by wave eq. part (d = 0, F-analysis for semi-discrete system) • as before we can apply global R-extrapolation for higher order

39. Numerical illustration 2D test problem from Sommeijer & V., JCP ’07 Unit square, periodicity, 4-th order central in space, h = 1/100 CO2 Temporal errors only CO4E

40. Summary • Time integration of Maxwell’s equations containing a conduction term for the electric field • Wave (curl) terms explicit (if allowed) to avoid solving expensive saddle-type matrix systems • For nonstiff conduction, high order composition (beyond 2) is very efficient (but is not stable for stiff conduction) • For stiff conduction, high order extrapolation of the 2nd-order CO2 rule is very efficient • Extrapolation of a similar 2nd-order CO2 rule is also efficient for similar damped wave equations (coupled sound-heat flow)