General linear methods

(1)

General linear methods

for ordinary differential equations

John Butcher

The University of Auckland, New Zealand

Acknowledgements:

Will Wright, Zdzisław Jackiewicz, Helmut Podhaisky

(2)

Solving ordinary differential equations numerically is, even today, still a great challenge.

This applies especially to stiff differential equations and to closely related problems involving algebraic

constraints (DAEs).

Although the problem seems to be solved — there are already highly efficient codes based on Runge–Kutta methods and linear multistep methods — questions concerning methods that lie between the traditional families are still open.

I will talk today about some aspects of these so-called

“General linear methods”.

(3)

Introduction to stiff problems

We will write a standard initial value problem in the form y^′(x) = f(x, y(x)), y(x₀) = y₀.

The variable x will be referred to as “time” even if it does not represent physical time.

Consider the problem





y₁^′ (x) y₂^′ (x) y₃^′ (x)



 =





y₂(x)

−y₁(x)

−Ly₁(x) + y₂(x) + Ly₃(x)



 ,

with initial values y₁(0) = 0, y₂(0) = 1, y₃(0) = ǫ.

(5)

If ǫ = 0, the solution is y(x) =





sin(x) cos(x)

sin(x)



 .

Change the value of ǫ and the third component changes to

sin(x) + ǫ exp(Lx).

We will see how this works for values of L = −25 and ǫ = 2.

The solution is plotted for the interval [0, 0.85].

For the approximate solution computed by the Euler

method, we choose n = 16 steps and then decrease n to 10.

(6)

0 1 2

0 0.5 x

y₁

y₂ y₃

Exact solution for L = −25

(7)

b b b

b b b b b b b b b b b b b

b b b b b b b b b b b

b b

b

b b

b b b b b b b b b b

0 1 2

0 0.5 x

y₁

y₂ y₃

Approximate solution for L = −25 with n = 16 steps

(8)

b

0 1 2

0 0.5 x

y1

y2

y³

L = −25, n = 10

General linear methods for ordinary differential equations – p. 8/46

(9)

Note that there seems to be no difficulty approximating y₁ and y₂, even for quite large stepsizes.

For stepsizes less than h = 0.08, the approximation to y₃ converges to the same value as y₁, but not in precisely the same way as for the exact solution.

For stepsizes greater than h = 0.08 the approximations to y₃ are hopeless.

We can understand what happens better by studying the difference y₃ − y₁.

(10)

The value of y₃ − y₁ satisfies the equation y^′(x) = Ly(x),

with solution const · exp(Lx).

The exact solution to y^′ = Ly is multiplied by exp(z) (where z = hL) in each time step but what about the computed solution?

According to the formula for the Euler method,

y_n = y_n−1 + hf(t_n−1, y_n−1), the numerical solution is multiplied instead by

1 + z in each time step.

(11)

The fact that 1 + z is a poor approximation to exp(z) when z is a very negative number, is at the heart of the difficulty in solving stiff problems.

Sometimes the “implicit Euler method” is used for the solution of stiff problems because in this case the

approximation exp(z) ≈ 1 + z is replaced by exp(z) ≈ (1 − z)⁻¹.

To see what happens in the case of the contrived problem we have considered, when Euler is replaced by implicit Euler, we present three further figures.

(12)

b

b b b

b b b b b b b b b b

b b

b

b b b b b b b b b b b b b b

0 1 2

0 0.5 x

y₁

y₂ y₃

Solution computed by implicit Euler for L = −25 with n = 16 steps

(13)

b

b b b b

b b

b

b b

b

b b

b

0 1 2

0 0.5 x

y₁

y₂ y₃

(14)

b

b b

b

b b

b

0 1 2

0 0.5 x

y₁

y₂ y₃

(15)

It looks as though, even for n as low as 5, we get completely acceptable performance.

The key to the reason for this is that |(1 − z)⁻¹| is bounded by 1 for all z in the left half-plane.

This means that a purely linear problem, with constant coefficients, will have stable numerical behaviour

