Jan Verwer

1
On the Numerical Integration of
Maxwells Equations

• Jan Verwer

Centrum voor Wiskunde en Informatica
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
• Maxwells equations with a damping (conduction)
  term
• Aim methods of higher order (beyond standard
  two)

• Ill 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
Maxwells equations
sE 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 a 1

Note that t 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

a 1 and t , h decrease simultaneously

CO2
RK4
CO4
26
Accuracy versus work at t 0.5

a 1 and t , h decrease simultaneously Order
reduction for RK4 and CO4 due to time- dependent
BCs

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 gt 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 t-expansion
• Results for the 4th-order extrapolation
• (3 times more expensive than CO2)

30
Accuracy versus work at t 0.5

a 108 and t , h decrease simultaneously Same
order behaviour for non-stiff case (small a)

CO2
CO4E
31
Accuracy versus work at t 0.5

a 108 and t , h decrease simultaneously Same
order behaviour for non-stiff case (small a)

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 10p (zero bndry) and L
p (nonzero bndry) (using 4th-order in space)
34
Accuracy versus work at t p for L 10p


t , h decrease simultaneously
CO2
no order reduction
CO4E
4
35
Accuracy versus work at t p for L p


t , h decrease simultaneously
CO2
tiny order reduction
3
CO4E
4
36
Another damped wave equation systemcoupled
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 Maxwells 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)