Validated Methods for Initial Value Problems for Ordinary - - PowerPoint PPT Presentation

Validated Methods for Initial Value Problems for Ordinary Differential Equations

Ken Jackson

Computer Science Department University of Toronto

Fields Institute Workshop on Hybrid Methodologies for Symbolic-Numeric Computation 19 November 2011

Principal Contributors

Ned Nedialkov, McMaster University Wayne Hayes, University of California at Irvine Markus Neher, Universit¨ at Karlsruhe George Corliss, Marquette University John Pryce, Cardiff University

1

Outline

1

Preliminaries Standard Versus Validated Numerical Methods Interval Arithmetic Automatic Differentiation

2

Validated Numerical Methods for IVPs for ODEs Overview Taylor Series Methods for IVPs for ODEs The Wrapping Effect Controlling the Wrapping Effect Analyzing the Methods for Controlling the Wrapping Effect

2

Standard Versus Validated Numerical Methods

Standard numerical methods compute approximations to a mathematically correct result. Validated numerical methods (also called interval methods) compute bounds that are guaranteed to contain the mathematically correct result. Validated methods use interval arithmetic to bound rounding errors errors from numerical approximations errors from the model Middle ground between standard numerical methods and symbolic computation.

3

Interval Arithmetic (Moore 1966)

An interval [a] is defined by two endpoints a ¯ ≤ ¯ a: [a] = [a ¯, ¯ a] = {x | a ¯ ≤ x ≤ ¯ a} . Example: 1/3 = 0.3333333... ∈ [0.333, 0.334]. Interval-arithmetic operations: If [a] and [b] are intervals and ∗ ∈ {+, −, ×, /}, then [a] ∗ [b] = {x ∗ y | x ∈ [a] , y ∈ [b]} , 0 / ∈ [b] when ∗ = /, which can be written in the equivalent form: [a] + [b] =

• a

¯ + b ¯, ¯ a + ¯ b

• [a] − [b]

=

• a

¯ − ¯ b, ¯ a − b ¯

• [a] [b]

=

• min
• a

¯b ¯, a ¯ ¯ b, ¯ ab ¯, ¯ a¯ b

• , max
• a

¯b ¯, a ¯ ¯ b, ¯ ab ¯, ¯ a¯ b

• [a] / [b]

= [a ¯, ¯ a]

• 1/¯

b, 1/b ¯

• ,

0 / ∈ [b] .

4

Some Properties of Interval Arithmetic

Interval arithmetic is Inclusion Monotone: If [a1] ⊆ [a2] and [b1] ⊆ [b2], then [a1] ∗ [b1] ⊆ [a2] ∗ [b2] , ∗ ∈ {+, −, ×, /} . The dependency problem: a source of overestimations [x] − [x] = 0 The distributive law does not hold in general: another source of overestimations Instead, we have sub-distributivity [a] ([b] + [c]) ⊆ [a] [b] + [a] [c] with equality in some special cases only. For example, if a ∈ R, then a([b] + [c]) = a [b] + a [c] .

5

Some Extensions

Interval Vectors and Matrices are vectors and matrices with interval components. the arithmetic operations are defined by the same formulas as in the scalar case, except that scalars are replaced by intervals. [A] ⊆ [B] ⇐ ⇒ [aij] ⊆ [bij] ∀i, j Some Definitions width: w ([a]) = ¯ a − a ¯ = ⇒ w ([a] + [b]) = w ([a]) + w ([b]) midpoint: m([a]) = (¯ a + a ¯)/2 Width and midpoint apply componentwise to vectors and matrices. If T is a point matrix and [a] is an interval vector, then w (T [a]) = |T|w ([a])

6

Interval-Valued Functions

Range of f over an interval vector [a] is R (f ; [a]) = {f (x) | x ∈ [a]} . Interval-arithmetic evaluation of f on [a], f ([a]), is obtained by replacing reals with intervals. Thus, R (f ; [a]) ⊆ f ([a]). Mean-value form: If f : Rn → R is continuously differentiable on D ⊆ Rn, then for [a] ⊆ D and any y and ˆ y ∈ [a] f (y) = f (ˆ y) + fy(η)(y − ˆ y) for some η ∈ [a]. Therefore R (f ; [a]) ⊆ fM([a] , ˆ y) ≡ f (ˆ y) + fy([a])([a] − ˆ y). Often fM([a] , ˆ y) is a better enclosure of R (f ; [a]) than f ([a]).

