Lecture 4: Numerical solution of ordinary differential equations - PowerPoint PPT Presentation

Numerical solution of ODEs • Second derivative can be found by differentiating the equation with respect to t : d 2 x dt 2 = d dt f ( t , x ) = ∂ f ∂ t ( t , x ) + ∂ f ∂ x ( t , x ) dx dt . • Second order Taylor method � ∂ f � x k +1 = x k + (∆ t ) f ( t k , x k ) + (∆ t ) 2 ∂ t ( t k , x k ) + ∂ f ∂ x ( t k , x k ) f ( t k , x k ) ( ∗ ) . 2 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Proposition: • Suppose that f ∈ C 2 . • Then ( ∗ ): of second order. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Proof: • f ∈ C 2 ⇒ x ∈ C 3 . • ⇒ truncation error T k given by T k (∆ t ) = (∆ t ) 2 d 3 x dt 3 ( τ ) , 6 for some τ ∈ [ t k , t k +1 ] and so, ( ∗ ): of second order. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Drawbacks of higher order Taylor methods: (i) Owing to their dependence upon the partial derivatives of f , f needs to be smooth; (ii) Efficient evaluation of the terms in the Taylor approximation and avoidance of round off errors. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Integral equation method • Avoid the complications inherent in a direct Taylor expansion. • x ( t ) coincides with the solution to the integral equation � t x ( t ) = x 0 + f ( s , x ( s )) ds , t ∈ [0 , T ] . 0 Starting at the discretization point t k instead of 0, and integrating until time t = t k +1 gives � t k +1 ( ∗∗ ) x ( t k +1 ) = x ( t k ) + f ( s , x ( s )) ds . t k • Implicitly computes the value of the solution at the subsequent discretization point. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Compare formula ( ∗∗ ) with the explicit Euler method x k +1 = x k + (∆ t ) f ( t k , x k ) . • ⇒ Approximation of the integral by � t k +1 f ( s , x ( s )) ds ≈ (∆ t ) f ( t k , x ( t k )) . t k • Left endpoint rule for numerical integration. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Left endpoint rule for numerical integration: • Left endpoint rule : not an especially accurate method of numerical integration. • Better methods include the Trapezoid rule: Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Numerical integration formulas for continuous functions. (i) Trapezoidal rule : � � � t k +1 g ( s ) ds ≈ ∆ t g ( t k +1 ) + g ( t k ) ; 2 t k (ii) Simpson’s rule : � � � t k +1 g ( s ) ds ≈ ∆ t g ( t k +1 ) + 4 g ( t k + t k +1 ) + g ( t k ) ; 6 2 t k (iii) Trapezoidal rule: exact for polynomials of order one ; Simpson’s rule: exact for polynomials of second order . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Use the more accurate Trapezoidal approximation � t k +1 � � f ( s , x ( s )) ds ≈ (∆ t ) f ( t k , x ( t k )) + f ( t k +1 , x ( t k +1 )) . 2 t k • Trapezoidal scheme: � � x k +1 = x k + (∆ t ) f ( t k , x k ) + f ( t k +1 , x k +1 ) . 2 • Trapezoidal scheme: implicit numerical method. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Proposition: • Suppose that f ∈ C 2 and (∆ t ) C f ( ∗ ∗ ∗ ) < 1; 2 C f : Lipschitz constant for f in x . • Trapezoidal scheme: convergent and of second order. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Proof: • Consistency: � � Φ( t , x , ∆ t ) := 1 f ( t , x ) + f ( t + ∆ t , x + (∆ t )Φ( t , x , ∆ t )) . 2 • ∆ t = 0. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Stability: • � � � ≤ C f | x − y | � Φ( t , x , ∆ t ) − Φ( t , y , ∆ t ) � � +∆ t � Φ( t , x , ∆ t ) − Φ( t , y , ∆ t ) � . 2 C f • ⇒ �� � � 1 − (∆ t ) C f � Φ( t , x , ∆ t ) − Φ( t , y , ∆ t ) � ≤ C f | x − y | . 2 • ⇒ Stability holds with C f C Φ = , 1 − (∆ t ) C f 2 provided that ∆ t satisfies ( ∗ ∗ ∗ ). Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Second order scheme: • By the mean-value theorem, x ( t k +1 ) − x ( t k ) T k (∆ t ) = ∆ t � � − 1 f ( t k , x ( t k )) + f ( t k +1 , x ( t k +1 )) 2 12(∆ t ) 2 d 3 x − 1 = dt 3 ( τ ) , for some τ ∈ [ t k , t k +1 ] ⇒ second order scheme, provided that f ∈ C 2 (and consequently x ∈ C 3 ). Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • An alternative scheme: replace x k +1 by x k + (∆ t ) f ( t k , x k ). • ⇒ Improved Euler scheme: � � x k +1 = x k + (∆ t ) f ( t k , x k ) + f ( t k +1 , x k + ( ∆t ) f ( t k , x k )) . 2 • Proposition: Improved Euler scheme: convergent and of second order. • Improved Euler scheme: performs comparably to the Trapezoidal scheme, and significantly better than the Euler scheme. • Alternative numerical approximations to the integral equation ⇒ a range of numerical solution schemes. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Midpoint rule : � t k +1 f ( s , x ( s )) ds ≈ (∆ t ) f ( t k + ∆ t 2 , x ( t k + ∆ t 2 )) . t k • Midpoint rule: same order of accuracy as the trapezoid rule. 2 ) by x k + ∆ t • Midpoint scheme: approximate x ( t k + ∆ t 2 f ( t k , x k ), � � t k + ∆ t 2 , x k + ∆ t x k +1 = x k + (∆ t ) f 2 f ( t k , x k ) . • Midpoint scheme: of second order. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Example of linear systems • Consider the linear system of ODEs  dx  dt = Ax ( t ) , t ∈ [0 , + ∞ [ ,  x (0) = x 0 ∈ R d . • A ∈ M d ( C ): independent of t . • DEFINITION: • A one-step numerical scheme for solving the linear system of ODEs: stable if there exists a positive constant C 0 s.t. | x k +1 | ≤ C 0 | x 0 | for all k ∈ N . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Consider the following schemes: (i) Explicit Euler’s scheme: x k +1 = x k + (∆ t ) Ax k ; (ii) Implicit Euler’s scheme: x k +1 = x k + (∆ t ) Ax k +1 ; (iii) Trapezoidal scheme: � � x k +1 = x k + (∆ t ) Ax k + Ax k +1 , 2 with k ∈ N , and x 0 = x 0 . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Proposition: Suppose that ℜ λ j < 0 for all j . The following results hold: (i) Explicit Euler scheme: stable for ∆ t small enough; (ii) Implicit Euler scheme: unconditionally stable; (iii) Trapezoidal scheme: unconditionally stable. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Proof: • Consider the explicit Euler scheme. By a change of basis, x k +1 = ( I + ∆ t ( D + N )) k � x 0 , � x k = Cx k . where � x 0 ∈ E j , then • If � min { k , d } � x k = k (1 + ∆ t λ j ) k − l (∆ t ) l N l � C l x 0 , � l =0 C l k : binomial coefficient. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • If | 1 + (∆ t ) λ j | < 1, then � x k : bounded. x 0 s.t. | � • If | 1 + (∆ t ) λ j | > 1, then one can find � x k | → + ∞ (exponentially) as k → + ∞ . x 0 s.t. N � x 0 � = 0 , N 2 � x 0 = 0, • If | 1 + (∆ t ) λ j | = 1 and N � = 0, then for all � x k = (1 + (∆ t ) λ j ) k � x 0 + (1 + (∆ t ) λ j ) k − 1 k ∆ tN � x 0 � goes to infinity as k → + ∞ . • Stability condition | 1 + (∆ t ) λ j | < 1 ⇔ ∆ t < − 2 ℜ λ j | λ j | 2 , holds for ∆ t small enough. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Implicit Euler scheme: x k +1 = ( I − ∆ t ( D + N )) − k � x 0 . � • All the eigenvalues of the matrix ( I − ∆ t ( D + N )) − 1 : of modulus strictly smaller than 1. • ⇒ Implicit Euler scheme: unconditionally stable. • Trapezoidal scheme: x k +1 = ( I − (∆ t ) ( D + N )) − k ( I + (∆ t ) x 0 . ( D + N )) k � � 2 2 • Stability condition: | 1 + (∆ t ) λ j | < | 1 − (∆ t ) λ j | , 2 2 holds for all ∆ t > 0 since ℜ λ j < 0. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • REMARK: Explicit and implicit Euler schemes: of order one; Trapezoidal scheme: of order two. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Runge-Kutta methods: • By far the most popular and powerful general-purpose numerical methods for integrating ODEs. • Idea behind: evaluate f at carefully chosen values of its arguments, t and x , in order to create an accurate approximation (as accurate as a higher-order Taylor expansion) of x ( t + ∆ t ) without evaluating derivatives of f . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Runge-Kutta schemes: derived by matching multivariable Taylor series expansions of f ( t , x ) with the Taylor series expansion of x ( t + ∆ t ). • To find the right values of t and x at which to evaluate f : • Take a Taylor expansion of f evaluated at these (unknown) values; • Match the resulting numerical scheme to a Taylor series expansion of x ( t + ∆ t ) around t . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Generalization of Taylor’s theorem to functions of two variables: THEOREM: • f ( t , x ) ∈ C n +1 ([0 , T ] × R ). Let ( t 0 , x 0 ) ∈ [0 , T ] × R . • There exist t 0 ≤ τ ≤ t , x 0 ≤ ξ ≤ x , s.t. f ( t , x ) = P n ( t , x ) + R n ( t , x ) , • P n ( t , x ): n th Taylor polynomial of f around ( t 0 , x 0 ); • R n ( t , x ): remainder term associated with P n ( t , x ). Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • � � ( t − t 0 ) ∂ f ∂ t ( t 0 , x 0 ) + ( x − x 0 ) ∂ f P n ( t , x ) = f ( t 0 , x 0 ) + ∂ x ( t 0 , x 0 ) � ( t − t 0 ) 2 ∂ 2 f ∂ t 2 ( t 0 , x 0 ) + ( t − t 0 )( x − x 0 ) ∂ 2 f + ∂ t ∂ x ( t 0 , x 0 ) 2 � +( x − x 0 ) 2 ∂ 2 f ∂ x 2 ( t 0 , x 0 ) 2 � 1 � n � ∂ n f C n j ( t − t 0 ) n − j ( x − x 0 ) j . . . + ∂ t n − j ∂ x j ( t 0 , x 0 ) ; n ! j =0 • � n +1 ∂ n +1 f 1 C n +1 ( t − t 0 ) n +1 − j ( x − x 0 ) j R n ( t , x ) = ∂ t n +1 − j ∂ x j ( τ, ξ ) . j ( n + 1)! j =0 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Illustration: obtain a second-order accurate method (truncation error O ((∆ t ) 2 )). • Match x + ∆ tf ( t , x ) + (∆ t ) 2 � ∂ f � + (∆ t ) 3 d 2 ∂ t ( t , x ) + ∂ f ∂ x ( t , x ) f ( t , x ) dt 2 [ f ( τ, x )] 2 6 to x + (∆ t ) f ( t + α 1 , x + β 1 ) , τ ∈ [ t , t + ∆ t ] and α 1 and β 1 : to be found. • Match + (∆ t ) 2 d 2 � ∂ f � f ( t , x ) + (∆ t ) ∂ t ( t , x ) + ∂ f ∂ x ( t , x ) f ( t , x ) dt 2 [ f ( t , x )] 2 6 with f ( t + α 1 , x + β 1 ) at least up to terms of the order of O (∆ t ). Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Multivariable version of Taylor’s theorem to f , ∂ x ( t , x ) + α 2 ∂ 2 f f ( t + α 1 , x + β 1 ) = f ( t , x ) + α 1 ∂ f ∂ t ( t , x ) + β 1 ∂ f 1 ∂ t 2 ( τ, ξ ) 2 + α 1 β 1 ∂ 2 f ∂ t ∂ x ( τ, ξ ) + β 2 ∂ 2 f 1 ∂ x 2 ( τ, ξ ) , 2 t ≤ τ ≤ t + α 1 and x ≤ ξ ≤ x + β 1 . • ⇒ α 1 = ∆ t β 1 = ∆ t and 2 f ( t , x ) . 2 • ⇒ Resulting numerical scheme: explicit midpoint method: the simplest example of a Runge-Kutta method of second order. • Improved Euler method: also another often-used Runge-Kutta method. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • General Runge-Kutta method: � m x k +1 = x k + ∆ t c i f ( t i , k , x i , k ) , i =1 m : number of terms in the method. • Each t i , k denotes a point in [ t k , t k +1 ]. • Second argument x i , k ≈ x ( t i , k ) can be viewed as an approximation to the solution at the point t i , k . • To construct an n th order Runge-Kutta method, we need to take at least m ≥ n terms. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Best-known Runge-Kutta method: fourth-order Runge-Kutta method , which uses four evaluations of f during each step.  κ 1 := f ( t k , x k ) ,       2 , x k + ∆ t  κ 2 := f ( t k + ∆ t 2 κ 1 ) ,      2 , x k + ∆ t κ 3 := f ( t k + ∆ t 2 κ 2 ) ,    κ 4 := f ( t k +1 , x k + ∆ t κ 3 ) ,         x k +1 = x k + (∆ t )  ( κ 1 + 2 κ 2 + 2 κ 3 + κ 4 ) . 6 • Values of f at the midpoint in time: given four times as much weight as values at the endpoints t k and t k +1 (similar to Simpson’s rule from numerical integration). Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Construction of Runge-Kutta methods: • Construct Runge-Kutta methods by generalizing collocation methods. • Discuss their consistency, stability, and order. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Collocation methods: • P m : space of real polynomials of degree ≤ m . • Interpolating polynomial: • Given a set of m distinct quadrature points c 1 < c 2 < . . . < c m in R , and corresponding data g 1 , . . . , g m ; • There exists a unique polynomial, P ( t ) ∈ P m − 1 s.t. P ( c i ) = g i , i = 1 , . . . , m . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • DEFINITION: • Define the i th Lagrange interpolating polynomial l i ( t ), i = 1 , . . . , m , for the set of quadrature points { c j } by m � t − c j l i ( t ) := . c i − c j j � = i , j =1 • Set of Lagrange interpolating polynomials: form a basis of P m − 1 ; • Interpolating polynomial P corresponding to the data { g j } given by � m P ( t ) := g i l i ( t ) . i =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Consider a smooth function g on [0 , 1]. • Approximate the integral of g on [0 , 1] by exactly integrating the Lagrange interpolating polynomial of order m − 1 based on m quadrature points 0 ≤ c 1 < c 2 < . . . < c m ≤ 1. • Data: values of g at the quadrature points g i = g ( c i ), i = 1 , . . . , m . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Define the weights � 1 b i = l i ( s ) ds . 0 • Quadrature formula: � 1 � 1 m m � � g ( s ) ds ≈ g i l i ( s ) ds = b i g ( c i ) . 0 0 i =1 i =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • f : smooth function on [0 , T ]; t k = k ∆ t for k = 0 , . . . , K = T / (∆ t ): discretization points in [0 , T ]. � t k +1 • f ( s ) ds can be approximated by t k � t k +1 � 1 � m f ( s ) ds = (∆ t ) f ( t k + ∆ t τ ) d τ ≈ (∆ t ) b i f ( t k + (∆ t ) c i ) . t k 0 i =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • x : polynomial of degree m satisfying  x (0) = x 0 ,  dx  dt ( c i ∆ t ) = F i , F i ∈ R , i = 1 , . . . , m . • Lagrange interpolation formula ⇒ for t in the first time-step interval [0 , ∆ t ], � m dx F i l i ( t dt ( t ) = ∆ t ) . i =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Integrating over the intervals [0 , c i ∆ t ] ⇒ � c i � m � m x ( c i ∆ t ) = x 0 + (∆ t ) F j l j ( s ) ds = x 0 + (∆ t ) a ij F j , 0 j =1 j =1 for i = 1 , . . . , m , with � c i a ij := l j ( s ) ds . 0 • Integrating over [0 , ∆ t ] ⇒ � 1 � m � m x (∆ t ) = x 0 + (∆ t ) F i l i ( s ) ds = x 0 + (∆ t ) b i F i . 0 i =1 i =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Writing dx / dt = f ( x ( t )), on the first time step interval [0 , ∆ t ],  m �    F i = f ( x 0 + (∆ t ) a ij F j ) , i = 1 , . . . , m ,    j =1  � m     x (∆ t ) = x 0 + (∆ t ) b i F i .  i =1 • Similarly, we have on [ t k , t k +1 ]  m �    F i , k = f ( x ( t k ) + (∆ t ) a ij F j , k ) , i = 1 , . . . , m ,    j =1  � m    x ( t k +1 ) = x ( t k ) + (∆ t ) b i F i , k .   i =1 • In the collocation method : one first solves the coupled nonlinear system to obtain F i , k , i = 1 , . . . , m , and then computes x ( t k +1 ) from x ( t k ). Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • REMARK: • m � t l − 1 = c l − 1 l i ( t ) , t ∈ [0 , 1] , l = 1 , . . . , m , i i =1 • ⇒ m � = 1 b i c l − 1 l , l = 1 , . . . , m , i i =1 and m � = c l a ij c l − 1 i i , l = 1 , . . . , m . l , j j =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Runge-Kutta methods as generalized collocation methods • In the collocation method, the coefficients b i and a ij : defined by certain integrals of the Lagrange interpolating polynomials associated with a chosen set of quadrature nodes c i , i = 1 , . . . , m . • Natural generalization of collocation methods: obtained by allowing the coefficients c i , b i , and a ij to take arbitrary values, not necessary related to quadrature formulas. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • No longer assume the c i to be distinct. • However, assume that m � c i = a ij , i = 1 , . . . , m . j =1 • ⇒ Class of Runge-Kutta methods for solving the ODE,  � m  F i , k = f ( t i , k , x k + (∆ t )   a ij F j , k ) ,    j =1  m �   x k +1 = x k + (∆ t )  b i F i , k ,   i =1 t i , k = t k + c i ∆ t , or equivalently,  m �   x i , k = x k + (∆ t )  a ij f ( t j , k , x j , k ) ,    j =1  � m   x k +1 = x k + (∆ t )  b i f ( t i , k , x i , k ) .   i =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Let κ j := f ( t + c j ∆ t , x j ); • Define Φ by  m �    x i = x + (∆ t ) a ij κ j ,    j =1  � m    Φ( t , x , ∆ t ) = b i f ( t + c i ∆ t , x i ) .   i =1 • ⇒ One step method. • If a ij = 0 for j ≥ i ⇒ scheme: explicit. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • EXAMPLES: • Explicit Euler’s method and Trapezoidal scheme: Runge-Kutta methods. • Explicit Euler’s method: m = 1 , b 1 = 1 , a 11 = 0. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Trapezoidal scheme: m = 2 , b 1 = b 2 = 1 / 2 , a 11 = a 12 = 0 , a 21 = a 22 = 1 / 2. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Fourth-order Runge-Kutta method: m = 4 , c 1 = 0 , c 2 = c 3 = 1 / 2 , c 4 = 1 , b 1 = 1 / 6 , b 2 = b 3 = 1 / 3 , b 4 = 1 / 6 , a 21 = a 32 = 1 / 2 , a 43 = 1 , and all the other a ij entries are zero. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Consistency, stability, convergence, and order of Runge-Kutta methods • Runge-Kutta scheme: consistent iff � m b j = 1 . j =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Stability: • | A | = ( | a ij | ) m i , j =1 . • Spectral radius ρ ( | A | ) of the matrix | A | : ρ ( | A | ) := max {| λ j | , λ j : eigenvalue of | A |} . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • THEOREM: • C f : Lipschitz constant for f . • Suppose (∆ t ) C f ρ ( | A | ) < 1 . • Then the Runge-Kutta method : stable. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • PROOF: • � � m � Φ( t , x , ∆ t ) − Φ( t , y , ∆ t ) = b i f ( t + c i ∆ t , x i ) − f ( t + c i ∆ t , y i ) , i =1 with m � x i = x + (∆ t ) a ij f ( t + c j ∆ t , x j ) , j =1 and m � y i = y + (∆ t ) a ij f ( t + c j ∆ t , y j ) . j =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • ⇒ � � m � x i − y i = x − y + (∆ t ) a ij f ( t + c j ∆ t , x j ) − f ( t + c j ∆ t , y j ) . j =1 • ⇒ For i = 1 , . . . , m , � m | x i − y i | ≤ | x − y | + (∆ t ) C f | a ij || x j − y j | . j =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • X and Y :     | x 1 − y 1 | | x − y |  .   .  . . X = and Y =  .    . . | x m − y m | | x − y | • X ≤ Y + (∆ t ) C f | A | X , ⇒ X ≤ ( I − (∆ t ) C f | A | ) − 1 Y , provided that (∆ t ) C f ρ ( | A | ) < 1 . • ⇒ stability of the Runge-Kutta scheme. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Dahlquist-Lax equivalence theorem ⇒ Runge-Kutta scheme: convergent provided that � m j =1 b j = 1 and (∆ t ) C f ρ ( | A | ) < 1 hold. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Order of the Runge-Kutta scheme: compute the order as ∆ t → 0 of the truncation error T k (∆ t ) = x ( t k +1 ) − x ( t k ) − Φ( t k , x ( t k ) , ∆ t ) . ∆ t • Write � m � m T k (∆ t ) = x ( t k +1 ) − x ( t k ) − b i f ( t k + c i ∆ t , x ( t k ) + ∆ t a ij κ j ) . ∆ t i =1 j =1 • Suppose that f : smooth enough ⇒ m � f ( t k + c i ∆ t , x ( t k ) + ∆ t a ij κ j ) j =1 � � � c i ∂ f a ij κ j ) ∂ f = f ( t k , x ( t k )) + ∆ t ∂ t ( t k , x ( t k )) + ( ∂ x ( t k , x ( t k )) j =1 + O ((∆ t ) 2 ) . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • � � a ij κ j = ( a ij ) f ( t k , x ( t k )) + O (∆ t ) = c i f ( t k , x ( t k )) + O (∆ t ) . j =1 j =1 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • ⇒ � m f ( t k + c i ∆ t , x ( t k ) + ∆ t a ij κ j ) j =1 � ∂ f � ∂ t ( t k , x ( t k )) + ∂ f = f ( t k , x ( t k )) + ∆ tc i ∂ x ( t k , x ( t k )) f ( t k , x ( t k )) + O ((∆ t ) 2 ) . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • THEOREM: • Assume that f : smooth enough. • Then the Runge-Kutta scheme: of order 2 provided that the conditions � m b j = 1 j =1 and m � b i c i = 1 2 i =1 hold. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Higher-order Taylor expansions ⇒ • THEOREM: • Assume that f : smooth enough. • Then the Runge-Kutta scheme: of order 3 provided that the conditions m � b j = 1 , j =1 m � b i c i = 1 2 , i =1 and � m � m � m i = 1 b i a ij c j = 1 b i c 2 3 , 6 i =1 i =1 j =1 hold. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Of Order 4 provided that in addition m m m m m � � � � � i = 1 b i c i a ij c j = 1 j = 1 b i c 3 b i c i a ij c 2 4 , 8 , 12 , i =1 i =1 j =1 i =1 j =1 � m � m � m b i a ij a jl c l = 1 24 i =1 j =1 l =1 hold. • The (fourth-order) Runge-Kutta scheme: of order 4. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Multi-step methods • Runge-Kutta methods: improvement over Euler’s methods in terms of accuracy, but achieved by investing additional computational effort. • The fourth-order Runge-Kutta method involves four function evaluations per step. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • For comparison, by considering three consecutive points t k − 1 , t k , t k +1 , integrating the differential equation between t k − 1 and t k +1 , and applying Simpson’s rule to approximate the resulting integral yields � t k +1 x ( t k +1 ) = x ( t k − 1 ) + f ( s , x ( s )) ds t k − 1 � � ≈ x ( t k − 1 ) + (∆ t ) f ( t k − 1 , x ( t k − 1 )) + 4 f ( t k , x ( t k )) + f ( t k +1 , x ( t k +1 )) , 3 ⇒ � � x k +1 = x k − 1 + (∆ t ) f ( t k − 1 , x k − 1 ) + 4 f ( t k , x k ) + f ( t k +1 , x k +1 ) . 3 • Need two preceding values, x k and x k − 1 in order to calculate x k +1 : two-step method . • In contrast with the one-step methods: only a single value of x k required to compute the next approximation x k +1 . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • General n -step method: n n � � α j x k + j = (∆ t ) β j f ( t k + j , x k + j ) , j =0 j =0 α j and β j : real constants and α n � = 0. • If β n = 0, then x k + n : obtained explicitly from previous values of x j and f ( t j , x j ) ⇒ n -step method: explicit . Otherwise, the n -step method: implicit . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • EXAMPLE: (i) Two-step Adams-Bashforth method: explicit two-step method � � x k +2 = x k +1 + (∆ t ) 3 f ( t k +1 , x k +1 ) − f ( t k , x k ) ; 2 (ii) Three-step Adams-Bashforth method: explicit three-step method � � x k +3 = x k +2 +(∆ t ) 23 f ( t k +2 , x k +2 ) − 16 f ( t k +1 , x k +1 )+ f ( t k , x k ) ; 12 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs (iii) Four-step Adams-Bashforth method: explicit four-step method � x k +4 = x k +3 + (∆ t ) 55 f ( t k +3 , x k +3 ) − 59 f ( t k +2 , x k +2 ) 24 � +37 f ( t k +1 , x k +1 ) − 9 f ( t k , x k ) ; (iv) Two-step Adams-Moulton method: implicit two-step method � � x k +2 = x k +1 + (∆ t ) 5 f ( t k +2 , x k +2 ) + 8 f ( t k +1 , x k +1 ) + f ( t k , x k ) ; 12 (v) Three-step Adams-Moulton method: implicit three-step method � � x k +3 = x k +2 +(∆ t ) 9 f ( t k +3 , x k +3 )+19 f ( t k +2 , x k +2 ) − 5 f ( t k +1 , x k +1 ) − 9 f ( t k , x k ) . 24 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Construction of linear multi-step methods • Suppose that x k , k ∈ N : sequence of real numbers. • Shift operator E , forward difference operator ∆ + and backward difference operator ∆ − : E : x k �→ x k +1 , ∆ + : x k �→ x k +1 − x k , ∆ − : x k �→ x k − x k − 1 . • ∆ + = E − I and ∆ − = I − E − 1 ⇒ for any n ∈ N , � n ( E − I ) n = ( − 1) j C n j E n − j j =0 and n � ( I − E − 1 ) n = ( − 1) j C n j E − j . j =0 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • ⇒ � n + x k = ∆ n ( − 1) j C n j x k + n − j j =0 and � n − x k = ∆ n ( − 1) j C n j x k − j . j =0 Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • y ( t ) ∈ C ∞ ( R ); t k = k ∆ t , ∆ t > 0. • Taylor series ⇒ for any s ∈ N , � + ∞ � � � � 1 l !( s ∆ t ∂ e s (∆ t ) ∂ E s y ( t k ) = y ( t k + s ∆ t ) = ∂ t ) l y ∂ t y ( t k ) = ( t k ) , l =0 • ⇒ E s = e s (∆ t ) ∂ ∂ t . • Formally, (∆ t ) ∂ ∂ t = ln E = − ln( I − ∆ − ) = ∆ − + 1 − + 1 2∆ 2 3∆ 3 − + . . . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • x ( t ): solution of ODE: � � ∆ − + 1 − + 1 2∆ 2 3∆ 3 (∆ t ) f ( t k , x ( t k )) = − + . . . x ( t k ) . • Successive truncation of the infinite series ⇒ x k − x k − 1 = (∆ t ) f ( t k , x k ) , 3 2 x k − 2 x k − 1 + 1 2 x k − 2 = (∆ t ) f ( t k , x k ) , 11 6 x k − 3 x k − 1 + 3 2 x k − 2 − 1 3 x k − 3 = (∆ t ) f ( t k , x k ) , and so on. • Class of implicit multi-step methods: backward differentiation formulas. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Similarly, E − 1 ((∆ t ) ∂ ∂ t ) = (∆ t ) ∂ ∂ t E − 1 = − ( I − ∆ − ) ln( I − ∆ − ) . • ⇒ ((∆ t ) ∂ ∂ t ) = − E ( I − ∆ − ) ln( I − ∆ − ) = − ( I − ∆ − ) ln( I − ∆ − ) E . • ⇒ � � ∆ − − 1 − − 1 2∆ 2 6∆ 3 (∆ t ) f ( t k , x ( t k )) = − + . . . x ( t k +1 ) . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Successive truncation of the infinite series ⇒ explicit numerical schemes: x k +1 − x k = (∆ t ) f ( t k , x k ) , 1 2 x k +1 − 1 2 x k − 1 = (∆ t ) f ( t k , x k ) , 1 3 x k +1 + 1 2 x k − x k − 1 + 1 6 x k − 2 = (∆ t ) f ( t k , x k ) , . . . • The first of these numerical scheme: explicit Euler method, while the second: explicit mid-point method. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Construct further classes of multi-step methods: • For y ∈ C ∞ , � t k D − 1 y ( t k ) = y ( t 0 ) + y ( s ) ds , t 0 and � t k +1 ( E − I ) D − 1 y ( t k ) = y ( s ) ds . t k • ( E − I ) D − 1 = ∆ + D − 1 = E ∆ − D − 1 = (∆ t ) E ∆ − ((∆ t ) D ) − 1 , Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • ⇒ � � − 1 . ( E − I ) D − 1 = − (∆ t ) E ∆ − ln( I − ∆ − ) Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • � (∆ t ) DE − 1 � − 1 . ( E − I ) D − 1 = E ∆ − D − 1 = ∆ − ED − 1 = ∆ − ( DE − 1 ) − 1 = (∆ t )∆ − • ⇒ � � − 1 ( E − I ) D − 1 = − (∆ t )∆ − ( I − ∆ − ) ln( I − ∆ − ) . Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • � t k +1 f ( s , x ( s )) ds = ( E − I ) D − 1 f ( t k , x ( t k )) , x ( t k +1 ) − x ( t k ) = t k • ⇒  � � − 1 f ( t k , x ( t k ))  − (∆ t )∆ − ( I − ∆ − ) ln( I − ∆ − ) x ( t k +1 ) − x ( t k ) = � � − 1 f ( t k , x ( t k )) .  − (∆ t ) E ∆ − ln( I − ∆ − ) Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Expand ln( I − ∆ − ) into a Taylor series on the right-hand side ⇒ � � I + 1 2∆ − + 5 − + 3 12∆ 2 8∆ 3 x ( t k +1 ) − x ( t k ) = (∆ t ) − + . . . f ( t k , x ( t k )) and � � I − 1 2∆ − − 1 − − 1 12∆ 2 24∆ 3 x ( t k +1 ) − x ( t k ) = (∆ t ) − + . . . f ( t k +1 , x ( t k +1 )) . • Successive truncations ⇒ families of (explicit) Adams-Bashforth methods and of (implicit) Adams-Moulton methods. Numerical methods for ODEs Habib Ammari

Numerical solution of ODEs • Consistency, stability, and convergence • Introduce the concepts of consistency, stability, and convergence for analyzing linear multi-step methods. Numerical methods for ODEs Habib Ammari