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 Preliminaries 1 Standard Versus Validated Numerical Methods Interval Arithmetic Automatic Differentiation Validated Numerical Methods for IVPs for ODEs 2 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 + ¯ [ a ] + [ b ] = ¯ + b a ¯ , ¯ b � � ¯ − ¯ [ a ] − [ b ] = b , ¯ a − b a � � ¯ � � �� ¯ a ¯ ¯ a ¯ [ a ] [ b ] = min a ¯ b ¯ , a b , ¯ ab ¯ , ¯ b , max a ¯ b ¯ , a b , ¯ ab ¯ , ¯ b � ¯ � ¯ 1 / ¯ [ a ] / [ b ] = [ a ¯ , ¯ a ] b , 1 / b , 0 / ∈ [ b ] . ¯ 4

Some Properties of Interval Arithmetic Interval arithmetic is Inclusion Monotone: If [ a 1 ] ⊆ [ a 2 ] and [ b 1 ] ⊆ [ b 2 ], then [ a 1 ] ∗ [ b 1 ] ⊆ [ a 2 ] ∗ [ b 2 ] , ∗ ∈ { + , − , × , / } . 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 ] ⇐ ⇒ [ a ij ] ⊆ [ b ij ] ∀ 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 : R n → R is continuously differentiable on D ⊆ R n , then for [ a ] ⊆ D and any y and ˆ y ∈ [ a ] f ( y ) = f (ˆ y ) + f y ( η )( y − ˆ y ) for some η ∈ [ a ]. Therefore R ( f ; [ a ]) ⊆ f M ([ a ] , ˆ y ) ≡ f (ˆ y ) + f y ([ a ])([ a ] − ˆ y ) . Often f M ([ a ] , ˆ y ) is a better enclosure of R ( f ; [ a ]) than f ([ a ]). 7

Example � 1 � Consider f ( x ) = x − x 2 for x ∈ [ a ] = 4 , 3 . 4 The range of f ( x ) is R ( f ; [ a ]) = 1 16 [3 , 4]. � 1 � = 1 16 [3 , 4] 16 . w For f ( x ) = x − x 2 , 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] � 1 � w 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] � 1 � = 1 w 16 [1 , 9] 2 . Mean-value evaluation gives: f M ([ 1 4 , 3 4 ] , 1 2 ) = f ( 1 2 ) + f ′ ([ 1 4 , 3 4 ]) × ([ 1 4 , 3 4 ] − 1 2 ) = 1 16 [2 , 6] � 1 � = 1 16 [2 , 6] 4 . w 8

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

A Better Approach: Interval Newton’s Method Suppose f (ˆ x ) = 0 and ˆ x ∈ [ x n ]. Let x n be some point in [ x n ]. (Often a good choice is x n = m ([ x n ]).) Then, for some η ∈ [ x n ], 0 = f (ˆ x ) = f ( x n ) + f ′ ( η )(ˆ x − x n ) Therefore ˆ = x n − f ( x n ) / f ′ ( η ) x x n − f ( x n ) / f ′ ([ x n ]) x ˆ ∈ � � x n − f ( x n ) / f ′ ([ x n ]) x ˆ ∈ [ x n ] ∩ Interval Newton’s Method: � � x n − f ( x n ) / f ′ ([ x n ]) [ x n +1 ] = [ x n ] ∩ Note ˆ x ∈ [ x n +1 ] ⊆ [ x n ]. Under suitable conditions, w ([ x n ]) decreases quadratically with n . 10

Example Find the positive root of f ( x ) = x 2 − 2. √ The solution is 2 = 1 . 4142135 ... . For this example, Interval Newton’s Method is � � x n − ( x 2 [ x n +1 ] = [ x n ] ∩ n − 2) / (2[ x n ]) . Starting with [ x 0 ] = [1 , 2] and applying interval Newton’s method, we get [ x 1 ] = [1 . 375 , 1 . 4375] [ x 2 ] = [1 . 4140625 , 1 . 41441761363636] [ x 3 ] = [1 . 41421355929452 , 1 . 41421356594718] [ x 4 ] = [1 . 41421356237309 , 1 . 4142135623731] . 11

