Solving Systems of Ordinary Differential Equations

Systems of Ordinary Differential Equations Outline Solving Systems of Ordinary Differential Equations • Review basics of numerical solutions of ordinary differential equations • Error control in simple methods discussed previously • Systems of differential equations • Reducing higher-order equations into a system of first-order equations • Writing and using software for solving systems of ordinary differential equations 2 Schedule for April and May Monday Tuesday Wednesday Thursday Friday 1 4 14 15 21 22 28 PA 6 5 PA 7 29 6 2 3 Spring Break 16 Quiz 7 & 17 PA5 23 24 30 Quiz 8 7 Prog Exam • Solve initial value problem, dy/dx = f(x,y) (known equation) with y(x0) = y0 – Use a finite difference grid: xi+1 – xi = hi – Replace derivative by finite-difference approximation: dy/dx  (yi+1 – yi) / (xi+1 – xi) = (yi+1 – yi) / hi – Derive a formula to compute favg the average value of f(x,y) between xi and xi+1 – Replace dy/dx = f(x,y) by (yi+1 – yi) / hi = favg – Repeatedly compute yi+1 = yi + hi favg 16 35 May 1 8 Review Numerical Approach 2 9 12 Final 3 4 Review Notation and Order • • • • • • • • Review Simple Methods xi is independent variable yi is numerical solution at x = xi fi is derivative found from xi yi: fi = f(xi, yi) y(xi) is the exact value of y at x = xi f(xi,y(xi)) is the exact derivative e1 = y(x1) – y1 = local truncation error Ej = y(xj) – yj is global truncation error If e is O(hn), then E is O(hn-1) 5 • Euler: yi+1 = yi + hifi = yi + hi f(xi, yi) • Huen’s method y 0 i 1 yi 1  yi  hi 1 f ( xi , yi ) xi 1  xi  hi 1   h y  yi01  hi 1 f ( xi 1 , yi01 )  yi  i 1 f ( xi , yi )  f ( xi 1 , yi01 )  i 2 2 • Modified Euler method h  yi  1  yi   i 1  f ( xi , yi ) 2  2  yi 1  yi  hi 1 f ( xi  1 , yi  1 ) 2 xi  1  xi  2 hi 1 2 2 6 1 Systems of Ordinary Differential Equations Review 4th Order Runge-Kutta Error (at Error x = 1) versus stepforsize forODE ODE solvers versus Step Size Simple Solvers 1.E+00 • Uses four derivative evaluations per step k  2k 2  2k 3  k 4 y i 1  y i  1 6 k1  hi 1 f ( xi , y i ) 1.E-02 xi 1  xi  hi 1 Error 1.E-07 Euler Huen Modified Euler Runge-Kutta 1.E-08 1.E-11 1.E-12 1.E-13 1.E-14 1.E-15 1.E-16 1.E-07 1.E-02 1.E-03 1.E-04 Euler (N = 10) Euler (N = 100) Euler N = 1,000) Huen (N =10) Huen (N = 100) Huen (N = 1,000) RK4 (H = 10) RK4 (N = 100) 1.E-05 1.E-06 1.E-07 1.E-08 1.E-09 1.E-10 1.E-11 1.E-12 1.E-13 1.E-14 0.8 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 Step size, h 8 Error Control 1.E-01 Error 1.E-04 1.E-10 1.E+00 0.6 with y  y0  1 at x  x0  0 1.E-09 Error Propagation in Solutions of ODEs 0.4 dy  e- x- y dx 1.E-06 7 0.2 1.E-03 1.E-05 h k k 2  hi 1 f  xi  i 1 , y i  1  2 2  h k 3  hi 1 f  xi  i 1 , y i  k 2  2 2  k 4  hi 1 f ( xi  hi 1 , y i  k 3 ) 0 1.E-01 1 • How do we choose h to maintain desired accuracy? • Want to obtain a result with some desired small global error • Can just repeat calculations with smaller h until two results are sufficiently close • Can use algorithms that estimate error and adjust step size during the calculation based on the error 9 10 x Runge-Kutta Error Control Dormand-Prince Method • Take step with two formulas of different orders using the same function evaluations in each formula • Use the difference between the two formulas as an estimate of the error • Use the error estimate to control the step size based on a user-input desired error • Some methods in final slides not shown • Uses fourth- and fifth-order expressions to find error estimate • Modifies step size by comparing error estimate with desired error • Used in MATLAB solver ode45 which is in next programming assignment • Algorithm equations shown on next two v1.11 of ode45 code at slides 11 http://users.powernet.co.uk/kienzle/octave/ matcompat/scripts/ode_v1.11/ode45.m 12 2 Systems of Ordinary Differential Equations Dormand-Prince Equations Dormand-Prince Method II 44 56 32 − + 45 15 9 44 − 168 + 160 = 45 36 4 = = 45 5 e is desired error per step hopt is optimum step size Toshinori Kimura, “On Dormand-Prince Method http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf Toshinori Kimura, “On Dormand-Prince Method http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf 13 14 Next Steps Basic Idea • We have shown some algorithms for solving the initial value problem, dy/dx = f(x,y) with y = y0 at x = x0 • Need to solve higher order equations • Can show than any nth-order equation can be expressed as a system of n firstorder equations • Then have to consider solving systems of first order equations • Any nth order ODE can be written as a system of n first-order ODEs • For a system of two or more ODEs with varying orders: – Break the kth ODE with order nk into nk firstorder ODEs; repeat for all equations – The result will be a system of first order ODEs • The total number of ODEs will be n1+ n2 + … = the sum of the orders of each original ODE 15 Higher order equations Higher Order Equations II • Look at spring-mass-damper equation as example: md2y/dt2 + Cdy/dt + ky = F – F is applied force which may be a known function of t, y and dy/dt (may have F = 0) – Define a velocity, v = dy/dt – dv/dt = d2y/dt2 – Our second order equation becomes mdv/dt + Cv + ky = F giving two ODEs dv C  F - v - ky  f v t , y, v  dt m dy  v  f y t , y, v  dt 16 dv  f v t , y, v  dt dy  f y t , y, v  dt General Notation dy  f t , y  dt dyk  f k t , y1  y N  dt • The two first order equations we found look like our initial value problem – This assumes that we have an initial condition for velocity, v • Will cover situations where this is not true (known as boundary value problems) later – We can apply all algorithms developed for single first-order equations to any system of simultaneous first-order equations 17 18 3 Systems of Ordinary Differential Equations Solving Simultaneous ODEs Runge-Kutta for ODE System • Apply same algorithms used for single ODEs – y(n) is vector of dependent variables at x = xn – k(1), k(2), k(3), and k(4), are vectors containing intermediate Runge-Kutta results – f is a vector containing the derivatives – k(1) = hf = hf(xn, y(n)) – k(2) = hf(xn + h/2, y(n) + k(1)/2) – k(3) = hf(xn + h/2, y(n) + k(2)/2) – k(4) = hf(xn + h, y(n) + k(3)) – y(n+1) = (k(1) + 2k(2) + 2k(3) + k(4))/6 – Must apply each part of each algorithm step to all equations in system before going on to next step – Key is having consistent x and y values in determination of fk(x,y) – All yk values in y must be available at the same x point when computing the fk – E.g., in Runge-Kutta we must evaluate k1 for all equations before finding k2 19 20 ODE System by RK4 ODE System by RK4 II • dy/dx = –y + z and dz/dx = y – z with y(0) = 1 and z(0) = -1 with h = .1 • Details of first step from y0 to y1 • k(1)y = h[-y + z] = 0.1[-1 + (-1)] = -.2 • k(1)z = h[y - z] = 0.1[1 - (-1)] = .2 • k(2)y = h[-(y+ k(1)y/2) + z + k(1)z/2] = 0.1[ -(1 + -0.2/2) + (-1 + .2/2)] = -.18 • k(2)z = h[(y+ k(1)y/2) –(z + k(1)z/2)] = 0.1[(1 + -0.2)/2 - (-1 + .2/2)] = .18 • k(3)y = h[-(y+ k(2)y/2) + z + k(2)z/2] = 0.1[ -(1 + -0.18/2) + (-1 + .18/2)] = -.182 • k(3)z = h[(y+ k(2)y/2) –(z + k(2)z/2)] = 0.1[(1 + -0.18)/2 - (-1 + .18/2)] = .182 • k(4)y = h[-(y+ k(3)y) + z + k(3)z] = 0.1[ -(1 + -0.182) + (-1 + .182)] = -.1636 • k(4)z = h[(y+ k(3)y) –(z + k(3)z)] = 0.1[ (1 + -0.182) - (-1 + .182)] = .1636 21 22 ODE System by RK4 III Numerical Software for ODEs • yi+1 = yi + (k(1)y + 2k(2)y + 2k(3)y + k(4)y )/6 = 1 + [ (–.2) + 2(–.18) + 2(–.182) + (–.1636)]/6 = .8187 (here i = 0) • zi+1 = zi + (k(1)z + 2k(2)z + 2k(3)z + k(4)z )/6 = –1 + [(.2) + 2(.18) + 2(.182) + (.1636)]/6 = –.8187 • Continue in this fashion until desired final x value is reached – Note all km computed before any km+1 • Usually written to solve a system of N equations, but will work for N = 1 • User has to code a subroutine or function to compute the f array – Input variables are x and y; f is output – Some codes have one dimensional parameter array to pass additional information from main program into the function that computes derivatives • No x dependence for f in this example 23 24 4 Systems of Ordinary Differential Equations Derivative Subroutine Example • Visual Basic sub code for system of ODEs at right is shown below dy1  - y1  y2 - y3e 2 x dx dy3 dy2  -2 y12  -3 y1 y2 dx dx Sub fsub( x as Double, y() as Double, _ f() as Double ) f(1) = -y(1) + Sqr(y(2)) - y(3)*Exp(2*x) f(2) = -2 * y(1)^2 f(3) = -3 * y(1) * y(2) End Sub 25 Derivative Function Example dy 1  - y1  y2 - y3e 2 x • Visual Basic dx function code for dy3 dy2  -2 y12  -3 y1 y2 system of ODEs dx dx Function ff( x As Double, y() As Double, _ N As Long) As Variant Dim fd() As Double; reDim fd(1 to N) fd(1) = -y(1) + Sqr(y(2)) - y(3)*Exp(2*x) fd(2) = -2 * y(1)^2 fd(3) = -3 * y(1) * y(2) : ff = fd End Function In main program use f = ff(x,y,N) where f is an array of same size 26 Derivative Subroutine Example Derivative Subroutine Example • MATLAB function for system of ODEs at right is shown below • C++ code for system of ODEs at right is shown below dy1  - y1  y2 - y3e 2 x dx dy3 dy2  -2 y12  -3 y1 y2 dx dx dy1  - y1  y2 - y3e 2 x dx dy3 dy2  -2 y12  -3 y1 y2 dx dx function f = test( x, y) f = zeros(3,1); f(1) =-y(1) +Sqrt(y(2)) - y(3)*Exp(2*x); f(2) = -2 * y(1)^2; f(3) = -3 * y(1) * y(2); end %semicolons prevent output 27 void fsub(double x, double y[], double f[]) { f[1] = -y[1] + sqrt(y[2]) - y[3]*exp(2*x); f[2] = -2 * y[1] * y[1]; f[3] = -3 * y[1] * y[2]; } 28 Derivative Subroutine Example Assignment Seven and Last • Fortran code for system of ODEs at right is shown below dy1  - y1  y2 - y3e 2 x dx dy3 dy2  -2 y12  -3 y1 y2 dx dx subroutine fsub( x, y, f ) real(KIND=8) x, y(:), f(:) f(1) = -y(1) + sqrt(y(2)) -y(3)*exp(2*x) f(2) = -2 * y(1)**2 f(3) = -3 * y(1) * y(2) end subroutine fsub 29 • Write routine in MATLAB and in VBA to solve general system of N ODEs with following form: dym/dt = fm(t,y), m = 1, N • Apply to N = 3 system with known exact solutions same as on previous slides dy1  - y1  y2 - y3e 2 x dx dy3 dy2  -2 y12  -3 y1 y2 dx dx y1  e - t y 2  e - 2t y 3  e - 3t 30 5 Systems of Ordinary Differential Equations Assignment Seven Algorithms • Code Huen algorithm in MATLAB • Modified Euler algorithm in VBA y m0 ,i 1  y m,i  hf m (t i , y i ) y m ,i y m,i 1  t i 1  t i  h  y m0 ,i 1  hf m (t i 1 , y 0m,i 1 ) 2 Runge-Kutta Error Control • Following slides not shown in class • Present other approach to error control in Runge-Kutta algorithms h f m (t i , y i ) 2  y m,i  hf m (t i  h 2 , y i 1 2 ) t i 1  t i  h – Older method: Compare results for two steps at step size h with one step of 2h y m,i 1 2  y m,i  y m,i 1 • Requires 12 derivative function evaluations over both steps • Use MATLAB function ode45 for solving same set of equations and compare work for same accuracy – Runge-Kutta-Fehlberg: early method to use algorithms of different order with same function evaluations 31 32 Runge-Kutta Error Control Runge-Kutta Error Control II • Control error by doing integration with h and 2h all along integration • y2h – yh = D is measure of truncation error • User specifies DD, the desired error – Integration with 2h step requires 3 additional function evaluations per 2 steps – Analyze local truncation error, which is O(h5), for both steps – Many ways to specify this, single value, relative values, relative to increments for y in one step • Since error scales as h5, we can adjust step size such that hnew = hold|DD/D|1/5 • Typically use safety factor to avoid making hnew too large y ( x  2h)  yh  Ah5  Bh6   y ( x  2h)  y2 h  A2h   B2h    0  y2 h - yh  25 - 1 Ah5  O(h6 ) 5  6  33 34 Runge-Kutta-Fehlberg Runge-Kutta-Fehlberg II • Uses two equations to compute yn+1, one has O(h5), the other O(h6) error • Requires six derivative evaluations per step (same evaluations used for both equations) • The error estimate can be used for step size control based on an overall 5th order error • Cask-Karp version and Runge-KuttaVerner use same idea 35 • • • • • • • • See page 663 in Rao text for algorithm Typical formula components below yn+1= yn + (16k1/135 + 6656k2/12825 … k3 = hf(xn + 3h/8, yn + 3k1/32 + 9k2/32) Error = k1/360 - 128k3/4275 … hnew = hold || holdEDesired/Error ||1/4 EDesired is set by user RKF45 code by Watts and Shampine 36 6 Systems of Ordinary Differential Equations Increments, k, for RKF Final RKF Step Results k1  hf x j , y j  h k   k 2  hf  x j  , y j  1  4 4  3h 3k  9k 2   k3  hf  x j  , y j  1  8 32   1932k1 - 7200k 2  7296k3  12h  k 4  hf  x j  , yj   13 2197   3680k3 845k 4  439k1  k5  hf  x j  h, y j  - 8k 2   216 513 4104   3544k3 1859k 4 11k5  h 8k  k6  hf  x j  , y j - 1  2k 2   2 27 2565 4104 40   37 • Proposed value of yi+1 is yi 1  yi  16k1 6656k3 28561k4 9k5 2k6    135 12825 56430 60 55 • Estimated error is e i 1  • • • • 16k1 6656k3 28561k4 9k5 2k6    135 12825 56430 60 55 If emin< ei+1 < emax accept yi+1; keep h If ei+1 < emin accept yi+1 and get new h If ei+1 > emax get new h and redo step hnew = hold(ε*hold/|ei+1|)¼ ε* = desired error 38 7