whenever the same is true for the exact solution.

Methods with this property are said to be “A-stable”.

Although by no means a guarantee that it will solve all stiff problems adequately, A-stability is an important attribute and a useful certification.

(16)

Why do we want a new numerical method?

We want methods that are efficient for stiff problems Hence we cannot use highly implicit Runge-Kutta methods

Hence we want to avoid non- A-stable BDF methods We want methods that have high stage order

Otherwise order reduction can occur for stiff problems

Otherwise it is difficult to estimate accuracy Otherwise it is difficult to interpolate for dense output

(17)

We want methods for which asymptotically correct error estimates can be found

This means having accurate information available in the step to make this possible

High stage order makes this easier

We want methods for which changes to adjacent orders can be considered and evaluated

This means being able to evaluate, for example, h^p⁺²y⁽^p⁺²⁾ as well as h^p⁺¹y⁽^p⁺¹⁾

(p denotes the order of the method)

Evaluating the error for lower order alternative methods is usually easy

(18)

What are General Linear Methods?

General linear methods are the natural generalization of Runge–Kutta (multi-stage) methods and linear multistep (multivalue) methods.

This can be illustrated in a diagram.

General Linear Methods

Runge-Kutta Linear Multistep

Euler

(19)

We will consider r-value, s-stage methods, where r = 1 for a Runge–Kutta method

and

s = 1 for a linear multistep method.

Each step of the computation takes as input a certain

number (r) of items of data and generates for output the same number of items.

The output items correspond to the input items but advanced through one time step (h).

Within the step a certain number (s) of stages of computation are performed.

Denote by p the order of the method and by q the

“stage-order”.

(20)

At the start of step number n, denote the input items by y_i^[ⁿ^−1], i = 1, 2, . . . , r.

Denote the stages computed in the step and the stage derivatives by Yi and Fi respectively, i = 1, 2, . . . , s. For a compact notation, write

y^[ⁿ^−1] =







y₁^[n^−1]

y₂^[n−1]

...

yr^[n^−1]







, y^[ⁿ^] =







y₁^[n]

y₂^[n]

...

yr^[n]







, Y =





 Y₁ Y₂

...

Y_s







, F =





 F₁ F₂

...

F_s







(21)

The stages are computed by the formulae Yi =

s

X

j=1

aijhFj +

r

X

j=1

uijy_j^[n−1], i = 1, 2, . . . , s

and the output approximations by the formulae y_i^[ⁿ^] =

s

X

j=1

b_ijhF_j +

r

X

j=1

v_ijy_j^[ⁿ^−1], i = 1, 2, . . . , r

These can be written together using a partitioned matrix Y

y^[n]

=

A U B V

hF y^[n−1]

(22)

What we want from General Linear Methods

Convergence (consistency and stability) Local accuracy

Global accuracy (understand propagation of errors) A-stability for stiff problems

RK stability (behaves like Runge-Kutta method)

(23)

Propagation of errors

Let S denote a “starting method”, that is a mapping from R^N to R^rN and a corresponding finishing method

F : R^rN → R^N such that F ◦ S = id

The order of accuracy of a multivalue method is defined in terms of the diagram

E

S S

M O(h^p⁺¹)

(h = stepsize)

(24)

By duplicating this diagram over many steps, global error estimates may be found

E E E

S S S S S

M M M

O(h^p) F

O(h^p)

(25)

Algebraic analysis of order

We recall the way order is defined for Runge-Kutta methods.

We consider the solution of an autonomous differential equation y^′(x) = f(y(x)).

Let T denote the set of rooted trees

. . . To each of these is associated an “elementary

differential” F (t) defined in terms of the function f.

Also associated with each tree is its symmetry σ(t) and its density γ(t). In the following table, r(t) is the

“order” (number of vertices) of t.

(26)

r(t) t σ(t) γ(t) F (t) F (t)ⁱ

1 1 1 f fⁱ

