Solving Ordinary Differential Equations Sanzheng Qiao Department of - PowerPoint PPT Presentation

IVP Software Summary Solving Ordinary Differential Equations Sanzheng Qiao Department of Computing and Software McMaster University March, 2014

IVP Software Summary Outline Initial Value Problem 1 Euler’s Method Runge-Kutta Methods Multistep Methods Implicit Methods Hybrid Method Software Packages 2

IVP Software Summary Outline Initial Value Problem 1 Euler’s Method Runge-Kutta Methods Multistep Methods Implicit Methods Hybrid Method Software Packages 2

IVP Software Summary Problem setting Initial Value Problem (first order) find y ( t ) such that y ′ = f ( y , t ) , given initial value y ( t 0 ) . Usually we assume t 0 = 0

IVP Software Summary Problem setting Initial Value Problem (first order) find y ( t ) such that y ′ = f ( y , t ) , given initial value y ( t 0 ) . Usually we assume t 0 = 0 Generalization 1: system of first order ODEs: y is a vector and f a vector-valued function. Example � y ′ 1 = f 1 ( y 1 , y 2 , t ) y ′ 2 = f 2 ( y 1 , y 2 , t ) or in vector notations: y ′ = f ( y , t )

IVP Software Summary Problem setting (cont.) Generalization 2: high order equation u ′′ = g ( u , u ′ , t ) . Let y 1 = u u ′ y 2 = and transform the above into the following system of first order ODEs: � y ′ 1 = y 2 y ′ 2 = g ( y 1 , y 2 , t ) Initial conditions: u ( 0 ) and u ′ ( 0 )

IVP Software Summary Solution family A differential equation has a family of solutions, each corresponds to an initial value. 3.5 3 2.5 2 1.5 1 0.5 0 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 y ′ = − y , solution family y ( t ) = Ce − t , y ( 0 ) = C .

IVP Software Summary Euler’s method We consider the initial value problem: y ′ = f ( y , t ) , y ( t 0 ) = y 0 Numerical solution: find approximations y n ≈ y ( t n ) , for n = 1 , 2 , ... Note: y 0 = y ( t 0 ) (initial value)

IVP Software Summary Euler’s method We consider the initial value problem: y ′ = f ( y , t ) , y ( t 0 ) = y 0 Numerical solution: find approximations y n ≈ y ( t n ) , for n = 1 , 2 , ... Note: y 0 = y ( t 0 ) (initial value) A k -step method: Compute y n + 1 using y n − k + 1 , ..., y n − 1 , y n .

IVP Software Summary Euler’s method (cont.) A single-step method: Euler’s method. f ( y 0 , t 0 ) = y ′ ( t 0 ) ≈ y ( t 1 ) − y ( t 0 ) , h 0 where h 0 = t 1 − t 0 . The first step: y 1 = y 0 + h 0 f ( y 0 , t 0 )

IVP Software Summary Euler’s method (cont.) A single-step method: Euler’s method. f ( y 0 , t 0 ) = y ′ ( t 0 ) ≈ y ( t 1 ) − y ( t 0 ) , h 0 where h 0 = t 1 − t 0 . The first step: y 1 = y 0 + h 0 f ( y 0 , t 0 ) Euler’s method y n + 1 = y n + h n f ( y n , t n ) Produces: y 0 = y ( t 0 ) , y 1 ≈ y ( t 1 ) , y 2 ≈ y ( t 2 ) , ...

IVP Software Summary Example y ′ = − y , y ( 0 ) = 1 . 0. (Solution y = e − t )

IVP Software Summary Example y ′ = − y , y ( 0 ) = 1 . 0. (Solution y = e − t ) h = 0 . 4 Step 1: y 1 = y 0 − hy 0 = 1 . 0 − 0 . 4 × 1 . 0 = 0 . 6 ( ≈ y ( 0 . 4 ) = e − 0 . 4 ≈ 0 . 6703) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

IVP Software Summary Example u 1 ( t ) = 0 . 6 e − t + 0 . 4 ≈ 0 . 8951 e − t in the solution family. u ′ 1 = − u 1 , u 1 ( 0 ) ≈ 0 . 8951 ( u 1 ( 0 . 4 ) = 0 . 6) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

IVP Software Summary Example Step 2: y 2 = y 1 − hy 1 = 0 . 6 − 0 . 4 × 0 . 6 = 0 . 36

IVP Software Summary Example Step 2: y 2 = y 1 − hy 1 = 0 . 6 − 0 . 4 × 0 . 6 = 0 . 36 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

IVP Software Summary Example u 2 ( t ) = 0 . 36 e − t + 0 . 8 ≈ 0 . 8012 e − t in the solution family. u ′ 2 = − u 2 , u 2 ( 0 ) ≈ 0 . 8012 ( u 2 ( 0 . 8 ) = 0 . 36) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

IVP Software Summary Euler’s method In general u ′ n ( t ) = f ( u n ( t ) , t ) , in the solution family u n ( t n ) = y n , passing ( t n , y n ) u n ( t n + 1 ) ≈ u n ( t n ) + h n u ′ n ( t n ) = y n + h n f ( u n ( t n ) , t n ) = y n + h n f ( y n , t n ) = y n + 1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