Automatic Differentiation Denote the i -th Taylor coefficient of u ( t ) evaluated at some point t j by ( u j ) i = u ( i ) ( t j ) . i ! If we know the Taylor coefficients of u ( t ) and v ( t ) at t j , then ( u j ± v j ) i = ( u j ) i ± ( v j ) i i � ( u j v j ) i = ( u j ) r ( v j ) i − r r =0 � � � u j � � u j � � i 1 = ( u j ) i − ( v j ) r v j v j v j i i − r r =1 12

Application of Automatic Differentiation to ODEs If y ′ ( t ) = f ( y ) and y ( t j ) = y j , then f [0] ( y j ) = y j , ( y j ) 0 = f [ i ] ( y j ) = 1 ( y j ) i = i ( f ( y j )) i − 1 , for i ≥ 1 . We require at most O ( k 2 ) work to compute ( y j ) 1 , ( y j ) 2 , . . . , ( y j ) k . Note also we can replace the scalars above by intervals. � � � � u i v i For example, if ( u j ) i ∈ and ( v j ) i ∈ , then j j � � � � � � u i v i ( u j ± v j ) i ∈ ( u j ± v j ) i = ± j j 13

Example Consider the IVP for the ODE y ′ = f 1 ( y 1 , y 2 ) = y 1 y 2 , y 1 (0) = 1 , 1 y ′ = f 2 ( y 1 , y 2 ) = y 1 + y 2 , y 2 (0) = 2 . 2 Compute the Taylor series coefficients at t = 0. ( y 1 ) 0 = 1 ( y 2 ) 0 = 2 ( y 1 ) 1 = f 1 ( y 1 , y 2 ) = y 1 y 2 = 2 ( y 2 ) 1 = f 2 ( y 1 , y 2 ) = y 1 + y 2 = 3 14

Example continued 1 2 { f 1 ( y 1 , y 2 ) } 1 = 1 ( y 1 ) 2 = 2 ( y 1 y 2 ) 1 1 = 2 { ( y 1 ) 0 ( y 2 ) 1 + ( y 1 ) 1 ( y 2 ) 0 } = 7 / 2 1 2 { f 2 ( y 1 , y 2 ) } 1 = 1 ( y 2 ) 2 = 2 ( y 1 + y 2 ) 1 1 = 2 { ( y 1 ) 1 + ( y 2 ) 1 } = 5 / 2 1 3 { f 1 ( y 1 , y 2 ) } 2 = 1 ( y 1 ) 3 = 3 ( y 1 y 2 ) 2 1 = 3 { ( y 1 ) 0 ( y 2 ) 2 + ( y 1 ) 1 ( y 2 ) 1 + ( y 1 ) 2 ( y 2 ) 0 } = 31 / 6 3 { f 2 ( y 1 , y 2 ) } 2 = 1 1 ( y 2 ) 3 = 3 ( y 1 + y 2 ) 2 1 = 3 { ( y 1 ) 2 + ( y 2 ) 2 } = 2 15

Validated Numerical Methods for IVPs for ODEs For f : R n → R n , we consider the IVP for ODEs y ′ ( t ) = f ( y ( t )) , y ( t 0 ) = y 0 We denote the solution by y ( t ; t 0 , y 0 ). For an interval [ y 0 ], y ( t ; t 0 , [ y 0 ]) = { y ( t ; t 0 , y 0 ) | y 0 ∈ [ y 0 ] } . y Our goal is to find interval vectors [ y 0 ] [ y j ] that are guaranteed to con- tain y ( t j ; t 0 , [ y 0 ]) at the points [ y 1 ] [ y 2 ] t 0 < t 1 < t 2 < · · · < t N = T . [ y 3 ] True solutions t 0 t 1 t 2 t 3 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 can be used to prove theorems 1 ensure critical calculations are accurate 2 helpful in testing standard numerical method 3 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 [ y j ] at t j , such that y ( t j ; t 0 , [ y 0 ]) ⊆ [ y j ]. We advance the solution by using two algorithms. Algorithm I validates existence and uniqueness of the solution on [ t j , t j +1 ] and computes an a priori enclosure [˜ y j ] such that y ( t ; t j , [ y j ]) ⊆ [˜ y j ] for t ∈ [ t j , t j +1 ] . Algorithm II computes a tighter enclosure [ y j +1 ] of y ( t j +1 ; t 0 , [ y 0 ]). [˜ y j ] y ( t ; t j , [ y j ]) [ y j +1 ] [ y j ] t j t j +1 18