2 1 2 f^′f f_jⁱf^j 3 2 3 f^′′(f, f) f_jkⁱ f^jf^k 3 1 6 f^′f^′f f_jⁱf_k^jf^k 4 6 4 f^′′′(f, f, f) f_jklⁱ f^jf^kf^l 4 1 8 f^′′(f, f^′f) f_jkⁱ f^jf_l^kf^l 4 2 12 f^′f^′′(f, f) f_jⁱf_kl^j f^kf^l 4 1 24 f^′f^′f^′f f_jⁱf_k^jf_l^kf^l

(27)

We will look at a single example tree in more detail.

t = r(t) = 6

σ(t) = 4

γ(t) = 18

6

1 1 3

1 1

F (t) = f^′′′(f,f,f^′′(f,f))

′′′

f f f^′′

f f

(28)

The solution to the standard autonomous initial-value problem

y^′(x) = f(y(x)), y(x₀) = y₀ has the formal Taylor series

y(x₀ + h) = y₀ + X

t∈T

h^r(t)

σ(t)γ(t)F (t)(y₀)

A Runge-Kutta method has the same series but with

1/γ(t) replaced by a sequence of “elementary weights”

characteristic of the specific Runge-Kutta method.

Expressions of this type are usually referred to as B-series.

(29)

Let G denote the set of mappings T ^# → R where T ^# = {∅} ∪ T (and ∅ is the empty tree).

Also write G₀ and G₁ as subsets for which ∅ 7→ 0 or

∅ 7→ 1 respectively.

A mapping G₁ × G → G can be defined which represents compositions of various operations.

In particular if this is restricted to G₁ × G₁ → G₁ it becomes a group operation.

Let E denote the mapping t 7→ 1/γ(t) and D the

mapping which takes every tree to zero except the unique order 1 tree, which maps to one.

(30)

A general linear method defined from the matrices [A, U, B, V ] has order p if there exists ξ ∈ G^r and η ∈ G^s₁, such that

η = AηD + U ξ,

Eξ = BηD + V ξ, (*)

where pre-multiplication by E and post-multiplication by D are scalar times vector products and (*) is to hold only up to trees of order p.

(31)

High order and stage order

If we consider only methods with q = p, then simpler criteria can be found.

Let exp(cz) denote the component-by-component

exponential then order p and stage-order q ≥ p − 1 is equivalent to

exp(cz) = zA exp(cz) + U w(z) + O(z^q+1), exp(z)w(z) = zB exp(cz) + V w(z) + O(z^p+1), where the power series

w(z) = w₀ + w₁z + · · · + w_pz^p represents an input approximation

y^[0] =

p

X

i=0

wihⁱy⁽ⁱ⁾(x₀).

(32)

GL Methods and RK stability

A General Linear Method is “Runge–Kutta stable” if its stability matrix

M(z) = V + zB(I − zA)⁻¹U has only one non-zero eigenvalue.

Armed with this property we should expect attraction to the manifold S(R^N), and stable adherence to this

manifold.

This means that the method starts acting like a one-step method.

(33)

Doubly companion matrices

Matrices like the following are “companion matrices” for the polynomial

zⁿ + α₁zⁿ⁻¹ + · · · + α_n or

zⁿ + β₁zⁿ⁻¹ + · · · + βn, respectively:







−α₁−α₂−α₃· · · −α_n−1−α_n 1 0 0 · · · 0 0 0 1 0 · · · 0 0 ... ... ... ... ...

0 0 0 · · · 0 0 0 0 0 · · · 1 0





 ,







0 0 0 · · · 0 −β_n 1 0 0 · · · 0 −β_n−1 0 1 0 · · · 0 −β_n−2

... ... ... ... ...

0 0 0 · · · 0 −β₂ 0 0 0 · · · 1 −β₁







(34)

Their characteristic polynomials can be found from det(I − zA) = α(z) or β(z), respectively, where,

α(z) = 1+α₁z+· · ·+αnzⁿ, β(z) = 1+β₁z+· · ·+βnzⁿ. A matrix with both α and β terms:

X =







