Round-off Error Analysis of Explicit One-Step Numerical Integration - - PowerPoint PPT Presentation

round off error analysis of explicit one step numerical
SMART_READER_LITE
LIVE PREVIEW

Round-off Error Analysis of Explicit One-Step Numerical Integration - - PowerPoint PPT Presentation

Round-off Error Analysis of Explicit One-Step Numerical Integration Methods 24th IEEE Symposium on Computer Arithmetic Sylvie Boldo 1 Florian Faissole 1 Alexandre Chapoutot 2 1 Inria - LRI, Univ. Paris-Sud et CNRS - Univ. Paris-Saclay 2 U2IS,


slide-1
SLIDE 1

Round-off Error Analysis of Explicit One-Step Numerical Integration Methods

24th IEEE Symposium on Computer Arithmetic⋆ Sylvie Boldo1 Florian Faissole1 Alexandre Chapoutot2

1Inria - LRI, Univ. Paris-Sud et CNRS - Univ. Paris-Saclay 2U2IS, ´

ENSTA ParisTech ⋆ we thank the IEEE for registration bursary

slide-2
SLIDE 2

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Table of contents

1

Motivations and numerical methods

2

Roundoff errors of RK methods Local roundoff errors Global roundoff errors of classical methods

3

Conclusion and perspectives

2 / 34

slide-3
SLIDE 3

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Ordinary differential equations (ODEs)

y′(t) = f(y,t). Exact resolution is hard ⇒ numerical methods.

3 / 34

slide-4
SLIDE 4

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

4 / 34

slide-5
SLIDE 5

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

5 / 34

slide-6
SLIDE 6

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

6 / 34

slide-7
SLIDE 7

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

7 / 34

slide-8
SLIDE 8

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

8 / 34

slide-9
SLIDE 9

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

9 / 34

slide-10
SLIDE 10

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

10 / 34

slide-11
SLIDE 11

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Numerical integration

11 / 34

slide-12
SLIDE 12

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Euler method

k1 = hf(tn,yn) yn+1 = yn + k1 + O(h2)

12 / 34

slide-13
SLIDE 13

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

RK2 method

k1 = h × f(tn,yn) k2 = h × f(tn + h

2 ,yn + k1 2 )

yn+1 = yn + k2 + O(h3)

13 / 34

slide-14
SLIDE 14

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Working assumptions

  • FP arithmetic;
  • neither underflow nor overflow,
  • radix 2 double precision,

u = 2−53

14 / 34

slide-15
SLIDE 15

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Working assumptions

  • FP arithmetic;
  • neither underflow nor overflow,
  • radix 2 double precision,

u = 2−53

  • ODEs;
  • first-order,
  • linear,

y′ = λy

  • y ∶ R → R,

14 / 34

slide-16
SLIDE 16

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Working assumptions

  • FP arithmetic;
  • neither underflow nor overflow,
  • radix 2 double precision,

u = 2−53

  • ODEs;
  • first-order,
  • linear,

y′ = λy

  • y ∶ R → R,
  • Methods.
  • explicit,
  • one step,
  • constant step.

14 / 34

slide-17
SLIDE 17

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

RK methods on linear problems: linear stability

