Document Sample

Introduction to Numerical Solutions of Ordinary Differential Equations Larry Caretto Mechanical Engineering 501AB Seminar in Engineering Analysis October 29, 2003 Outline • Review last class and homework • Midterm Exam November 4 covers material up to and including last week’s lecture • Overivew of numerical solutions – Initial value problems in first-order equations – Systems of first order equations and initial value problems in higher order equations – Boundary value problems – Stiff systems and eigenvalues 2 Review Last Class • Systems of equations – Combine to higher order equation – Matrix approaches – Reduction of order • Laplace transforms for ODEs – Transform ODE from y(t) to Y(s) – Solve for Y(s) – Rearrange solution – Transform back to y(t) 3 Simultaneous Solution • Differentiate first (y1) equation to obtain derivative(s) of y2 that are present in second equation • Solve first equation for y2 • Substitute equations for y2 and its deriv- atives into the result of the first step • Result is a higher order equation that can be solved (if possible) by usual methods 4 Matrix Differential Equations II dy • Matrix components in Ay r dt y1 dy1 y y1 dt a11 a12 a13 a1n y dy2 a a23 a2 n 2 21 a22 2 dt y3 dy d y3 dy3 a31 a32 a33 a3n y dt A dt dt yn dyn an1 an 2 an 3 ann yn dt r r1 r2 r3 rn T 5 dy Solving Ay r dt • Assume that A(n x n) has n linearly independent eigenvectors • Eigenvectors are columns of matrix, X • Define a new vector s = X-1y (y = Xs) • Substitute y = Xs into the matrix differential equation • Obtain equation for s using L = X-1AX that gives independent equations 6 ds Terms in Solution of Ls p dt l1t e p1 0 0 0 C1 l1 C 0 e l2 t 0 0 p2 2 l 0 e l3t C3 E(t ) 0 0 C 1 p 2 L p 3l 3 0 ln t e 0 0 Cn p nl • With these definitions, i s = C elit + p /l n i i i becomes s = EC + L-1p (E(0) = I) 7 Apply Initial Conditions on y(0) • Get y = Xs = XEC + XL-1p • At t = 0, E = I, and y = y0 = initial y components, giving y0 = XIC + XL-1p • Premultiply by X-1 to obtain X-1y0 = X-1XC + X-1XL-1p = C + L-1p • Constant vector, C = X-1y0 - L-1p • Result: y = XE [X-1y0 - L-1p] + XL-1p • Homogenous(p = 0): y = XEX-1y0 8 Reduction of Order • An nth order equation can be written as a system of n first order equations • We can write a general nonlinear nth order equation as shown below dny dy d n1 y f x, y, ,, n1 dx n dx dx • Redefine y as z0 and define the following derivatives dy dz0 d 2 y d 2 z0 dz1 z1 2 2 z2 dx dx dx dx dx 9 Reduction of Order II • Continue this definition up to zn-1 d n y d n z0 d n 1 z1 d n 2 z 2 dzn 2 dzn 1 n n n 1 n2 2 dx dx dx dx dx dx • z0 to zn-1 are n variables that satisfy simultaneous, first-order ODEs dy n d n1 y dzn1 dx n f x, z0 , z1 ,, zn 1 dy n f x, y, ,, n1 dx dx dx dzk zk 1 k 0,1,, n 2 dx 10 Reduction of Order III • The main application of reduction of order is to numerical methods • Numerical solutions of ODEs develop methods to solve a system of first order equations • Higher order equations are solved by converting them to a system of first order equations 11 Simple Laplace Transforms f(t) F(s) f(t) F(s) ctn cn!/sn+1 ceatsin wt cw ctx cG(x+1)/sx+1 ( s a) 2 w 2 ceat c/(s – a) ceatcos wt c( s a) csin wt cw/(s2 + w2) ( s a) 2 w 2 ccos wt cs/(s2 + w2) Additional transforms csinh wt cw/(s2 - w2) in Table 5.9, pp 297- 299 of Kreyszig ccosh wt cs/(s2 - w2) 12 Transforms of Derivatives • The main formulas for differential equations are the Laplace transforms of derivatives • These have transform of desired function, L [f(t)] = F(s) and the initial conditions, f(0), f’(0), etc. [ f ' (t )] sF (s) f (0) [ f ' ' (t )] s 2 F ( s) sf (0) f ' (0) [ f ( n ) (t )] s n F ( s) s n 1 f (0) s n 2 f ' (0) sf ( n 2) (0) f ( n 1) (0) 13 Solving Differential Equations • Transform all terms in the differential equation to get an algebraic equation – For a differential equation in y(t) we get the transforms Y(s) = L [y(t)] – Similar notation for other transformed functions in the equation R(s) = L [r(t)] • Solve the algebraic equation for Y(s) • Obtain the inverse transform for Y(s) from tables to get y(t) 14 Partial Fractions • Method to convert fraction with several factors in denominator into sum of individual factors (in denominator) • Necessary to modify complicated solution for Y(s) into simpler functions whose transforms are known • Special rules for repeated factors and for complex conjugates 15 Shifting Theorems • Sometimes required for transform table – First theorem is defined for function f(t) with transform F(s) – F(s - a) has inverse of e-atf(t) – Example Y(s) = [2(s+1)+4]/[(s + 1)2 + 1] – For F(s) = s/(s2 + w2); f(t) = cos wt and f(t) = cos wt for F(s) = w/(s2 + w2) – Here Y(s+1) = 2(s+1)+4]/[(s + 1)2 + 1] + 4/[(s + 1)2 + 1] so y(t) = e-t times Y(s) inverse – y(t) = e-t[2 cos t + 4 sin t] 16 Shifting Theorems II • Second shifting theorem – Applies to e-asF(s), where F(s) is known transform of a function f(t) – Inverse transform is f(t – a) u(t – a) where u(t – a) is the unit step function – If we have e-ass/(s2 + w2), we recognize s/(s2 + w2) as F(s) for f(t) = cos wt – Thus e-ass/(s2 + w2) is the Laplace transform for cos[wt(t – a)] u(t – a) 17 Numerical Analysis Problems • Numerical solution of algebraic equations and eigenvalue problems • Solution of one or more nonlinear algebraic equations f(x) = 0 • Linear and nonlinear optimization • Constructing interpolating polynomials • Numerical quadrature • Numerical differentiation • Numerical differential equations 18 Interpolation • Start with N data pairs xi, yi • Find a function (polynomial) that can be used for interpolation • Basic rule: the interpolation polynomial must fit all points exactly • Denote the polynomial as p(x) • The basic rule is that p(xi) = yi • Many different forms 19 Newton Polynomials • p(x) = a0 + a1(x – x0) + a2(x – x0)(x – x1) + a3(x – x0)(x – x1)(x – x2) + … + an-1(x – x0)(x – x1)(x – x2) … (x – xn-2) • Terms with factors of x – xi are zero when x = xi – Use this and rule that p(xi) = yi to find ai • a0 = y0, a1 = (y1 – y0) / (x1 – x0) • y2 = a0 + a1(x2 – x0) + a2(x2 – x0)(x2 – x1) – Solve for a2 using results for a0 and a1 20 Newton Polynomials II • y2 = a0 + a1(x2 – x0) + a2(x2 – x0)(x2 – x1) y1 y0 y2 y0 ( x2 x0 ) y2 a0 a1 ( x2 x0 ) x1 x0 a2 ( x2 x0 )(x2 x1 ) ( x2 x0 )(x2 x1 ) • y2 = a0 + a1(x2 – x0) + a2(x2 – x0)(x2 – x1) • Data determine coefficients • Develop scheme known as divided difference table to compute ak 21 Divided Difference Table x0 y0 a0 y1 y0 F0 a1 x1 x0 F1 F0 a2 x1 y1 S0 x2 x0 y2 y1 S1 S 0 F1 T0 x2 x1 x3 x0 F2 F1 x2 y2 S1 a3 x3 x1 y3 y 2 F2 x3 x2 x3 y3 22 Divided Difference Example 0 0 a0 10 0 F0 1 a1 10 0 3 1 10 10 S0 .1 a2 20 0 40 10 .2 .15 1 F1 3 T0 20 10 30 0 600 20 40 63 a3 S1 .15 30 10 100 40 F2 6 30 20 30 100 23 Divided Difference Example II • Divided difference table gives a0 = 0, a1 = 1, a2 = .1, and a3 = 1/600 • Polynomial p(x) = a0 + a1(x – x0) + a2(x – x0)(x – x1) + a3(x – x0)(x – x1)(x – x2) = 0 + 1(x – 0) + 0.1(x – 0)(x – 10) + (1/600)(x – 0)(x – 10)(x – 20) = x + 0.1x(x – 10) + (1/600)x(x – 10)(x – 20) • Check p(30) = 30 + .1(30)(20) + (1/600) (30)(20)(10) = 30 + 60 + 10 = 100 (correct) 24 Constant Step Size • Divided differences work for equal or unequal step size in x • If Dx = h is a constant we have simpler results – Fk = Dyk/h = (yk+1 – yk)/h – Sk = D2yk/h2 = (yk+2 – 2yk-1 + yk)/h2 – Tk = D3yk/h3 = (yk+3 – 3yk+2 + 3yk+1 – yk)/h3 – Dnyk is called the nth forward difference – Can also define backwards and central differences 25 Interpolation Approaches • When we have N data points how do we interpolate among them? – Order N-1 polynomial not good choice – Use piecewise polynomials of lower order (linear or quadratic) – Can match first and or higher derivatives where piecewise polynomials join – Cubic splines are piecewise cubic polynomials that match first and second derivatives (as well as values) 26 Cubic Spline Interpolation 0.7 0.6 0.5 y values Known f' 0.4 Natural No Knot 0.3 Data 0.2 0.1 0 0 1 2 3 4 5 6 x values 27 Newton Interpolating Polynomial 5 4 3 Polynomial Y Values Data 2 1 0 -1 0 1 2 3 4 5 6 X Values 28 Polynomial Applications • Data interpolation • Approximation functions in numerical quadrature and solution of ODEs • Basis functions for finite element methods • Can obtain equations for numerical differentiation • Statistical curve fitting (not discussed here) usually used in practice 29 Derivative Expressions • Obtain from differentiating interpolation polynomials or from Taylor series • Series expansion for f(x) about x = a df 1 d2 f 1 d3 f f ( x) f (a ) ( x a) ( x - a) 2 ( x - a)3 .... dx x a 2! dx2 xa 3! dx3 xa • Note: d0f/dx0 = f 1 dn f f ( x) ( x - a) n and 0! = 1 n 0 n ! dx n x a • What is error from truncating series? 30 Truncation Error • If we truncate series after m terms m 1 dn f 1 dn f f ( x) ( x - a) n ( x - a) n n 0 n ! dx n x a n m 1 n ! dx n xa Terms used Truncation error, em • Can write truncation error as single term at unknown location (derivation based on the theorem of the mean) 1 dn f 1 d m1 f em ( x - a) n ( x - a) m1 n m 1 n ! dx n xa (m 1) ! dx m1 x 31 Derivative Expressions • Look at finite-difference grid with equal spacing: h = Dx so xi = x0 + ih • Taylor series about x = xi gives f(xi + kh) = f[x0 + (i+k)h] = fi+k in terms of f(xi) = fi df 1 d2 f 1 d3 f f ( xi kh) f ( xi ) kh (kh) 2 (kh)3 ..... dx x xi 2! dx 2 x xi 3! dx3 x xi • Compact derivative notation df d2 f dn f fi ' fi 2 '' ... f i n n dx x xi dx x xi dx x xi 32 Derivative Expressions II • Combine all definitions for compact series notation df 1 d2 f 1 d3 f f ( xi kh) f ( xi ) kh (kh) 2 (kh)3 ..... dx x xi 2! dx 2 x xi 3! dx3 x xi f i '' (kh) 2 f i ''' (kh)3 fik f i f i ' kh ..... 2! 3! • Use this formula to get expansions for various grid locations about x = xi and use results to get derivative expressions 33 Derivative Expressions III • Apply general f i '' (kh) 2 f i ''' (kh)3 fik f i f i ' kh .. equation for k 2! 3! = 1 and k = –1 f i '' h 2 f i ''' h 3 f i 1 f i f i ' h ..... 2! 3! f i '' h 2 f i ''' h 3 f i 1 f i f i ' h ..... 2! 3! f i 1 f i f i '' h f i ''' h f i 1 f i fi ' ..... Ah Forward h 2! 3! h f i f i 1 f i '' h f i ''' h f i f i 1 fi ' ..... Ah h 2! 3! h Backward 34 Derivative Expressions IV • Subtract fi+1 and fi-1 expressions f i '' h 2 f i ''' h 3 f i 1 f i f i ' h ..... 2! 3! f i '' h 2 f i ''' h 3 f i 1 f i f i ' h ..... 2! 3! 2 f i '' ' h 3 2 f i ' '' h 5 f i 1 f i 1 2 f i ' h ..... 3! 5! f i 1 f i 1 f i ''' h 2 f i ''''' h 4 f i 1 f i 1 fi ' ..... Ah 2 h 3! 5! 2h • Result called central difference expression 35 Order of the Error • Forward and backward derivative have error term that is proportional to h • Central difference error is proportional to h2 • Error proportional to hn called nth order • Reducing step size by a factor of a reduces nth order error by an n h2 e 2 e1 h 1 36 Order of the Error Notation • Write the error term for nth error term as O(hn) – Big oh notation, O, denotes order – Recognizes that factor multiplying hn may change slightly with h First order forward First order backward f i 1 f i f i f i 1 fi ' O ( h) fi ' O ( h) h h Second order central f i 1 f i 1 fi ' O( h 2 ) 2h 37 Higher Order Derivatives • Add fi+1 and fi-1 expressions f i '' h 2 f i ''' h 3 f i 1 f i f i ' h ..... 2! 3! f i '' h 2 f i ''' h 3 f i 1 f i f i ' h ..... 2! 3! f i ' ' h 2 2 f i '' '' h 4 2 f i '' '''' h 6 f i 1 f i 1 2 f i 2 ..... 2! 4! 6! f i 1 f i 1 2 f i f i ''' h 2 f i ''''' h 4 f i 1 f i 1 2 f i fi '' h2 3! 5! ..... h2 O h2 • f’’ is second-order, central difference 38 Higher Order Directional • We can get higher order derivative expressions at the expense of more computations • Get second order forward and backward derivative expressions from fi+2 and fi-2, respectively • Combine with previous expressions for fi+1 and fi-1 to eliminate first order error term 39 Specific Taylor Series • General f i '' (kh) 2 f i ''' (kh)3 f i k f i f i ' kh ..... equation 2! 3! f i '' h 2 f i '' ' h 3 fi 2 fi 2 fi 'h 4 8 ..... • k=2 2! 3! f i '' h 2 f i ''' h 3 f i 2 f i 2 f i 'h 4 8 ..... • k = -2 2! 3! f i '' h 2 f i ' '' h 3 • k=3 f i 3 f i 3 f i ' h 9 27 ..... 2! 3! • k = -3 f i '' h 2 f i '' ' h 3 f i 3 f i 3 f i ' h 9 27 ..... 2! 3! 40 Second Order Forward • Subtract 4fi+1 from fi+2 to eliminate h2 term '' h 2 ''' h 3 f i 2 4 f i 1 f i 2 f i h 4 f i ' 8 fi ... 2 6 '' h 2 ''' h 3 4 f i f i h f i ' fi .. 2 6 h 3Second f i 2 4 f i 1 3 f i 2hf i ' 4 f i ''' ... order 6 error f i 2 4 f i 1 3 f i '' ' h 2 fi ' fi ... 2h 3 41 Second Order Backwards • Add 4fi-1 to –fi-2 to eliminate h2 term '' h 2 ''' h 3 f i 2 4 f i 1 f i 2 f i h 4 f i ' 8 fi ... 2 6 '' h 2 ''' h 3 4 f i f i h f i ' fi .. 2 6 h Second 3 f i 2 4 f i 1 3 f i 2hfi 4 f i ' ... ' '' 6 order error f i 2 4 f i 1 3 f i ''' h 2 fi ' fi ... 2h 3 42 Other Derivative Expressions • Can continue in this fashion – Write Taylor series for fi+1, fi-1, fi+2, fi-2, fi+3, fi-3, etc. – Create linear combinations with factors that eliminate desired terms – Eliminate fi term to obtain central difference – Keep only terms in fk with k i for forward difference expressions – Keep only terms in fk with k i for forward difference expressions – See results on page 271 of Hoffman 43 Order of Error Examples • Table 2-1 in notes shows error in first derivative for ex around x = 1 – Using first- and second-order forward and second-order central differences – Step h = 0.4, 0.2, and 0.1 – Error ratio for doubling step size • 4.01 to 4.02 for central differences • 2.07 to 2.15 for first-order forward differences • 4.32 to 4.69 for second-order forward log e 2 e log(e ) log(e ) n 1 2 1 log 2 log( h2 ) log( h1 ) h h 44 1 Roundoff Error • Possible in derivative expressions from subtracting close differences • Example f(x) = ex: f’(x) (ex+h – ex-h)/(2h) and error at x = 1 is (e1+h – e1-h)/(2h) - e 3.004166 2.722815 E 2.718282 4.5 x10 3 2(0.1) Second order error 2.7185536702 2.7180100139 E 2.7182818284 59 4.5 x10 9 2(0.0001 ) 2.7182821002 8724 2.7182815566 0388 E 2.718281828 5.9 x10 9 2(0.0000001 ) 45 Figure 2-1. Effect of Step Size on Error 1.E+01 1.E+00 1.E-01 1.E-02 1.E-03 1.E-04 Error 1.E-05 1.E-06 1.E-07 1.E-08 1.E-09 1.E-10 1.E-11 1.E-17 1.E-15 1.E-13 1.E-11 1.E-09 1.E-07 1.E-05 1.E-03 1.E-01 Step Size 46 Numerical ODE Solutions • The initial value problem • Euler method as a prototype for the general algorithm • Local and global errors • More accurate methods • Step-size control for error control • Applications to systems of equations • Higher-order equations 47 The Initial Value Problem • dy/dx = f(x,y) (known) with y(x0) = y0 • Basic numerical approach – Use a finite difference grid: xi+1 – xi = h – Replace derivative by finite-difference approximation: dy/dx (yi+1 – yi) / (xi+1 – xi) = (yi+1 – yi) / h – 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) / h = favg – Repeatedly compute yi+1 = yi + h favg 48 Notation • xi is the value of the independent variable at point i on the grid – Determined from the user-selected value of step size (or the series of hi values) – Can always specify exactly the independent variable’s value, xi • yi is the value of the numerical solution at the point where x = xi • fi is derivative value found from xi and the numerical value, yi. I.e., fi = f(xi, yi). 49 More Notation • y(xi) is the exact value of y at x = xi – Usually not known but notation is used in error analysis of algorithms • f(xi,y(xi)) is the exact value of the derivative at x = xi • e1 = y(x1) – y1 is the local truncation error – This is error for one step of algorithm starting from known initial condition 50 Local versus Global Error • At the initial condition we know the solution, y, exactly • First step introduces some error • Remaining steps have single step error plus previous accumulated error • Ej = y(xj) – yj is global truncation error – Difference between numerical and exact solution after several steps – This is the error we want to control 51 Euler’s Method • Simplest algorithm, example used for error analysis, not for practical use • Define favg = f(xi, yi) • Euler’s method algorithm is yi+1 = yi + hifi = yi + hi f(xi, yi) • Example dy/dx = x + y, y = 0 at x = 0 • Choose h = 0.1 • We have x0 = 0, y0 = 0, f0 = x0 + y0 = 0 x1 = x0 + h = 0.1, y1 = y0 + hf0 = 0 52 Euler Example Continued • Next step is from x1 = 0.1 to x2 = 0.2 • f1 = x1 + y1 = 0.1 + 0 = 0.1 • y2 = y1 + hf1 = 0 + (.1)(.1) = .01 • For dy/dx = x + y, we know the exact solution is (x0 + y0 + 1)ex-x0 – x – 1 • For x0 = y0 = 0, y = ex – x – 1 • Look at application of Euler algorithm for a few steps and compute the error 53 Euler Example xi yi fi y(xi) E(xi) 0 0 0 0 0 .1 0 .1 .0052 .0052 .2 .01 .21 .0214 .0114 .3 .031 .331 .04986 .01886 .4 .0641 .4641 .091825 .027725 54 Euler Example Plot 0.8 0.7 0.6 0.5 0.4 y Exact Numerical 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 x 55 Error Propagation • Behavior of Euler algorithm is typical of all algorithms for numerical solutions • Error grows at each step • We usually do not know this global error, but we would like to control it • Look at local error for Euler algorithm • Then discuss general relationship between local and global error 56 Taylor Series to Get Error • Expand y(x) in Taylor series about x = a dy 1 d2y 1 d3y y ( x) y (a ) ( x a) ( x - a) 2 ( x - a) 3 .... dx x a 2! dx 2 x a 3! dx 3 xa • Look at one step from known initial condition, a = x0, to x0 + h so x – a = h dy 1 d2y 2 1 d3y 3 y ( x 0 h) y ( x 0 ) h 2 h 3 h .... dx 0 2! dx 0 3! dx 0 • In ODE notation, dy/dx|0 = f(x0, y(x0)) 57 Local Euler Error • Result of Taylor series on last chart 1 d2y 2 1 d3y 3 y( x0 h) y( x0 ) hf ( x0 , y ( x0 )) 2 h 3 h .... 2! dx 0 3! dx 0 Euler Algorithm Truncation Error • This is only the Euler algorithm for the first step when we know f(x0,y(x0)) • This gives the local truncation error • Local truncation error for Euler algorithm is second order 58 Global Error • We will show that a local error of order n, has a global error of order n-1 • To show this consider the global error at x = x0 + kh after k algorithm steps – Is approximately k times the local error – If local error is O(hn), approximate global error after k steps is k O(hn) kAhn – A new step size, h/r, takes kr steps to get to the same x value 59 Global Error Concluded • Compare error for same x = kh with step sizes h and h/r • Ex=kh(h) kAhn Ex kh h r kr h n r 1 • Ex=kh(h/r) = krA(h/r) n k h n n 1 Ex kh (h) r • When we reduce the step size by a factor of 1/r we reduce the error by a factor of 1/rn-1; this is the behavior of an algorithm whose error is order n-1 60 Euler Local and Global Error • Previously showed Euler algorithm to have second order local error • Should have first order global error • Results for previous Euler example at x = 1 with different step sizes Step size First step Final error h = 0.1 5.17x10-3 1.25 x10-1 h = 0.01 5.02 x10-5 1.35 x10-2 h = 0.001 5.00 x10-7 1.36 x10-3 61 Better Algorithms • Seek high accuracy with low computational work • Could improve Euler accuracy by cutting step size, but this is not efficient • Use other algorithms that have higher order errors • Runge-Kutta methods typically used – This is a class of methods that use several function evaluation methods per step 62 Second-order Runge Kutta • Huen’s method yi01 yi hi 1 f ( xi , yi ) xi 1 xi hi 1 yi 1 yi hi 1 2 f ( xi , yi ) f ( xi 1 , yi01 ) yi yi01 hi 1 f ( xi 1 , yi01 ) 2 • Modified Euler method hi 1 hi 1 yi 1 yi f ( xi , y i ) xi 1 xi 2 2 2 2 y i 1 y i hi 1 f ( xi 1 , y i 1 ) 2 2 63 Fourth-order Runge Kutta • Uses four derivative evaluations per step k1 2 k 2 2 k 3 k 4 y i 1 y i xi 1 xi hi 1 6 k1 hi 1 f ( xi , y i ) k 2 hi 1 f xi i 1 , y i 1 h k 2 2 x hi 1 , y k 2 k 3 hi 1 f i 2 i 2 k 4 hi 1 f ( xi hi 1 , y i k 3 ) 64 Comparison of Methods • Look at Euler, Heun, Modified Euler and fourth-order Runge-Kutta • Solve dy/dx = e-y-x with y(0) = 1 • Compare numerical values to exact solution y = ln( ey0 + e-x0 – e-x) • Look at errors in the methods at x = 1 as a function of step size • Compare error propagation (increase in error as x increases) 65 Error versus Step Size for Simple ODE Solvers 1.E+00 dy e x y with y y0 1 at x x0 0 1.E-01 1.E-02 1.E-03 dx 1.E-04 1.E-05 1.E-06 1.E-07 Euler Error 1.E-08 Huen 1.E-09 Modified Euler Runge-Kutta 1.E-10 1.E-11 1.E-12 1.E-13 1.E-14 1.E-15 1.E-16 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 Step size, h 66 Error Propagation in Solutions of ODEs 1.E+00 1.E-01 1.E-02 1.E-03 1.E-04 Euler (N = 10) 1.E-05 Euler (N = 100) 1.E-06 Euler N = 1,000) Error Huen (N =10) 1.E-07 Huen (N = 100) 1.E-08 Huen (N = 1,000) 1.E-09 RK4 (H = 10) RK4 (N = 100) 1.E-10 1.E-11 1.E-12 1.E-13 1.E-14 0 0.2 0.4 0.6 0.8 1 67 x Some Basic Concepts • A numerical method is convergent with the solution of the ODE if the numerical solution approaches the actual solution as h 0 (with increase in numerical precision at smaller h) – Mainly a theoretical concept; algorithms in use are convergent • Stability refers to the ability of a numerical algorithm to damp any errors introduced during the solution 68 More on Stability • Finite difference equations in numerical algorithms, when iterated, may be numerically increase without bound • Stability usually is obtained by keeping step size h small, sometimes smaller than the h required for accuracy • For most ODEs stability is not a problem, but it is for stiff systems of ODEs and for partial differential equations 69 Euler Solution of Decaying Exponential 5 4.5 4 y(0) = 1 y(0) = 2 3.5 Euler algorithm errors move y(0) = 3 numerical solution of dy/dx = -y 3 y(0) = 4 to solution for anohter initial condition. Error bounded by y(0) = 5 y 2.5 decreasing exponential Euler 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 x 70 Euler Solution for Increasing Exponential 100 y(0) = 1 90 y(0) = 2 80 Euler algorithm errors move y(0) = 3 numerical solution of dy/dx = -y 70 y(0) = 4 to solution for anohter initial 60 y(0) = 5 condition. Error becomes large for increasing exponential Euler y 50 40 30 20 10 0 0 0.5 1 1.5 2 2.5 3 x 71 Error Control • How do we choose h? • 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 72 Runge-Kutta-Fehlberg • 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 different coefficients 73 Runge-Kutta-Fehlberg II • See page 377 in Hoffman 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|EDesired/Error|1/5 • EDesired is set by user • RKF45 code by Watts and Shampine 74 Other Error Controls • Do fourth-order Runge-Kutta two times, with step sizes different by a factor of two and compare results • Runge-Kutta-Verner is similar to Runge- Kutta-Fehlberg, but uses higher order expressions 75 Systems of Equations • Any problem with one or more ODEs of any order can be reduced to a system for first order ODEs • Last week we reduced a system of two second order equations to a system of four first order equations • In this process the sum of the orders is constant • Example from last week 76 Reduction of Order Example d 2 y1 k1 k2 k2 d 2 y2 k2 k3 k 2 2 y1 y2 0 2 y1 y2 0 dt m1 m1 dt m2 m2 • Define variables y3 = dy1/dt and y4 = dy2/dt • Then dy3/dt = d2y1/dt2 and dy4/dt = d2y2/dt2 • Have four simultaneous first-order ODEs dy1 dy2 y3 f1 ( x, y ) y 4 f 2 ( x, y ) dt dt dy3 k1 k 2 k2 y1 y 2 f 3 ( x, y ) dt m1 m1 dy4 k 2 k3 k 2 y1 y 2 f 4 ( x, y ) dt m2 m2 77 Solving Simultaneous ODEs • Apply same algorithms used for single ODEs – Must apply each step and substep to all equations in system – Key is consistent determination of fi(x,y) – All yi values in y must be available at the same x point when computing the fi – E.g., in Runge-Kutta we must evaluate k1 for all equations before finding k2 78 Runge-Kutta for ODE System – y(i) is vector of dependent variables at x = xi – 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(xi, y(i)) – k(2) = hf(xi + h/2, y(i) + k(1)/2) – k(3) = hf(xi + h/2, y(i) + k(2)/2) – k(4) = hf(xi + h, y(i) + k(3)) – y(i+1) = (k(1) + 2k(2) + 2k(3) + k(4))/6 79 ODE System by RK4 • dy/dx = -y + z and dz/dx = y – z with y(0) = 1 and z(0) = -1 with h = .1 • 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 80 ODE System by RK4 II • 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 81 ODE System by RK4 III • 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 • 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 • No x dependence for f in this example 82 Numerical Software for ODEs • 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 83 Derivative Subroutine Example • Visual Basic code dy1 y1 y2 y3e 2 x for system of dx ODEs at right is dy2 dy3 2 y12 3 y1 y2 shown below 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 84 Derivative Subroutine Example • Fortran code for dy1 y1 y2 y3e 2 x system of ODEs dx at right is shown dy2 dy3 2 y12 3 y1 y2 below 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 85 Derivative Subroutine Example • C++ code for dy1 y1 y2 y3e 2 x system of ODEs dx at right is shown dy2 dy3 2 y12 3 y1 y2 below dx dx 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]; } 86

DOCUMENT INFO

Shared By:

Categories:

Tags:

Stats:

views: | 6 |

posted: | 3/13/2012 |

language: | English |

pages: | 86 |

OTHER DOCS BY ewghwehws

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.