Validated Solution of Initial Value Problems for ODEs with Interval Parameters Youdong Lin and Mark A. Stadtherr Department of Chemical and Biomolecular Engineering, University of Notre Dame, Notre Dame, IN, USA 2nd NSF Workshop on Reliable Engineering Computing, Georgia Tech University, Savannah, GA, February 22–24, 2006
Outline • Motivating Problem: Bioreactor Simulation • Background • Overview of Method • Examples and Results – Lotka-Volterra Problem – Lorenz Problem – Double Pendulum Problem – Bioreactor Problem • Concluding Remarks 2
Motivating Example – Bioreactor Simulation • In a bioreactor, microbial growth may be described by ˙ X = ( µ − αD ) X S = D ( S i − S ) − kµX, ˙ where X and S are concentrations of biomass and substrate, respectively. • The growth rate µ may be given by µ m S µ = (Monod Law) K S + S or µ m S µ = (Haldane Law) K S + S + K I S 2 3
Motivating Example – Bioreactor Simulation • Problem data Value Units Value Units day − 1 α µ m 0.5 - [1.19, 1.21] k K S 10.53 g S/ g X [7.09, 7.11] g S/l day − 1 (g S/l) − 1 D K I 0.36 [0.49, 0.51] S i X 0 5.7 g S/l [0.82, 0.84] g X/l S 0 0.80 g S/l • Three parameters ( µ m , K S and K I ) and one initial state ( X 0 ) are uncertain and given by intervals. • Problem: Determine a validated enclosure of all possible solutions to this ODE system. • Issue: Standard tools for validated solution of ODEs are designed to deal with interval-valued initial states, not interval-valued parameters. 4
Problem Definition • Consider ˙ x = f ( x , θ ) x ( t 0 ) = x 0 ∈ X 0 θ ∈ Θ x = state vector ( m variables) θ = parameter vector ( p parameters) X 0 = interval enclosure of x 0 Θ = interval enclosure of θ • Consider time steps h j = t j +1 − t j , j = 0 , . . . , N − 1 • Notation: x ( t ; t j , x j , θ ) denotes a solution of ˙ x = f ( x , θ ) for the initial condition x = x j at t = t j and x ( t ; t j , X j , Θ ) is the set of solutions x ( t ; t j , X j , Θ ) = { x ( t ; t j , x j , θ ) | x j ∈ X j , θ ∈ Θ } • Problem: Determine enclosures X j of the state variables at each time t j , j = 1 , . . . , N , such that x ( t j ; t 0 , X 0 , Θ ) ⊆ X j 5
Background – Interval Taylor Series • In an Taylor series expansion of x ( t ) with respect to t , the coefficients can be obtained recursively in terms of ˙ x ( t ) = f ( x , θ ) using f [0] = x f [1] = f ( x , θ ) � � ∂ f [ i − 1] 1 f [ i ] = ( x , θ ) , i ≥ 2 . f i ∂ x • Values of these coefficients can be easily generated using automatic differentiation techniques. • For an interval Taylor series (ITS), the coefficients F [ i ] are interval enclosures of f [ i ] . 6
Background – Taylor Models • Taylor Model T f = ( p f , R f ) : Bounds a function f ( x ) over X using a q -th order Taylor polynomial p f and an interval remainder bound R f , usually from a truncated Taylor series. q i ! [( x − x 0 ) · ▽ ] i f ( x 0 ) 1 � p f = i =0 ( q +1)! [( x − x 0 ) · ▽ ] q +1 F [ x 0 + ( x − x 0 ) ζ ] 1 R f = where, x 0 ∈ X ; ζ ∈ [0 , 1] [ g · ▽ ] k = ∂ k j 1 ! ··· j m ! g j 1 k ! 1 · · · g j m � m ∂x j 1 1 ··· ∂x jm m j 1+ ··· + jm = k 0 ≤ j 1 , ··· ,jm ≤ k • Store and operate on coefficients of p f only. Floating point errors are accumulated in R f . 7
Background – Taylor Model Operations • Taylor model of f ± g f ± g ∈ ( p f , R f ) ± ( p g , R g ) = ( p f ± p g , R f ± R g ) T f ± g = ( p f ± g , R f ± g ) = ( p f ± p g , R f ± R g ) • Taylor model of f × g f × g ∈ ( p f , R f ) × ( p g , R g ) ⊆ p f × p g + p f × R g + p g × R f + R f × R g Split p f × p g into q -th order part p f × g and higher-order terms p e . Then T f × g = ( p f × g , R f × g ) R f × g = B ( p e ) + B ( p f ) × R g + B ( p g ) × R f + R f × R g B ( p ) indicates an interval bound on the function p . • Reciprocal operation and intrinsic functions can also be defined. 8
Background – Interval IVPs • Consider standard ODE system (non-parametric) x = f ( x ) ˙ x ( t 0 ) = x 0 ∈ X 0 • “Standard” approach (step j + 1 ): Assuming X j is known, then – Phase 1: Compute a coarse enclosure ˜ X j and prove existance and uniqueness. Use fixed point iteration with Picard operator using high-order interval Taylor series. – Phase 2: Refine the coarse enclosure to obtain X j +1 . Use high-order interval Taylor series with Taylor coefficients bounded using mean value theorem. Reduce wrapping effect using QR-factorization approach. • Implementations include AWA and VNODE. 9
Method for Parametric ODEs • Consider again parametric ODE system x = f ( x , θ ) ˙ x ( t 0 ) = x 0 ∈ X 0 θ ∈ Θ • To apply standard methods, can treat parameters as additional state variables with zero derivative (Lohner, 1988) • Our method for parametric ODE system: Assuming X j is known, then – Phase 1: Same as “standard” approach. Compute a coarse enclosure ˜ X j and prove existance and uniqueness. Use fixed point iteration with Picard operator using high-order interval Taylor series. – Phase 2: Refine the coarse enclosure to obtain X j +1 . Use Taylor models in terms of the uncertain parameters and initial states. • Implemented in VSPODE (Validating Solver for Parametric ODEs) (Lin and Stadtherr, 2005). 10
Method for Phase 2 • Represent uncertain initial states and parameters using Taylor models T x 0 and T θ , with components T x i 0 = ( m ( X i 0 ) + ( x i 0 − m ( X i 0 )) , [0 , 0]) , i = 1 , · · · , m T θ i = ( m (Θ i ) + ( θ i − m (Θ i )) , [0 , 0]) , i = 1 , · · · , p. • Compute Taylor models T f [ i ] for the interval Taylor series coefficients using Taylor model operations and obtain the polynomial part of T x j +1 . • Determine the remainder bound of T x j +1 by the mean value theorem and reduce the wrapping effect using a QR factorization approach, where the remainder is represented by R x j +1 = A j +1 V j +1 . • Compute the enclosure X j +1 = B ( T x j +1 ) by bounding over X 0 and Θ . 11
Examples and Results • Computations done with Intel Pentium 4 3.2GHz CPU on a Linux workstation. • For comparsions, VNODE was used, with interval parameters treated as additional state variables • VSPODE run using → q = 5 (order of Taylor model) → k = 17 (order of interval Taylor series) → QR • VNODE run using → k = 17 order interval Hermite-Obreschkoff → QR 12
Example 1. Lotka-Volterra Problem • ODE model is x 1 = θ 1 x 1 (1 − x 2 ) ˙ x 2 = θ 2 x 2 ( x 1 − 1) ˙ x 0 = (1 . 2 , 1 . 1) T θ 1 ∈ [2 . 99 , 3 . 01] θ 2 ∈ [0 . 99 , 1 . 01] • Integrate from t 0 = 0 to t N = 10 . • Constant step size of h = 0 . 1 used in both VSPODE and VNODE. 13
Example 1. Lotka-Volterra Problem 1.5 1.4 ← y 1, VNODE 1.3 1.2 1.1 ← y 2, VSPODE y 1 /y 2 1 ← y 1, VSPODE 0.9 0.8 ← y 2, VNODE 0.7 0.6 0.5 0 1 2 3 4 5 6 7 8 9 10 t (Eventual breakdown of VSPODE at t = 31 . 8 ) 14
Example 1. Lotka-Volterra Problem • To allow VNODE to integrate further: – Parameters intervals can be subdivided into equal-sized subintervals. – Apply VNODE to each parameter subinterval. – Final enclosure is the union of enclosures determined from each subinterval. • VNODE-NN indicates use of NN parameter subintervals. 15
Example 1. Lotka-Volterra Problem Final Enclosure ( t = 10 ) Method Width CPU time (s) VSPODE [ 1.120873, 1.173607 ] 0.052734 0.59 [ 0.875994, 0.893471 ] 0.017477 VNODE–16 [ 1.110859, 1.182814 ] 0.071955 1.42 [ 0.872528, 0.898407 ] 0.025879 VNODE–36 [ 1.116350, 1.177431 ] 0.061081 3.14 [ 0.874924, 0.895612 ] 0.020688 VNODE–64 [ 1.118151, 1.175692 ] 0.057541 5.59 [ 0.875651, 0.894736 ] 0.019085 VNODE–100 [ 1.118999, 1.174881 ] 0.055882 8.68 [ 0.875975, 0.894337 ] 0.018362 16
Example 2. Lorenz Problem • ODE model is x 1 = θ 1 ( x 2 − x 1 ) ˙ x 2 = x 1 ( θ 2 − x 3 ) − x 2 ˙ x 3 = x 1 x 2 − θ 3 x 3 ˙ x 0 = (10 , 10 , 10) T θ 1 ∈ 10 + [ − 0 . 01 , 0 . 01] θ 2 ∈ 28 + [ − 0 . 01 , 0 . 01] θ 3 ∈ 8 / 3 + [ − 0 . 01 , 0 . 01] • Integrate from t 0 = 0 to t N = 2 . • Constant step size of h = 0 . 01 used in both VSPODE and VNODE. 17
Example 2. Lorenz Problem 50 ← y 1, VNODE 40 ← y 3, VSPODE 30 20 y 1 /y 2 /y 3 ← y 2, VNODE 10 ← y 3, VNODE 0 ← y 1, VSPODE −10 ← y 2, VSPODE −20 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 t (Eventual breakdown of VSPODE at t = 2 . 8 ) 18
Example 2. Lorenz Problem Final Enclosure ( t = 2 ) Method Width CPU time (s) VSPODE [ -0.582033, -0.342358 ] 0.2397 2.66 [ -0.769513, -0.369357 ] 0.4002 [ 14.633803, 14.737535 ] 0.1037 VNODE–125 [ -8.663336, 7.988072 ] 16.6514 33.7 [ -10.060512, 8.797511 ] 18.8580 [ 9.031894, 21.106684 ] 12.0748 VNODE–512 [ -0.920184, 0.041287 ] 0.9615 141.5 [ -1.321734, 0.245595 ] 1.5673 [ 14.352124, 15.010891 ] 0.6588 VNODE–1000 [ -0.770156, -0.136139 ] 0.6340 263.1 [ -1.077794, -0.036474 ] 1.0413 [ 14.502030, 14.869122 ] 0.3671 19
Example 3. Double Pendulum Problem L 1 m 1 � 1 L 2 m 2 � 2 m 1 = m 2 = 1 kg L 1 = L 2 = 1 m 20
Recommend
More recommend