7

Example

Consider f (x) = x − x2 for x ∈ [a] = 1

4, 3 4

• .

The range of f (x) is R (f ; [a]) = 1

16 [3, 4].

w 1

16 [3, 4]

• = 1

16.

For f (x) = x − x2, interval-arithmetic evaluation gives: f ([ 1

4, 3 4]) = [ 1 4, 3 4] − [ 1 4, 3 4] × [ 1 4, 3 4] = 1 16[−5, 11]

w 1

16[−5, 11]

• = 1.

For f (x) = x(1 − x), interval-arithmetic evaluation gives: f ([ 1

4, 3 4]) = [ 1 4, 3 4] × (1 − [ 1 4, 3 4]) = 1 16[1, 9]

w 1

16[1, 9]

• = 1

2.

Mean-value evaluation gives: fM([ 1

4, 3 4], 1 2) = f ( 1 2) + f ′([ 1 4, 3 4]) × ([ 1 4, 3 4] − 1 2) = 1 16[2, 6]

w 1

16[2, 6]

• = 1

4.

8

A Bigger Example: Newton’s Method

Assume f : R → R. Newton’s method for finding a root of f is xn+1 = xn − f (xn)/f ′(xn) Replacing points by intervals, we get [xn+1] = [xn] − f ([xn])/f ′([xn]) But this doesn’t work, since w ([xn+1]) ≥ w ([xn]).

9

A Better Approach: Interval Newton’s Method

Suppose f (ˆ x) = 0 and ˆ x ∈ [xn]. Let xn be some point in [xn]. (Often a good choice is xn = m([xn]).) Then, for some η ∈ [xn], 0 = f (ˆ x) = f (xn) + f ′(η)(ˆ x − xn) Therefore ˆ x = xn − f (xn)/f ′(η) ˆ x ∈ xn − f (xn)/f ′([xn]) ˆ x ∈ [xn] ∩

• xn − f (xn)/f ′([xn])
• Interval Newton’s Method:

[xn+1] = [xn] ∩

• xn − f (xn)/f ′([xn])
• Note ˆ

x ∈ [xn+1] ⊆ [xn]. Under suitable conditions, w ([xn]) decreases quadratically with n.

10

Example

Find the positive root of f (x) = x2 − 2. The solution is √ 2 = 1.4142135.... For this example, Interval Newton’s Method is [xn+1] = [xn] ∩

• xn − (x2

n − 2)/(2[xn])

• .

Starting with [x0] = [1, 2] and applying interval Newton’s method, we get [x1] = [1.375, 1.4375] [x2] = [1.4140625, 1.41441761363636] [x3] = [1.41421355929452, 1.41421356594718] [x4] = [1.41421356237309, 1.4142135623731].

11

Automatic Differentiation

Denote the i-th Taylor coefficient of u(t) evaluated at some point tj by (uj)i = u(i)(tj) i! . If we know the Taylor coefficients of u(t) and v(t) at tj, then (uj ± vj)i = (uj)i ± (vj)i (ujvj)i =

i

• r=0

(uj)r (vj)i−r uj vj

• i

= 1 vj

• (uj)i −

i

• r=1

(vj)r uj vj

• i−r
• 12

Application of Automatic Differentiation to ODEs

If y′(t) = f (y) and y(tj) = yj, then (yj)0 = f [0](yj) = yj, (yj)i = f [i](yj) = 1 i (f (yj))i−1, for i ≥ 1. We require at most O(k2) work to compute (yj)1 , (yj)2 , . . . , (yj)k. Note also we can replace the scalars above by intervals. For example, if (uj)i ∈

• ui

j

• and (vj)i ∈
• vi

j

• , then

(uj ± vj)i ∈

• (uj ± vj)i
• =
• ui

j

• ±
• vi

j

• 13

Example

Consider the IVP for the ODE y′

1

= f1(y1, y2) = y1y2, y1(0) = 1, y′

2

= f2(y1, y2) = y1 + y2, y2(0) = 2. Compute the Taylor series coefficients at t = 0. (y1)0 = 1 (y2)0 = 2 (y1)1 = f1(y1, y2) = y1y2 = 2 (y2)1 = f2(y1, y2) = y1 + y2 = 3

14

Example continued

(y1)2 =

1 2 {f1(y1, y2)}1 = 1 2(y1y2)1

