Modelling Biochemical Reaction Networks Lecture 7: Numerical - PowerPoint PPT Presentation

Modelling Biochemical Reaction Networks Lecture 7: Numerical integration of ordinary differential equations Marc R. Roussel Department of Chemistry and Biochemistry

Recommended reading ◮ Fall, Marland, Wagner and Tyson, section 1.4.1

Ordinary differential equation initial value problems ◮ An ordinary differential equation (ODE) is an equation of the form d x dt = f ( x , t ) x is a vector describing the state of a system (e.g. the concentrations) and f ( x , t ) is a vector-valued function. ◮ In an initial value problem (IVP), we are provided with x at some initial time t 0 and want to get x for t > t 0 . ◮ A solution of an ODE IVP is a function x ( t ) satisfying the equation and the initial condition. ◮ Very few ODEs have analytic solutions, so we generally need to use approximate numerical integration methods.

Numerical integration of ordinary differential equations ◮ Basic idea: Approximate the derivative as d x dt ≈ ∆ x ∆ t = f ( x , t ) ◮ The continuous solution x ( t ) is replaced by values of x at a discrete set of times: ( t i , x i ) ( i = 1 , 2 , . . . ) with x 0 = x (0). ◮ Need to decide 1. what exactly we mean by ∆ x ; 2. how to choose ∆ t ; and 3. at what x and t we’re going to evaluate x .

Euler’s method ◮ Simplest possible method: ∆ x = x i +1 − x i ∆ t = t i +1 − t i is usually constant f evaluated at x i and t i ∆ x ∆ t = x i +1 − x i = f ( x i , t i ) ∆ t ∴ x i +1 = x i + f ( x i , t i )∆ t ◮ xppaut example: dx dt = − x , x (0) = 1

Runge-Kutta methods ◮ The problem with the Euler method is that it uses the derivative at a single point and extrapolates from there. ◮ Many numerical methods calculate the derivative using one or more intermediate points in order to obtain more refined estimates of the average derivative over one time step. ◮ Runge-Kutta methods impose the condition that the Taylor series of x ( t + h ) in powers of h = ∆ t should match its approximation by the numerical method. ◮ Runge-Kutta methods involve multiple stages of computation in which we use rates computed at previous points to estimate the position of intermediate points.

Two-stage Runge-Kutta method ◮ For simplicity, consider a scalar differential equation. ◮ Take one step of size h = ∆ t . ◮ Use one intermediate point: ◮ t intermed = t i + ch where c is a coefficient between 0 and 1 ◮ Use Euler’s method to estimate x intermed ≈ x i + chf ( x i , t i ) ◮ Blend the two derivatives at ( t i , x i ) and ( t intermed , x intermed ) to obtain the estimate of x i +1 : x i +1 = x i + h [ a 1 f ( x i , t i ) + a 2 f ( x i + chf ( x i , t i ) , t i + ch )] ◮ We need to choose values for the coefficients a 1 , a 2 and c .

Two-stage Runge-Kutta method ◮ Taylor expansion in powers of h of x i +1 = x ( t i + h ): x ( t i + h ) = x ( t i ) + hx ′ ( t i ) + h 2 2 x ′′ ( t i ) + O ( h 3 ) = x i + hf ( x i , t i ) + h 2 � d + O ( h 3 ) � dt x ′ ( t ) � 2 � ( t i , x i ) = x i + hf ( x i , t i ) + h 2 � d + O ( h 3 ) � dt f ( x , t ) � 2 � ( t i , x i ) = x i + hf ( x i , t i ) + h 2 � ∂ f dx dt + ∂ f � + O ( h 3 ) 2 ∂ x ∂ t ( t i , x i ) � � = x i + hf ( x i , t i ) + h 2 � � f ( x i , t i ) ∂ f + ∂ f � � � � 2 ∂ x ∂ t � � ( t i , x i ) ( t i , x i ) + O ( h 3 )

Two-stage Runge-Kutta method ◮ Taylor expansion in powers of h of the right-hand side (rhs) x i + h [ a 1 f ( x i , t i ) + a 2 f ( x i + chf ( x i , t i ) , t i + ch )]: � rhs = x i + h a 1 f ( x i , t i ) � �� � � f ( x i , t i ) + chf ( x i , t i ) ∂ f + ch ∂ f � � + a 2 � � ∂ x ∂ t � � ( t i , x i ) ( t i , x i ) + O ( h 3 ) � � � � f ( x i , t i ) ∂ f + ∂ f = x i + h ( a 1 + a 2 ) f ( x i , t i ) + a 2 ch 2 � � � � ∂ x ∂ t � � ( t i , x i ) ( t i , x i ) + O ( h 3 )

Two-stage Runge-Kutta method ◮ Now compare the two Taylor expansions: � � x i + hf ( x i , t i ) + h 2 � � f ( x i , t i ) ∂ f + ∂ f � � � � 2 ∂ x ∂ t � � ( t i , x i ) ( t i , x i ) � � � � f ( x i , t i ) ∂ f + ∂ f = x i + h ( a 1 + a 2 ) f ( x i , t i )+ a 2 ch 2 � � � � ∂ x ∂ t � � ( t i , x i ) ( t i , x i ) ◮ They are the same if a 1 + a 2 = 1 and a 2 c = 1 2 . ◮ Can write everything in terms of one parameter, say a 2 : 1 a 1 = 1 − a 2 and c = 2 a 2 ( a 2 � = 0).

Two-stage Runge-Kutta method ◮ We get a family of two-stage (and second-order, i.e. with an error O ( h 3 )) Runge-Kutta methods parameterized by a 2 : � x i +1 = x i + h (1 − a 2 ) f ( x i , t i ) � x i + 1 hf ( x i , t i ) , t i + 1 �� + a 2 f h 2 a 2 2 a 2 ◮ Examples: � x i + h 2 f ( x i , t i ) , t i + h � a 2 = 1: x i +1 = x i + hf 2 (Midpoint rule) a 2 = 1 2 : x i +1 = � 1 2 f ( x i , t i ) + 1 � x i + h 2 f ( x i + hf ( x i , t i ) , t i + h ) (Improved Euler)

Two-stage Runge-Kutta method dx Example: dt = x (1 − x ) , x (0) = 0 . 1 (a version of the logistic equation) ◮ For this ODE, f ( x , t ) = x (1 − x ). ◮ The general two-stage Runge-Kutta method is � x i +1 = x i + h (1 − a 2 ) f ( x i , t i ) � x i + 1 hf ( x i , t i ) , t i + 1 �� + a 2 f h 2 a 2 2 a 2 ◮ This is a map, a rule for calculating x i +1 from x i . ◮ Maps can be implemented in xppaut .