{ y0 ∈ R yn+1 = R(h,λ)yn (R polynomial in hλ)

15 / 34

slide-18
SLIDE 18

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

RK methods on linear problems: linear stability

{ y0 ∈ R yn+1 = R(h,λ)yn (R polynomial in hλ) Stable ⇔ ∣R(h,λ)∣ < 1:

15 / 34

slide-19
SLIDE 19

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

RK methods on linear problems: FP implementation

{ y0 ∈ R yn+1 = R(h,λ)yn { ̃ y0 ≃ y0 ̃ yn+1 = ̃ R(̃ h,̃ λ, ̃ yn)

16 / 34

slide-20
SLIDE 20

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

RK methods on linear problems: FP implementation

{ y0 ∈ R yn+1 = R(h,λ)yn { ̃ y0 ≃ y0 ̃ yn+1 = ̃ R(̃ h,̃ λ, ̃ yn) Euler:

  • R(h,λ) = 1 + hλ;
  • ̃

R(̃ h,̃ λ, ̃ yn) = ̃ yn ⊕ ̃ h ⊗ ̃ λ ⊗ ̃ yn.

16 / 34

slide-21
SLIDE 21

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

RK methods on linear problems: FP implementation

{ y0 ∈ R yn+1 = R(h,λ)yn { ̃ y0 ≃ y0 ̃ yn+1 = ̃ R(̃ h,̃ λ, ̃ yn) Euler:

  • R(h,λ) = 1 + hλ;
  • ̃

R(̃ h,̃ λ, ̃ yn) = ̃ yn ⊕ ̃ h ⊗ ̃ λ ⊗ ̃ yn. RK4:

  • R(h,λ) = 1 + hλ + 1

2(hλ)2 + 1 6(hλ)3 + 1 24(hλ)4;

  • ̃

R(̃ h,̃ λ, ̃ yn) =

̃ yn ⊕̃ h ⊘ 6 ⊗ ̃ λ ⊗ ̃ yn ⊕̃ h ⊘ 3 ⊗ ̃ λ ⊗ ̃ yn ⊕̃ h ⊗̃ h ⊘ 6 ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ yn ⊕̃ h ⊘ 3 ⊗ ̃ λ ⊗ ̃ yn⊗ ̃ h ⊗ ̃ h ⊘ 6 ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ yn ⊕ ̃ h ⊗ ̃ h ⊗ ̃ h ⊘ 12 ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ yn ⊕ ̃ h ⊘ 6 ⊗ ̃ λ ⊗ ̃ yn ⊕ ̃ h ⊗ ̃ h ⊘ 6 ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ yn ⊕ ̃ h ⊗ ̃ h ⊗ ̃ h ⊘ 12 ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ yn ⊕ ̃ h ⊗ ̃ h ⊗ ̃ h ⊗ ̃ h ⊘ 24 ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ λ ⊗ ̃ yn.

(> 60 flops!)

16 / 34

slide-22
SLIDE 22

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

State-of-the-art

Roundoff errors in numerical methods (N = nb of iterations):

  • probabilistic result: error in

√ N [Henrici,1963];

  • in practice (implicit RK): error in N [Hairer&al, 2008];
  • interval analysis [Bouissou-Martel, 2006];
  • numerical integration (fine-grained): Newton-Cotes,

Gauss-Legendre, ... [Fousse, 2006].

17 / 34

slide-23
SLIDE 23

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

State-of-the-art

Roundoff errors in numerical methods (N = nb of iterations):

  • probabilistic result: error in

√ N [Henrici,1963];

  • in practice (implicit RK): error in N [Hairer&al, 2008];
  • interval analysis [Bouissou-Martel, 2006];
  • numerical integration (fine-grained): Newton-Cotes,

Gauss-Legendre, ... [Fousse, 2006].

Our approach:

  • fined-grained analysis;
  • use of mathematical properties of the methods (stability).

17 / 34

slide-24
SLIDE 24

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Table of contents

1

Motivations and numerical methods

2

Roundoff errors of RK methods Local roundoff errors Global roundoff errors of classical methods

3

Conclusion and perspectives

18 / 34

slide-25
SLIDE 25

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Method error vs roundoff error

19 / 34

slide-26
SLIDE 26

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Method error vs roundoff error

20 / 34

slide-27
SLIDE 27

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Method error vs roundoff error

21 / 34

slide-28
SLIDE 28

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Local roundoff error vs global roundoff error

Local error: ε0 = ∣̃ y0 − y0∣ ∀n ∈ N∗, εn = ∣ ̃ R(̃ h,̃ λ, ̃ yn−1) − R(h,λ) ̃ yn−1∣. Global error: ∀n ∈ N, En = ̃ yn − yn.

22 / 34

slide-29
SLIDE 29

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

From local to global roundoff error

Local error: ε0 = ∣̃ y0 − y0∣ ∀n ∈ N∗, εn = ∣ ̃ R(̃ h,̃ λ, ̃ yn−1) − R(h,λ) ̃ yn−1∣. Global error: ∀n ∈ N, En = ̃ yn − yn. Theorem 1: Global absolute error of RK methods Let C ∈ R∗

+. Suppose ∀n ∈ N∗,εn ⩽ C∣ ̃

yn−1∣. Then, ∀n ∈ N, ∣En∣ ⩽ (C + ∣R(h,λ)∣)n (ε0 + n C∣y0∣ C + ∣R(h,λ)∣).

23 / 34

slide-30
SLIDE 30

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Relative roundoff errors

Relative error: ∣ ̃ yn − yn yn ∣ ⩽ (C + ∣R(h,λ)∣ ∣R(h,λ)∣ )

n

(ε0 + n C∣y0∣ C + ∣R(h,λ)∣).

24 / 34

slide-31
SLIDE 31

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Relative roundoff errors

Relative error: ∣ ̃ yn − yn yn ∣ ⩽ (C + ∣R(h,λ)∣ ∣R(h,λ)∣ )

n

(ε0 + n C∣y0∣ C + ∣R(h,λ)∣). If C ≪ ∣R(h,λ)∣, then: ∣ ̃ yn − yn yn ∣ ≲ ε0 + n C∣y0∣ ∣R(h,λ)∣. In practice (Euler, RK2, RK4): C ≤ 200u and 200u ≪ ∣R(h,λ)∣.

24 / 34

slide-32
SLIDE 32

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Table of contents

1

Motivations and numerical methods

2

Roundoff errors of RK methods Local roundoff errors Global roundoff errors of classical methods

3

Conclusion and perspectives

25 / 34

slide-33
SLIDE 33

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Technical lemma for local roundoff errors

Local error of stable Euler’s method (−2 ≤ hλ < 0): εn+1 = ∣ ̃ yn ⊕ (̃ h ⊗ ̃ λ ⊗ ̃ yn) − (1 + hλ)yn∣.

26 / 34

slide-34
SLIDE 34

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Technical lemma for local roundoff errors

Local error of stable Euler’s method (−2 ≤ hλ < 0): εn+1 = ∣ ̃ yn

  • X1

⊕(̃ h ⊗ ̃ λ

  • X2

⊗ ̃ yn

  • y

) − ( 1

  • α1

+ hλ

  • α2

) ̃ yn

  • y

∣. Lemma 2 Let y ∈ R. Let C1,C2 ∈ R+. Let α1,α2 ∈ R. Let X1,X2 ∈ F s.t. ∣X1 − α1y∣ ⩽ C1∣y∣ and ∣X2 − α2∣ ≤ C2. Then: ∣X1 ⊕ (X2 ⊗ y) − (α1 + α2)y∣ ⩽ ∣y∣(C1 + C2 + u(∣α1∣ + 2∣α2∣ + C1 + 2C2) +u2 (C2 + ∣α2∣)).

27 / 34

slide-35
SLIDE 35

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Technical lemma for local roundoff errors

Local error of stable Euler’s method (−2 ≤ hλ < 0): εn+1 = ∣ ̃ yn

  • X1

⊕(̃ h ⊗ ̃ λ

  • X2

⊗ ̃ yn

  • y

) − ( 1

  • α1

+ hλ

  • α2

) ̃ yn

  • y

∣. Lemma 2 Let y ∈ R. Let C1,C2 ∈ R+. Let α1,α2 ∈ R. Let X1,X2 ∈ F s.t. ∣X1 − α1y∣ ⩽ C1∣y∣ and ∣X2 − α2∣ ≤ C2. Then: ∣X1 ⊕ (X2 ⊗ y) − (α1 + α2)y∣ ⩽ ∣y∣(C1 + C2 + u(∣α1∣ + 2∣α2∣ + C1 + 2C2) +u2 (C2 + ∣α2∣)). ∣ ̃ yn

  • X1

− 1

  • α1

̃ yn

  • y

∣ ≤ 0

  • C1

∣ ̃ yn

  • y

∣, ∣ ̃ h ⊗ ̃ λ

  • X2

− hλ

  • α2

∣ ≤ 6u

  • C2

(Gappa [Melquiond]).

27 / 34

slide-36
SLIDE 36

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Local roundoff errors of higher-order methods

Stable RK4 method: C ̃ yn ○[̃ h 1

λ] 2u (Gappa) ○[ ̃ yn + ̃ h 1

λ ̃ yn] = ̃ yn ⊕ ̃ h ⊘ 6 ⊗ ̃ λ ⊗ ̃ yn 4u (Lemma 2) ○[̃ h 1

λ] 4u (Gappa) ○[̃ h̃ h 1

λ̃ λ] 12u (Gappa) ○[̃ h̃ h̃ h 1

12̃

λ̃ λ̃ λ] 28u (Gappa) ○[̃ h̃ h̃ h̃ h 1

24̃

λ̃ λ̃ λ̃ λ] 53u (Gappa) ... ... (7 × Lemma 2) ○[ ̃ yn + ⋅⋅⋅ + ̃ h̃ h̃ h̃ h 1

24̃

λ̃ λ̃ λ̃ λ ̃ yn] 194u (Lemma 2)

28 / 34

slide-37
SLIDE 37

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Local roundoff errors of classical methods

Lemma 3: Local error of Euler method Suppose −2 ⩽ hλ ⩽ −2−100 and 2−60 ⩽ h ⩽ 1. Then: ∀n ∈ N, εEuler

n+1

⩽ 11.01u ∣ ̃ yn∣. Lemma 4: Local error of RK2 method Suppose −2 ⩽ hλ ⩽ −2−100 and 2−60 ⩽ h ⩽ 1. Then: ∀n ∈ N, εRK2

n+1 ⩽ 28.01u∣ ̃

yn∣. Lemma 5: Local error of RK4 method Suppose −3 ⩽ hλ ⩽ −2−100 and 2−60 ⩽ h ⩽ 1. Then: ∀n ∈ N, εRK4

n+1 ⩽ 194u∣ ̃

yn∣.

29 / 34

slide-38
SLIDE 38

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Table of contents

1

Motivations and numerical methods

2

Roundoff errors of RK methods Local roundoff errors Global roundoff errors of classical methods

3

Conclusion and perspectives

30 / 34

slide-39
SLIDE 39

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Global roundoff errors of classical methods

Bounds on global errors:

  • Euler: C = 11.01u

∣En∣ ⩽ (11.01u + ∣R(h,λ)∣)n (ε0 + n 11.01u∣y0∣ 11.01u + ∣R(h,λ)∣);

  • RK2: C = 28.01u

∣En∣ ⩽ (28.01u + ∣R(h,λ)∣)n (ε0 + n 28.01u∣y0∣ 28.01u + ∣R(h,λ)∣);

  • RK4: C = 194u

∣En∣ ⩽ (194u + ∣R(h,λ)∣)n (ε0 + n 194u∣y0∣ 194u + ∣R(h,λ)∣). ⇒ no compensation but reasonable bounds .

31 / 34

slide-40
SLIDE 40

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Table of contents

1

Motivations and numerical methods

2

Roundoff errors of RK methods Local roundoff errors Global roundoff errors of classical methods

3

Conclusion and perspectives

32 / 34

slide-41
SLIDE 41

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

General methodology

yn+1 =

M

i=0

αiyn ̃ yn+1 = ⊕M

i=0 ̃

αi⊗ ̃ yn (αi = hλ, hλ 3 , hλ 6 , h4λ4 24 ...)

33 / 34

slide-42
SLIDE 42

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

General methodology

yn+1 =

M

i=0

αiyn ̃ yn+1 = ⊕M

i=0 ̃

αi⊗ ̃ yn (αi = hλ, hλ 3 , hλ 6 , h4λ4 24 ...) Methodology to bound roundoff errors:

33 / 34

slide-43
SLIDE 43

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

General methodology

yn+1 =

M

i=0

αiyn ̃ yn+1 = ⊕M

i=0 ̃

αi⊗ ̃ yn (αi = hλ, hλ 3 , hλ 6 , h4λ4 24 ...) Methodology to bound roundoff errors: 1) bound ε0 = ∣̃ y0 − y0∣;

33 / 34

slide-44
SLIDE 44

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

General methodology

yn+1 =

M

i=0

αiyn ̃ yn+1 = ⊕M

i=0 ̃

αi⊗ ̃ yn (αi = hλ, hλ 3 , hλ 6 , h4λ4 24 ...) Methodology to bound roundoff errors: 1) bound ε0 = ∣̃ y0 − y0∣; 2) bound the error on each term αi (Gappa + stability);

33 / 34

slide-45
SLIDE 45

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

General methodology

yn+1 =

M

i=0

αiyn ̃ yn+1 = ⊕M

i=0 ̃

αi⊗ ̃ yn (αi = hλ, hλ 3 , hλ 6 , h4λ4 24 ...) Methodology to bound roundoff errors: 1) bound ε0 = ∣̃ y0 − y0∣; 2) bound the error on each term αi (Gappa + stability); 3) bound local errors by M applications of Lemma 2;

33 / 34

slide-46
SLIDE 46

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

General methodology

yn+1 =

M

i=0

αiyn ̃ yn+1 = ⊕M

i=0 ̃

αi⊗ ̃ yn (αi = hλ, hλ 3 , hλ 6 , h4λ4 24 ...) Methodology to bound roundoff errors: 1) bound ε0 = ∣̃ y0 − y0∣; 2) bound the error on each term αi (Gappa + stability); 3) bound local errors by M applications of Lemma 2; 4) bound the global error by instantiating Theorem 1.

33 / 34

slide-47
SLIDE 47

Motivations and numerical methods Roundoff errors of RK methods Conclusion and perspectives

Conclusion and perspectives

Conclusion:

  • fine and mechanical analysis of roundoff errors;
  • results on useful and classical methods (Euler, RK2, RK4);
  • linear growth of roundoff errors;
  • no compensation (for explicit one-step methods).

Perspectives:

  • overflows and underflows;
  • formalization in Coq (proof assistant):
  • based on the Flocq library [Boldo-Melquiond],
  • based on the gappa and interval tactics [Melquiond],
  • more general ODEs: complex, matricial, non-linear;
  • more general methods: multi-step, variable step, implicit.

34 / 34