−α₁ −α₂ −α₃ · · · −αn−1 −αn − βn

1 0 0 · · · 0 −β_n−1

0 1 0 · · · 0 −βn−2

... ... ... ... ...

0 0 0 · · · 0 −β₂

0 0 0 · · · 1 −β₁





 ,

is known as a “doubly companion matrix” and has characteristic polynomial defined by

det(I − zX) = α(z)β(z) + O(zⁿ⁺¹)

(35)

Matrices Ψ⁻¹ and Ψ transforming X to Jordan canonical form are known.

In the special case of a single Jordan block with n-fold eigenvalue λ, we have

Ψ⁻¹ =







1 λ + α₁ λ² + α₁λ + α₂ · · · 0 1 2λ + α₁ · · ·

0 0 1 · · ·

... ... ... . . .





 ,

where row number i + 1 is formed from row number i by differentiating with respect to λ and dividing by i.

We have a similar expression for Ψ:

(36)

Ψ =







. . . ... ... ...

· · · 1 2λ + β₁ λ² + β₁λ + β₂

· · · 0 1 λ + β₁

· · · 0 0 1







The Jordan form is Ψ⁻¹XΨ = J + λI, where Jij =δi,j+1. That is

Ψ⁻¹XΨ =







λ 0 · · · 0 0 1 λ · · · 0 0 ... ... ... ...

0 0 · · · λ 0 0 0 · · · 1 λ







(37)

Inherent Runge–Kutta stability

Using doubly companion matrices, it is possible to construct GL methods possessing RK stability using

rational operations. The methods constructed in this way are said to possess “Inherent Runge–Kutta Stability”.

Apart from exceptional cases, (in which certain matrices are singular), we characterize the method with

r = s = p + 1 = q + 1 by the parameters λ single eigenvalue of A

c₁, c₂, . . . , c_s stage abscissae Error constant

β₁, β₂, . . . , β_p elements in last column of s × s doubly companion matrix X

Information on the structure of V

(38)

Consider only methods for which the step n outputs approximate the “Nordsieck vector”:







y₁^[n]

y₂^[ⁿ^] y₃^[n]

...

y_p^[ⁿ₊₁^]







≈







y(xn) hy^′(x_n) h²y^′′(xn)

...

h^py^(p)(x_n)







For such methods, V has the form V =

1 v^T 0 V˙

(39)

Such a method has the IRKS property if a doubly

companion matrix X exists so that for some vector ξ, BA = XB, BU = XV − V X + e₁ξ^T , ρ( ˙V ) = 0 It will be shown that, for such methods, the stability

matrix satisfies

M(z) ∼ V + ze₁ξ^T (I − zX)⁻¹

which has all except one of its eigenvalues zero. The non-zero eigenvalue has the role of stability function

R(z) = N(z) (1 − λz)^s

(40)

Proof of stability result

(I − zX)M(z)(I − zX)⁻¹

= (V − zXV + z(B − z

XB = BA

XB)(I − zA)⁻¹U)(I − zX)⁻¹

= (V − zXV + zBU)(I − zX)⁻¹

BU = XV − V X + e₁ξ^T

= V + ze₁ξ^T (I − zX)⁻¹

This has a single non-zero eigenvalue equal to

R(z) = 1 + zξ^T (I − zX)⁻¹e₁ = det(I + z(e₁ξ^T − X)) det(I − zX)

(41)

Example methods

The following third order method is explicit and suitable for the solution of non-stiff problems

"

AU BV

#

=







0 0 0 0 1 ¹₄ ₃₂¹ ₃₈₄¹

−₁₈₈₅¹⁷⁶ 0 0 0 1 ²²³⁷₃₇₇₀ ₁₅₀₈₀²²³⁷ ₉₀₄₈₀²¹⁴⁹

−³³⁵⁶²⁴₃₁₁₀₂₅ ²⁹₅₅ 0 0 1 ^1619591_1244100 ²⁶⁰⁰²⁷₉₀₄₈₀₀ _39811200^1517801