IVP Software Summary Euler’s method Starting with t 0 and y 0 = y ( t 0 ) , as we proceed, we jump from one solution in the family to another. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

IVP Software Summary Errors Two sources of errors: discretization error and roundoff error. Discretization error : caused by the method used, independent of the computer used and the program implementing the method.

IVP Software Summary Errors Two sources of errors: discretization error and roundoff error. Discretization error : caused by the method used, independent of the computer used and the program implementing the method. Two types of discretization error: Global error: e n = y n − y ( t n ) Local error: the error in one step

IVP Software Summary Local error Consider t n as the starting point and the approximation y n at t n as the initial value , if u n ( t ) is the solution of u ′ n = f ( u n , t ) , u n ( t n ) = y n then the local error is d n = y n + 1 − u n ( t n + 1 )

IVP Software Summary Example Step 1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Local error d 0 = y 1 − y ( t 1 ) = 0 . 6 − e − 0 . 4 ≈ − 0 . 0703. Global error e 1 is the same as the local error d 0 .

IVP Software Summary Example Step 2 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Local error d 1 = y 2 − u 1 ( t 2 ) = 0 . 36 − u 1 ( 0 . 8 ) ≈ − 0 . 0422. Global error e 2 = y 2 − y ( t 2 ) = 0 . 36 − e − 0 . 8 ≈ − 0 . 0893.

IVP Software Summary Example 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 d 0 = y 1 − y ( t 1 ) = 0 . 6 − e − 0 . 4 ≈ − 0 . 0703 d 1 = y 2 − u 1 ( t 2 ) = 0 . 36 − u 1 ( 0 . 8 ) ≈ − 0 . 0422 e 2 = y 2 − y ( t 2 ) = 0 . 36 − e − 0 . 8 ≈ − 0 . 0893

IVP Software Summary Stability Relation between global error e n and local error d n If the differential equation is unstable, N − 1 � | e N | > | d n | n = 0 If the differential equation is stable, N − 1 � | e N | ≤ | d n | n = 0 In this case, � N − 1 n = 0 | d n | is an upper bound for the global error | e N | .

IVP Software Summary Example In the previous example: Local errors | d 0 | = 0 . 0703 and | d 1 | = 0 . 0422 Global error | e 2 | = 0 . 0893 | e 2 | < | d 0 | + | d 1 |

IVP Software Summary Example In the previous example: Local errors | d 0 | = 0 . 0703 and | d 1 | = 0 . 0422 Global error | e 2 | = 0 . 0893 | e 2 | < | d 0 | + | d 1 | More generally, y ′ = α y , solution family y = Ce α t . Stable when α < 0.

IVP Software Summary Accuracy A measurement for the accuracy of a method An order p method: | d n | ≤ Ch p + 1 ( or O ( h p + 1 )) n n C : independent of n and h n .

IVP Software Summary Example: Euler’s method y n + 1 = y n + h n f ( y n , t n ) Local solution u n ( t ) u ′ n ( t ) = f ( u n ( t ) , t ) , u n ( t n ) = y n Taylor expansion at t n : u n ( t ) = u n ( t n ) + ( t − t n ) u ′ ( t n ) + O (( t − t n ) 2 ) Since y n = u n ( t n ) and u ′ ( t n ) = f ( y n , t n ) , we get u n ( t n + 1 ) = y n + h n f ( y n , t n ) + O ( h 2 n ) Local error d n = y n + 1 − u n ( t n + 1 ) = O ( h 2 n ) Euler’s method is a first order method ( p = 1)

IVP Software Summary Accuracy (cont.) Consider the interval [ t 0 , t N ] and partition t 0 , t 1 , ..., t N . Roughly, the global error N − 1 � | d n | ≈ N · O ( h p + 1 ) ≈ ( t N − t 0 ) · O ( h p ) | e N | ≈ n = 0 at the final point t N is roughly O ( h p ) for a method of order p .

IVP Software Summary Accuracy (cont.) Consider the interval [ t 0 , t N ] and partition t 0 , t 1 , ..., t N . Roughly, the global error N − 1 � | d n | ≈ N · O ( h p + 1 ) ≈ ( t N − t 0 ) · O ( h p ) | e N | ≈ n = 0 at the final point t N is roughly O ( h p ) for a method of order p . For a p th order method, if the subintervals h n are cut in half, then the average local error is reduced by a factor of 2 p + 1 , the global error is reduced by a factor of 2 p . (But double the number of steps, i.e., more work.)

IVP Software Summary Roundoff Error Each step of the Euler’s method y n + 1 = y n + h n f ( y n , t n ) + ǫ | ǫ | = O ( u ) . Total rounding error: N ǫ = b ǫ/ h ( b = t N − t 0 , fixed step size h ) Ch + ǫ � � total error ≈ b h

IVP Software Summary Roundoff Error Each step of the Euler’s method y n + 1 = y n + h n f ( y n , t n ) + ǫ | ǫ | = O ( u ) . Total rounding error: N ǫ = b ǫ/ h ( b = t N − t 0 , fixed step size h ) Ch + ǫ � � total error ≈ b h Remarks If h is too small, the roundoff error is large If h is too large, the discretization error is large