=

1 2 {(y1)0(y2)1 + (y1)1(y2)0}

= 7/2 (y2)2 =

1 2 {f2(y1, y2)}1 = 1 2(y1 + y2)1

=

1 2 {(y1)1 + (y2)1}

= 5/2 (y1)3 =

1 3 {f1(y1, y2)}2 = 1 3(y1y2)2

=

1 3 {(y1)0(y2)2 + (y1)1(y2)1 + (y1)2(y2)0}

= 31/6 (y2)3 =

1 3 {f2(y1, y2)}2 = 1 3(y1 + y2)2

=

1 3 {(y1)2 + (y2)2}

= 2

15

Validated Numerical Methods for IVPs for ODEs

For f : Rn → Rn, we consider the IVP for ODEs y′(t) = f (y(t)), y(t0) = y0 We denote the solution by y(t; t0, y0). For an interval [y0], y(t; t0, [y0]) = {y(t; t0, y0) | y0 ∈ [y0]}.

t0 t2 y [y2] [y3] [y1] t1 t3 t [y0] True solutions

Our goal is to find interval vectors [yj] that are guaranteed to con- tain y(tj; t0, [y0]) at the points t0 < t1 < t2 < · · · < tN = T.

16

Advantages and Disadvantages of Interval Methods for IVPs for ODEs

Advantages of interval methods over standard numerical methods:

1 ensure a unique solution exists 2 provide guaranteed bounds on the solution 1

can be used to prove theorems

2

ensure critical calculations are accurate

3

helpful in testing standard numerical method

3 can be efficient for problems with ranges of parameters

Disadvantages of interval methods:

1 computation is time consuming 2 harder to implement than standard methods 3 error bounds may be too large 17

Taylor Series Methods for IVPs for ODEs

Suppose that we have computed [yj] at tj, such that y(tj; t0, [y0]) ⊆ [yj]. We advance the solution by using two algorithms. Algorithm I validates existence and uniqueness of the solution on [tj, tj+1] and computes an a priori enclosure [˜ yj] such that y(t; tj, [yj]) ⊆ [˜ yj] for t ∈ [tj, tj+1] . Algorithm II computes a tighter enclosure [yj+1] of y(tj+1; t0, [y0]).

[yj+1] tj tj+1 [yj] [˜ yj] y(t; tj, [yj])

18

Algorithm I: Validating Existence and Uniqueness

First Order Enclosure (FOE) Method (Eijgenraam, 1981) If [yj] ⊆ [˜ yj] and the inclusion [yj] + [0, hj]f ([˜ yj]) ⊆ [˜ yj] holds, then there exists a unique solution y(t; tj, yj) ∈ [˜ y0

j ] ≡ [yj] + [0, hj]f ([˜

yj]) ⊆ [˜ yj] for t ∈ [tj, tj+1] and all yj ∈ [yj]. The proof is based on the Banach fixed-point theorem and the Picard-Lindel¨

• f operator

(Tu)(t) = yj + t

tj

f (u(s))ds.

19

Algorithm II: Computing a Tight Enclosure

Using [˜ yj], compute a tighter enclosure [yj+1] for y(tj+1; tj, [yj]) ⊆ [yj+1]. Basic approach: Taylor series + a remainder term. yj+1 = yj +

k−1

• i=1

f [i](yj)hi

j + f [k](y; tj, tj+1)hk j ,

for some y ∈ [˜ yj]. Hence [yj+1] = [yj] +

k−1

• i=1

f [i]([yj])hi

j + f [k]([˜

yj])hk

j ,

True Solutions

Problem: the width of [yj] always increases with j even if the true solutions contract. Can we do better?

20

Algorithm II: A Better Tight Enclosure

Use the mean-value form. For any yj, ˆ yj ∈ [yj], yj+1 = yj +

k−1

• i=1

f [i](yj)hi

j + f [k](y; tj, tj+1)hk j

∈ ˆ yj +

k−1

• i=1

f [i](ˆ yj)hi

j + f [k]([˜

yj])hk

j

+

• I +

k−1

• i=1

J(f [i]; [yj])hi

j

• ([yj] − ˆ

yj)

y(t; tj, ˆ yj) Enclosure of y(tj+1; tj,

• yj
• )

ˆ yj Enclosure of y(tj+1; tj, ˆ yj) y(t; tj,

• yj
• )

tj tj+1