−⁶⁷⁸⁴³₆₄₃₅ ³⁹⁵₃₃ −5 0 1 ²⁹⁴²⁸₆₄₃₅ ⁵²⁷₅₈₅ ₁₀₂₉₆₀⁴¹⁸¹⁹

0 0 0 1 0 0 0 0

82

33 −²⁷⁴₁₁ ¹⁷⁰₉ −⁴₃ 0 ⁴⁸²₉₉ 0 −¹⁶¹₂₆₄

−8 −12 ⁴⁰₃ −2 0 ²⁶₃ 0 0







(42)

The following fourth order method is implicit, L-stable, and suitable for the solution of stiff problems







1

4 0 0 0 0 1 ³₄ ¹₂ ¹₄ 0

−₅₄₂₇₂⁵¹³ ¹₄ 0 0 0 1 ²⁷⁶⁴⁹₅₄₂₇₂ ₂₇₁₃₆⁵⁶⁰¹ ₅₄₂₇₂¹⁵³⁹ −₆₇₈₄⁴⁵⁹

3706119

69088256 −₃₈₁₉⁴⁸⁸ ¹₄ 0 0 1 _207264768^15366379 _34544128⁷⁵⁶⁰⁵⁷ _69088256^1620299 −₄₅₄₅₂₈⁴⁸⁵⁴

32161061

197549232 −¹¹¹⁸¹⁴₂₃₂₉₅₉ ¹³⁴₁₈₃ ¹₄ 0 1−_197549232^32609017 _32924872⁹²⁹⁷⁵³ _32924872^4008881 _3465776¹⁷⁴⁹⁸¹

−_2948496¹³⁵⁴²⁵ −₁₀₄₃₁⁶⁴¹ ₁₈₃⁷³ ¹₂ ¹₄ 1 −_8845488³⁶⁷³¹³ −_1474248²²⁷²⁷ ₉₈₂₈₃₂⁴⁰⁹⁷⁹ ₂₅₈₆₄³²³

0 0 0 0 1 0 0 0 0 0

2255

2318 −⁴⁷¹²⁵₂₀₈₆₂ ⁴⁴⁷₁₂₂ −¹¹₄ ⁴₃ 0 −²⁸⁷⁴⁵₂₀₈₆₂ −₁₃₉₀₈¹⁹³⁷ ₁₈₅₄₄³⁵¹ ₉₇₆⁶⁵

12620

10431 −⁹⁶³⁸⁸₃₁₂₉₃ ³³⁶⁴₅₄₉ −¹⁰₃ ⁴₃ 0 −⁷⁰⁶³⁴₃₁₂₉₃ −₁₀₄₃₁²⁰⁵⁰ −₂₃₁₈¹⁸⁷ ¹¹³₃₆₆

414

1159 −²⁹⁹⁵⁴₃₁₂₉₃ ¹³⁰₆₁ −1 ¹₃ 0 −²⁷⁰⁵²₃₁₂₉₃ −₁₀₄₃₁¹¹³ −₄₆₃₆⁴⁹¹ ¹⁶¹₇₃₂







(43)

Implementation questions

Many implementation questions are similar to those for traditional methods but there are some new challenges.

We want variable order and stepsize and it is even a realistic aim to change between stiff and non-stiff methods automatically.

Because of the variable order and stepsize aims, we wish to be able to do the following:

Estimate the local truncation error of the current step Estimate the local truncation error of an alternative method of higher order

Change the stepsize with little cost and with little impact on stability

(44)

My final comment is that all these things are possible and many of the details are already known and understood.

What is not known is how to choose between different choices of the free parameters.

Identifying the best methods is the first step in the

construction of a competitive stiff and non-stiff solver.

(45)

References

The following recent publications contain references to earlier work:

Numerical Methods for Ordinary Differential Equations, J. C. Butcher, J. Wiley, Chichester, (2003)

General linear methods, J. C. Butcher, Acta Numerica 15 (2006), 157–256.

(46)

Gracias Thank you

Obrigado

General linear methods