• yj
• t

y

21

The Wrapping Effect

Moore’s example: y′

1 = y2,

y′

2 = −y1;

y0 ∈ [y0] . The solution is y(t) = A(t)y0, where y0 ∈ [y0] and A(t) =

• cos t

sin t − sin t cos t

• .
• 10
• 5

5 10

• 2

2 4 6 8 10 12

[y0] ∈ IR2 can be viewed as a paral- lelepiped. At each step the parallelepiped is rotated and has to be wrapped by another one. At t = 2π the blow up is by a factor of e2π ≈ 535, as the stepsize tends to zero.

22

A Method for Reducing the Wrapping Effect

The Direct Method, that we have considered so far, is of the form [yj+1] = uj+1 + [zj+1] + [Sj]([yj] − ˆ yj) On each step, represent the enclosure in the form yj − ˆ yj ∈ {Ajrj | rj ∈ [rj]} . Then proceed as follows

1 Set [r0] = [y0] − ˆ

y0, A0 = I

2 Compute [yj+1] = uj+1 + [zj+1] + ([Sj]Aj)[rj] 3 Let ˆ

yj+1 = uj+1 + m([zj+1]) (m(·) – midpoint)

4 Choose Aj+1 5 Propagate [rj+1] =

• A−1

j+1([Sj]Aj)

• [rj] + A−1

j+1([zj+1] − m([zj+1])).

The Parallelepiped Method: Aj+1 ∈ [Sj] Aj Lohner’s QR Method: Aj+1 = Qj+1, where Aj+1 ∈ [Sj] Aj and Qj+1Rj+1 = Aj+1Pj+1

23

An Illustration of Lohner’s QR method

Qy y y {r | r ∈ [r]} x Qx x y x y

• Qx
• Qy

x S [r] {Sr | r ∈ [r]} {Sr | r ∈ [r]} {Sr | r ∈ [r]}

QR = SP

24

An Example

Consider the problem y′

1 = y2,

y′

2 = −y1

for t ∈ [0, 10] with interval initial conditions y1(0) ∈ [0.9, 1.1], y2(0) ∈ [−0.1, 0.1] We integrate with a constant stepsize 0.4 and Lohner’s QR Method.

• 1.5
• 1
• 0.5

0.5 1 1.5 2 4 6 8 10 Tight Enclosure of Y(1) A priori Enclosure of Y(1)

• 1.5
• 1
• 0.5

0.5 1 1.5 2 4 6 8 10 Tight Enclosure of Y(2) A priori Enclosure of Y(2)

• 1.5
• 1
• 0.5

0.5 1 1.5

• 1.5
• 1
• 0.5

0.5 1 1.5 Phase Portrait

25

A Simple Test Problem for the Wrapping Effect

We analyze Interval Taylor Series (ITS) Methods with the Direct, Parallelepiped and Lohner’s QR Methods for controlling the wrapping effect when applied to the simple test problem y′ = By, y(0) = y0 and compare the enclosure each interval method produces to the error bound for the Point Taylor Series (PTS) Method. The Taylor polynomial T =

k

• i=0

(hB)i i! is central to all methods. For example, for the Point Taylor Series Method, yj+1 = Tyj.

26

A Summary of Our Analysis of the Wrapping Effect

For the PTS Method, the error growth in the computation depends on T k ≈ cρ(T)k, where ρ(T) = max{|λ| : λ an eigenvalue of T} is the spectral radius of T. Therefore, the error growth depends on ρ(T). For the ITS Method with the Direct Method, the error bound growth depends on ρ(|T|). Note, ρ(|T|) ≥ ρ(T) and often ρ(|T|) ≫ ρ(T).

27

A Summary of Our Analysis of the Wrapping Effect

The ITS Method with the Parallelepiped Method works will if all eigenvalues of T are of the same magnitude. If this is not the case, the Parallelepiped Method breaks down because the matrices Aj = T j become computationally singular. If the eigenvalues {λi : i = 1, . . . , n} of T satisfy |λ1| > |λ2| > |λ3| > · · · > |λn−1| > |λn| > 0 (1) then the growth of the error bounds for the ITS Method with Lohner’s QR Method is very similar to the growth of the error in the PTS Method. Even if the eigenvalues of T don’t satisfy (1), the ITS Method with Lohner’s QR Method usually performs very well, but we can’t prove it.

28

The End Thank You

29