COMPUTATIONAL METHODS MCEN 3030
Popular in Course
Miss Cloyd Cronin
verified elite notetaker
Popular in Mechanical Engineering
This 29 page Class Notes was uploaded by Deron Schmitt on Thursday October 29, 2015. The Class Notes belongs to MCEN 3030 at University of Colorado at Boulder taught by Staff in Fall. Since its upload, it has received 19 views. For similar materials see /class/232015/mcen-3030-university-of-colorado-at-boulder in Mechanical Engineering at University of Colorado at Boulder.
Reviews for COMPUTATIONAL METHODS
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/29/15
Chapter 6 Numerical Solution of Ordinary Differential Equations 61 Introduction Consider for the moment the single nonlinear rst order ordinary differential equation ODE dyi E 7 ft7y7 with initial value y0 yoi We will always assume the existence and uniqueness of the solution and also that ft y has continuous partial derivatives with respect to t and y of as high order as necessary The aim of all numerical methods for solution of this differential equation is to obtain the solution at time twirl tn h given the solution for 0 g t g tut We begin by considering the so called Taylor series methods Let s expand the solution at tn1 about the solution at tn 3 yn1ynhyy y i n From the differential equation we have y flttn7yn To evaluate the higher order derivatives we will use the chain rule dy df 0f 0fdy 7 y 7 dt dt 0t0ydt ffw 0 0 ym EltftffygtfEltftffyftt2fftyftfyff5f2fyy It is clear that the number of terms increases rapidly and the method is not very practical for higher than second order accuracy The method based on the rst two terms in the expansion is called the forward or explicit Euler method yn1 yn hftnynA 61 One simply starts from the initial condition yo and marches forward using the formula to obtain y1y2iw We will study the properties of this method extensively as it is a very simple method to analyze From the Taylor series it is apparent that the forward Euler method is second order accurate for one time step However as with the quadrature formulas the global error for advancing from the initial condition to the nal time tf is only rst order accuratei 39 40 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Among the more accurate methods that we will discuss are the Runge Kutta formulas With Runge Kutta methods the solution at time step tn is obtained in terms of yn ftn yn and ft y evaluated at intermediate steps between tn and n1 not including tn1i The higher accuracy is achieved because more information about f is provided due to intermediate evaluations of f This is in contrast to the Taylor series method where we provided higher derivatives of f at tni Higher accuracy can also be obtained by providing information at times t lt tni That is the formulas involve yn71yn2 i H and fn71 1272 i i n These methods are called multistep methods We will also distinguish between explicit and implicit methods The preceding methods were all expliciti The formulas that involve fty evaluated at tn belong to the class of implicit methodsi Since f may be a nonlinear function of y to obtain the solution at each time step implicit methods usually require the solution of nonlinear equations Although the computational cost per time step is high implicit methods offer the advantage of numerical stability which will be discussed lateri 62 Analytical Concepts Before embarking on a deeper study of numerical solutions of ordinary differential equations it is useful to classify ODE s into categories rst order vsi higher order single equation vsi system of equations autonomous vsi nonautonomous linear vsi nonlinear homogeneous vsi inhomogeneous initial value vsi boundary value We discuss the relationships among these categories from the analytical as well as the computational points of view Almost any ODE of order higher than rst can be reduced to a system of rst order ODE si If we have the nth order equation yo f Lyyy ymvyW D and let yiyi 1gt i12iun then we can write y yi1 yfty1y2iuyn i12uin71i In addition if the system is nonautonomous take yo t so that ye 1 Then we formally have the autonomous system y giyj i012iun j012iun 62 where goyj 1 Thus a higher order single equation is equivalent to a system of rst order equations as long as the highest derivative can be isolated Therefore our initial emphasis will be on rst order ODE7si In general the same computational methods can be used to solve linear and nonlinear problems However usually nonlinear problems require considerably more computational cost In many cases a nonlinear equation is solved iteratively by using a linearization of the equation Let us linearize the righthandside of the system of nonlinear equations 62 about some known solution 0239 mm 9239 yi 0939 yjiy m y Bi l Ai j l uw 63 NUMERICAL CONCEPTS 41 where repeated indeces imply summationi Then the linearized system is y Bi 142739ij where Aij is an n 1 gtlt n 1 constant matrix and Bi is a constant vector of length n 1 Now this linear system can be solved and using this solution we can subsequently linearize about this new solution thus calculating new values of Aij and Bi This procedure can be repeated in an iterative manner until the nonlinear system is satis ed within some prescribed accuracy The distinction between homogeneous and inhomogeneous problems is not important from the numer ical analysis point of View except in boundary value problems It turns out that inhomogeneous terms while important in producing the correct solution are relatively unimportant in numerical stability and accuracy analyses In any case by using the transformation 1139 yj Ajilek the nonhomogeneous differential equation can be converted to the homogeneous equation Aijzji Thus without loss of generality we will consider the linear homogeneous system yE Aijyr 6 3 Assuming that Aij has n1 distinct eigenvalues A i 0 12i i i n letting Aij be the diagonal matrix composed of such eigenvalues and 33 the matrix of eigenvectors so that A SikAklslg1 or A SgAlelj then if we let 2239 321ka we have that 22 3121 SilAkzyz SilAlelj Sf ym A232 The above system is an uncoupled set of equations 225xizi i012 ini So except for the introduction of stz ness to be discussed later by a higher order equation or by a system of equations the solution methods and their analyses are similar to those used on a single rstorder equationi Note that the solutions of the above linear system are ziekltlgttci i012iun and in terms of our original dependent variables yz39 239ka Sik AktCk7 or y OeW co SileM cl Sinekwcni The distinction between initial and boundary value problems is important However since some meth ods used to solve boundary value problems are based on those for initial value problems we concentrate rst on initial value problems 63 Numerical Concepts In analyzing different procedures used in solving ordinary as well as partial differential equations we make use of the following numerical concepts 42 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Accuracy The order of a nite difference approximation of a differential equation is the rate at which the global error of the nite difference solution approaches zero as the sizes of the grid spacings approach zero This corresponds to the global truncation error There are two types of truncation errors amplitude errors and phase errors Sti 39ness An integration dif culty that can only occur in systems of equations or equations of order larger than rst whenever the ratio between the largest and smallest eigenvalue of the system is arge Stability When applied to a differential equation that has a bounded solution a nite difference equa tion is stable if it produces a bounded solution and is unstable if it produces an unbounded solution It is quite possible that the numerical solution to a differential equation grows unbounded even though its exact solution is well behaved Of course there are cases for which the exact solution grows unbounded but for our discussion of stability we concentrate only on the cases in which the exact solution is bounded Given the system of ODE s 62 and a numerical method in stability analysis we seek conditions in terms of the parameters of the numerical method mainly the step size h for which the numerical solution remains bounded In this context we have three classes of numerical methods stable numerical scheme the numerical solution does not blow up with any choice of step size unstable numerical scheme the numerical solution blows up with any choice of step size conditionally stable numerical scheme numerical solution remains bounded with certain choices of step size Convergence A nite difference method is convergent if the solution of the nite difference equation approaches a limit as the sizes of the grid spacings go to zero Note that we have no guarantee in general that this limit corresponds to the exact solution of the differential equation unless we know the solution in which case there would be no need for a numerical approximation Consistency A nite difference equation is consistent with a differential equation if the difference between the nite difference equation and the differential equation ie the truncation error goes to zero as the sizes of the grid spacings go to zero independently An important result which connects the last three concepts is known as the LaxRichtmyer or Lax equivalence theorem Given a properly posed initial value problem and a nite di erence approximation to it that satis es the consistency condition then stability is the necessary and su icient condition for the numerical solution to converge to the analytical solution 64 Model Differential Problem For convenience and feasibility of analytical treatment and without much loss in generality stability analysis will be performed on the model initial value di erential problem y Ay y0 yo 64 where A is a constant The solution of this model problem is y eA yO In our treatment of equation 64 we allow A to be complex A A R iAI with A3 g 0 The generalization of complex A will allow us to readily apply the results of our analysis to systems of ordinary differential equations and partial differential equations For example consider the second order differential equation y Lazy 0 65 MODEL DIFFERENCE PROBLEM 43 The exact solution is sinusoidal y cl coswt Cg sinwti We can convert this second order equation to two rst order equations Ziniz W32 The eigenvalues of the 2 gtlt 2 matrix are A iiwi Diagonalizing the matrix 0 1 A lt o gt A SAS l with the matrix of its eigenvectors7 S leads to the uncoupled set of equations where zSil y1 W 7 and A is the diagonal matrix with eigenvalues of A on the diagonal The uncoupled differential equations for the components of z are 21 iwzl z iiwzzi This simple example illustrates that higher order differential equations or systems of rst order differential equations can reduce to uncoupled ODE s of the form 64 with a complex coef cient The imaginary part of the coef cient results in oscillatory solutions of the form eizmt7 and the real part dictates whether the solution grows or decaysi As pointed out earlier7 in our stability analysis of equation 64 we will be concerned only with cases where A has zero or negative real part 65 Model Difference Problem Analogous to the differential problem7 consider the single rst order linear model initial value di e39rence pm em yn10yn 7101W 7 where yo is given and a is complex in generali The solution of this problem is yn anyo Note that the solution remains bounded only if lal g 1 44 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Im bh Region ofstability for exact solution 0 Re7 h Figure 61 Region of stability for exact solution 66 Stability of Numerical Methods The connection between the exact solution and the difference solution is evident if we evaluate the exact solution at tn nh for n 01Hl where h gt 0 then yn E ytn EMWO eknhyo anyo where a 6 Since we take our exact solution to be bounded then lal lekhl g 1 gt Reh ARI L g or That is in the AR17A h plane the region of stability of the exact solution is the left halfplane as shown in Figure l i Explicit or forward Euler method Applying the forward Euler method yn1 yn hftmyn to the model problem leads to yn1 y Ahyn 1 MOW Thus the solution at time step 71 can be written as yn 1 Ah y0l Note that a 1 ML The numerical solution is stable bounded if lall1 w l1 ARhiA1hl 1 ARI02 AI1 g 1 Note that only a small portion of the left halfplane which is the region of stability of the analytical solution is the region of stability for the explicit Euler method This region is inside the circle 1 ARI L 9111 1 66 STABILITY OF NUMERICAL METHODS 45 Region ofstabilil y for explicit Euler Re M Figure 62 Region of stability for explicit Euler method which is shown in Figure 62 For any value of Ah in the left halfplane and outside this circle the numerical solution blows up while the exact solution decaysl The numerical method is conditionally stablel To have a stable numerical solution we must reduce the step size 11 so that Ah falls within the circle If A is real and negative then the maximum step size for stability is 2 0lthlt l W The circle is only tangent to the imaginary axisl Therefore the forward Euler method is always unstable irrespective of the step size for purely imaginary A If A is real and the numerical solution is unstable then we must have 1 Ahl gt 1 which means that 1 Ah is negative with magnitude greater than 1 Since 111 1 Ah y07 the numerical solution exhibits oscillations with change of sign at every time step This behavior of the numerical solution is a good indication of instability Implicit or backward Euler method The implicit Euler scheme is given by the following formula ynl y39n hftnlyynl 65 Note that in contrast to the explicit Euler scheme we cannot very easily obtain the solution at the next time step If f is nonlinear we must solve a nonlinear algebraic equation at each time step to obtain ynirli Therefore the computational cost per time step for this scheme is much higher than that for the explicit Euler schemel However as we shall see below the implicit Euler scheme has better stability properties Applying backward Euler to the model equation 64 we obtain ynl yn Ahynr Solving for yn17 we have yn1 17 Ah 1ym yn anyoy 46 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS where 1 017 The denominator is a complex number and can be written as the product of its modulus and phase factor 7 1 a A6297 where 2111pm Re17h2Im1ih2 likgh21h221 since AR g 0 and 7 Im17h 7 Mb Qiarg Ah tan Re ikh tan lithi The modulus of a is 39 1 2 39 2 0 lcos6 zs1n6l cos 9s1n 9 1S 1 A A A Thus the backward Euler scheme is unconditionally stable and the region of stability in the AR27 A172 plane is identical to that of the exact solution When the region of stability of a difference equation is identical to the region of stability of the differential equation the nite difference scheme is sometimes referred to as Astablei Unconditional stability is the usual characteristic of implicit methods The price is higher computational cost per time step We emphasize that numerical stability does not imply accuracyi A method can be stable but inaccuratei From the stability point of view our objective is to use the maximum step size 12 to reach the nal destination at time t tfi Large time steps translate to lower number of function evaluations and lower computational cost This may not be the optimum h for acceptable accuracy but it is optimum for stability 67 Numerical Accuracy Revisited We have shown that the numerical solution to the model problem y Ay 66 is of the form yn yet 67 The exact solution is ytn eMquot yo WWW 618 One can often determine the order of accuracy of a method by comparing the numerical and exact solutions i equations 67 and 618 That is we compare the ampli cation factor a with 1 1 1 Ah 2 2 3 3 4 4 1 Ah h h h e 2 6 24 For example the ampli cation factor of the explicit Euler scheme is a1h 6 8 TRAPEZOIDAL METHOD 47 and the ampli cation factor for the backward Euler scheme is 1 1 Ah A2112 A3113 17 0 Thus both methods are able to reproduce only up to the Ah term in the exponential expansioni Each method is second order accurate for one time step but globally rst orderi From now on we will call a method ozLh order if its ampli cation factor matches all terms up to and including the Aah xal term in the exponential expansioni Often the order of accuracy by itself is not very informative In particular in problems with oscillatory solutions one is interested in the phase and amplitude errors separatelyi To get a handle on this type of error analysis we consider the model equation with A purely imaginary y iwy y0 1A The exact solution is eim which is oscillatoryi The frequency of oscillations is w and the amplitude is 1 The numerical solution with the explicit Euler scheme is yn 039 7 where a 1 iwhi a is a complex number and can be written as 0146 where Am1w2h2m217 and 9 tan 1 whi It is clear that lal Z 1 which recon rms that the Euler method is unstable with purely imaginary X A measure of the phase error is obtained from comparison with the exact solution qbn 9 7 wh tan 1 wh 7 whi Using the power series for tan l 71 1 3 1 5 1 7 tan wh wh 7 wh wh 7 wh 3 5 7 we have 1 n77 ltwhgt3 which corresponds to a phase lag This is the phase error encountered at each time step The phase error after 71 time steps is lt1 nab 68 I rapezoidal Method The formal solution to the nonlinear differential equation with the condition ytn yn is t W yn ft 7ydt A in At ttn1 n1 yn1 yn ft y WA in 48 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Approximating the integral with the trapezoidal method leads to h ynl yn 5 flttnlyynl flttnvynl This is the trapezoidal method for the solution of ordinary differential equations Applying the trapezoidal method to the model equation yields h yn1 yn 5 AZn1 Ayn 7 1 Ali2 yn1 7 1 7 All2y Expanding the ampli cation factor 0 leads to 71Ah2 1 1 7 1 Ah 2h2 431 13 17 xii2 2 4 which indicates that the method is second order accurate The extra accuracy is obtained at virtually no extra cost over the backward Euler schemei Now we will examine the stability properties of the trapezoidal method by computing the modulus of a for complex 1 A3 i1 71 ARhQ iA1h2 17 ARhZ 7 mil2 Both numerator and denominator are complex and can be written as Aei g and Beio respectively where 1 2 1 2 A 1 1h h 2Rgtlt21gtv and 1 2 1 2 B lt17 ARhgt 51111 1 Thus A A 20970 229 0 Be Be or 7 A 1017 Bi Since we are only interested in cases where AR g 0 and for these cases A g B it follows that lal g 11 Thus the trapezoidal method is unconditionally stable which is expected since it is an implicit methodi Note however that for real and negative 1 lim 0 71 ham which implies that the numerical solution anyo oscillates between 1 and 71 from time step to time step but the solution will not blowup1 69 RUN GE K UTTA METHODS 49 Let us examine the accuracy of the trapezoidal method for oscillatory solutions A iwi In this case AR 0 A B and lal 1 Thus there is no amplitude error associated with the trapezoidal methodi Since a 76 then h a e226 9 tan 1 the phase error is given by 7 71 7 l i 8 i 8 q5n72tan lt2 wh72 gm1 24wh whi 12wh which is about four times better than the explicit Euler schemei 69 Runge Kutta Methods The order of accuracy of numerical methods can be increased if one supplies additional information about the function f RungeKutta methods introduce points between tn and tn17 and evaluate f at these intermediate points The additional function evaluations of course result in higher cost per time step but the accuracy is increased and as it turns out better stability properties are also obtaine i We begin by describing the general form of second order RungeKutta formulas for solving y ay 69 The solution at time step tn1 is obtained from ynl yn Ylkl 72 610 Where the functions 61 and 162 are de ned sequentially k1 hftmyn 611 kg hftnahyn k1 6 12 and a B 71 and 72 are constants to be determined These constants are determined to ensure the highest order accuracy for the method To establish the order of accuracy consider the Taylor series expansion of ytn1 hi yn1 ynhyn3ynu But y flttnvyn7 and using the chain rule we have already obtained y ft f fyv Where ft and fy are the partial derivatives of f with respect to t and y respectively Thus hi gnu yn hflttmyn 30 My a 613 Twodimensional Taylor series expansion of k2 equation 61 leads to 162 hftnyn Fa1ftquot k1fyquot 50 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Noting that 161 hftnyn and substituting in 6110 yields ynl y39n 71 YZV39Lfn 72011213 Yz hzfnfy 39 39 H 614 Comparison of 6113 and 6114 and matching coef cients of similar terms leads to l l 7172 17 72a 57 72B 5 These are three nonlinear equations for the four unknownsi Using a as a free parameter we have 1 1 1 Zay 04 71 72 2a With the constants chosen we have a oneparameter family of second order RungeKutta formulas 7 1 1 k 1 k 6 15 ynl 7 ya 2a 1 2a 27 where 761 hftnvyn7 k2 hftnahynak1i Thus we have a second order RungeKutta formula for each arbitrary value of 11 The choice of a 12 is made frequentlyi RungeKutta formulas are often presented in a different but equivalent formi For example the popular form of the second order RungeKutta formula 1 12 is presented in the following predictorcorrector format yZHz yn ftn yn forward Euler predictor halfstep yn1 yn hftn12yl21 midpoint rule corrector fullstep Applying the RungeKutta method in equation 6115 to the model equation y Ay results in 761 hynv k2 hyna2hynAh1ahyn so that 17 1 h 1 Ah1 aAh 1 Ah 1A2h2 ynl yn 2a yn 2a yn 2 y39n Thus we have a con rmation that the method is second order accuratei For stability we must have 101 g 1 where 1 a 1 Ah 5A2h2gt1 A convenient way to satisfy this restriction is to set 1 2 2 2396 1 Ah 5A h e and nd the complex roots Ah of this polynomial for different values of 91 Note that 101 1 for all values of 91 The resulting stability region is shown in Figure 6131 On the real axis the stability boundary is the same as that of the explicit Euler scheme however there is an improvement for complex A1 The method is also unstable for purely imaginary 1 iwi In this case lalA11iw4h4Z11 69 RUNGE KUTTA METHODS 4H HHHHHH HHHH HH HHHHHH H 7 FounhOrder 7 27 7 2 7 7 5 E 7 7 27 7 SecondOrder 7 74 l l l 74 73 7 Ream 71 o 1 Figure 63 Stability regions of Runge Kutta schemes The phase error is 7 1 wh 7 qbn 7 tan 1 7 w2h22gt whi But7 1 w h 7 122 144 tan lt17w2h22gt 7 whlt12wh4wh 3 1 12 2 14 4 73whlt12wh4whn l whgw3h3 Hence7 1 ngw3h3 which is only a factor of two better than the explicit Euler scheme7 but of opposite sign The most Widely used Runge Kutta method is the fourth order formula 1 l 1 th1 7 ya Eh EUW k3 gin Where k1 hf hugn k2 hf lttWlhynlk1gt 2 2 k3 hf lttWlhynlk2gt 2 2 k4 hf tnhynk3i 52 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Note that four function evaluations are required at each time step The fourth order RungeKutta scheme can also be written in the following predictorcorrector format y12 yn ftn yn forward Euler predictor halfstep yZ fHz yn ftn12yl2 backward Euler corrector halfstep yij l yn hftn12y12 midpoint rule predictor fullstep yn1 yn ftm yn 2ftnl27y12 2ftn12y12 ftn1yffl Simpson s rule corrector fullstep Applying the method to the model equation leads to y 1 1 Ah lxihi lABI LAM L4 y t 2 6 24 7 which con rms the fourth order accuracy of the method Again the stability diagram is obtained by nding the roots of the following fourth order polynomial with complex coef cients 1 1 1 1 Ah 421 12 3h3 4h4 291 2 6 24 e For this purpose a root nder for polynomials with complex coef cients is required The region of stability is shown in Figure 6131 It shows a signi cant improvement over the second order RungeKutta scheme In particular it has a large stability region on the imaginary axisi 610 Similarities Between Linear Differential and Difference Equa I tions A second order linear homogeneous equation is given as dzu du 7 0 a2 dtz a1 dt aou 7 22142 21141 20140 1 If we introduce the differential operator D E ddt then D2 dZdtz etc and BM um or Dkum um Then the equation can be rewritten in operator form PDu azD2 a1D 120 u 01 If we assume a solution of the form then we obtain the characteristic polynomial PM 2M2 2M a0 01 The roots A1 and A2 can be found and then one can write the solution as u ek tcl ekztczi For linear homogeneous difference equations the solution procedure is very similari Say we have azun2 alun1 aoun 01 611 MULTI STEP METHODS 53 If we introduce the displacement or shift operator d 1 d2 hil39 2 Eie t 1hdt2h dt2 we note that Eun un1 and Eikun unik utn j Now in operator form the equation is PEun azE2 alE a0 un 0 and just as for the ODE if we introduce a solution of the form un a c we obtain the characteristic polynomial Pa 1202 11039 a0 0 Again the roots 01 and 02 can be found and the solution written as un aflcl agczi 611 MultiStep Methods In the case of RungeKutta formulas higher order accuracy was obtained by several function evaluationsi Higher order accuracy can also be achieved if one uses data from prior to tn that is if the solution y andor f at tn1 794 are used This is another way of providing additional information about fl Methods that use information from prior to step 71 are called multi step schemes The apparent price for the higher order accuracy is the use of additional computer memory Multi step methods are not selfstartingi Usually another method such as explicit Euler is used to start the calculations for the first or the first few time steps The most general mstep method is m m Zaiynl7i 7 thiftnl7i73nl7i 07 i0 i0 and is defined by the number of steps m and the parameters 22 and bza Without loss of generality we can take 10 1 lf 0 0 the method is explicit otherwise it is implicit Note that multi step algorithms require only one new function evaluation for each step There are many multistep methodsi They typically come in families 7 each member of the family possessing a different order of accuracy Three of the more popular families are AdamsBashforth Adams Moulton and Geari AdamsBashforth The Kth order is obtained by setting m K b0 0 a1 71 and 22 0i 2 i i i Ki The remaining coef cients are chosen to reproduce the Taylor expansion up to order Ki It is a Kstep algorithm and explicit since 120 0 The formulas for orders one through six are presented in Table 61 and their corresponding regions of stability are shown in Figure 64 AdamsMoulton The Kth order is obtained by setting m K 7 1 a1 71 and 22 0i 7 2iHK 7 1 Again the remaining coef cients are chosen to reproduce Taylor expansion up of order Ki The Kth order is a K 7 1step algorithm and is implicit since b0 0 The formulas for orders one through six are presented in Table 62 and their corresponding regions of stability are shown in Figure 65 Gear The Kth order is obtained by setting m K and b 0i 1 i i i Ki The remaining coef cients are chosen to reproduce Taylor expansion up of order Ki The Kth order is a Kstep algorithm and is implicit since b0 0 The formulas for orders one through six are presented in Table 63 and their corresponding regions of stability are shown in Figure 66 54 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS ORDER First yn1 yn hfam yn Second ynl yn h3ftn7y n flttnilvynill Third yn1 yn h23ftn7yni16flttn717yn71 5ftn2yn72 Fourth ynl yn hl55f myn 59ftn717yn71 37ftn727yn72 9ftn737yn73l Fifth ynh 7 yn h1901flttmyn 7 2774flttn7hymgt 2616flttn7hyn72gt 71274ftn73yn73 251flttn4 yum Sixth ynh 7 yn himwmmyn 7 7923flttn1ynhgt 9982flttn7hyn72gt 72 ftn787yn73 2877ftn747yn74 475ftn757yn75l Table 61 Adams Bashforth algorithms A classical multi step method is the Leapfrog method yn1 yn71 2hflttnvyn This method is derived by applying the second order central difference formula for y in equation 69 Applying the Leapfrog scheme to the model equation leads to yn1 17171 why This is a secondorder difference equation for yni To solve it7 we assume a solution of the form yn anal Substitution in the difference equation leads to a 1 anil ZhAani Dividing by an 1 since a 07 we get a quadratic equation for a 0272hka710i The solution is am 7 Ah j WI L2 1i Having more than one root is characteristic of multistep methods For comparison with the exponential solution7 we expand the roots in powers of z 1 1 01 Ahx2h211Ah A2h27 A4h4 02 AhixVhZi 771Ah7 xih2 vh47n The rst root shows that the method is second order accurate The second root is spurious7 and often is the source of numerical problems Note that for h 7gt 07 the spurious root approaches 71 It is also 611 MULTI STEP METHODS ORDER First ynl yn hftnlvynl Second ynl yn hlftnlvynl flttnvynl Third yn1 yn 11 2h5ftn17yn1 8ftnvyn ftn717yn71l Fourth yn1 yn 21 4h9ftn17y n119ft nvy n7 5ftn717yn71 flttn727yn72 Fifth ynl yn h251ftn1yyn1 646ftnvyn 264ftn717yn71 106ftn727 11172 19ftn737 yn73 Sixth ynl yn hl475f n1yyn1gt1427flt n7yn 798ftn717yn71 482ftn727 11172 173ftn737yn73 27ftn747 yn74 Table 62 Adams Moulton algorithmsi ORDER First yn1 yn hftn17y l Second yn1 4yn ynil 2hftnlvynll Third yn1 11 118yn 7 9yn71 211172 6hftnlvynll Fourth yn1 21 548yn 7 3mm 162 7 3W4 12hftn1yn1l Fifth W1 7 soomz 7 300ym 200 7 mm 12yH 60hftnlv yulgtl Sixth W1 7 smyn 7 45mm 400 7 225yH 724 flogn75 60hflttnlv yulgtl Table 63 Gear algorithmsi 56 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS apparent that for A real and negative the spurious root has modulus greater than 1 which leads to instability Since the difference equation for yn is linear its general solution can be written as a linear combination of the roots 239e yn Ulla Uilcz That is the solution is composed of contributions from both the physical and the spurious roots The constants cl and Cg are obtained from the starting conditions yo and y1 y061627 y101 10262 Solving for cl and Cg leads to 7 y1702yo 7 01yo7y1 c1 7 Cg 7 i 0391 i 0392 0391 i 0392 Thus for the model problem if we choose yl alyo the spurious root is completely suppressed In general we can expect the starting scheme to play a role in determining the level of contribution of the spurious rooti Even if the spurious root is suppressed initially roundoff errors can restart it again In the case of the Leapfrog scheme the spurious root leads to oscillations from one step to the next Application of the Leapfrog scheme to the case where A lad is purely imaginary leads to 013 iwh j V17 Luz12f If lwhl g 1 then l012l1 In this case the Leapfrog scheme does not have any amplitude errori This is the main reason for the use of the Leapfrog method If lwhl gt 1 then 012l whi xwzh2 71 7 and the method is unstable Finally we present the widely used second order AdamsBashforth method This method can be easily derived by using the Taylor series expansion of yn1l yn1 yn hy hzy hgy i Substituting y own and a rst order approximation for yx y39n y39n7l 7 flttnvy11gt7 flttn7lyyn7l yn 7 T 001 7 f 005 leads to 3 1 3 yn1 yn ghf myn hftn7lyyn7l 001 l The Adams Bashforth method is second order accurate globally Applying the method to the model problem leads to a second order difference equation 3 1 3 1 yn1 yn Ekhyn 7 EAhZn7l 1 5 yn 7 Ahyn4i 611 MULTI STEP METHODS t i i i E 7 2 ltlt ltlt E r E 1 A i E i 2 quot Remh 2 quot Remh 1 1 E 2 ltlt ltlt E r E 3 4 w it w 2 quot Remh 2 quot Remh t i i i E 2 ltlt ltlt E E 5 1 i 1 i 2 quot Remh 2 quot Remh Figure 64 Regions of stability for the AdamsBashforth algorithms for orders one through six 58 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS 4 7 4 7 2 2 5 7 A 5 7 g J g 1 7 735 Re h 0390 3 5 7 O 735 Re h 00 3 395 2 7 7 2 7 7 7 E z ltlt ltlt i i i i i i i i i i i i i i 2 2 3 7 735 Re h 00 3 5 7 O 735 Re h DD 3 t 4 7 4 7 z F 2 ltlt ltlt g J g 5 74 7 I A 7 I 770 735 Re h 00 3 5 770 735 Re h 00 3 Figure 65 Regions of stability for the AdamsMoulton algorithms for orders one through six 611 MULTI STEP METHODS 59 Im7th Im7th Remh Remh Im7th Im7th Remh Remh Im7th Im7th Remh Remh Figure 66 Regions of stability for the Gear algorithms for orders one through six 60 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Assuming a solution of the form yn one results in a quadratic equation for a 3 1 2 a lt1Zkhgtagkh 0 1 3 9 1 hj 1 Ah 2h2 1 012 2 lt 2 4 gt Using the power series expansion 9 1 9 1 9 2 3 9 3 1 Ah A2h2 1 Ah A2h2 7 Ah A2h2 Ah A2h2 7 4 2 lt 4 gt 8 lt 4 48 4 we obtain which has roots 01 1 Ah Aih2 0013 1 ilii 3 azigkh 2Ah Oh1 The spurious root for the AdamsBashforth method appears to be less dangerous Observe that it approaches zero if h gt 01 The stability boundary of the AdamsBashforth method can be obtained easily by taking a e in the characteristic polynomial for 01 Then we easily obtain that 219 7 Ah 2e 3 1 3e 7 1 and by plotting Ah as a function of 9 for 0 g 9 g 27f in the complex Reh7Imh plane we nd that the stability region is oval shaped It crosses the real axis at 71 which is more limiting than the explicit Euler and second order RungeKutta methods It is also only tangent to the imaginary axis1 Thus strictly it is unstable for pure imaginary A but it turns out that the instability is very mild 612 System of First Order Ordinary Differential Equations Recall that higher order ODE s can be converted to a system of rst order ODE s1 Systems of ODE s also appear naturally in many physical situations such as chemical reactions among several species A system of ODE s can be written in the generic form Y Wm y0 yo 616 where y is a vector with elements yz and fty1y2y3 1 1 1 ym is a vector function with elements fi1 From the application point of view the numerical solution of a system of ODE s is a straight forward extension of the techniques for a single ODE1 For example application of the explicit Euler scheme to equation 6116 yields ygnm y hfz39lttngtvy1nvygnvy nwu735773 7 i172737m The right hand side can be calculated using data from the previous time step and each equation can be advanced forward1 From the conceptual point of view there is only one fundamental difference between the numerical solution of one ODE and that of a system1 This is the stiffness property that leads to some numerical 612 SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS 61 problems in systems but it is not an issue with single ODE si We shall discuss stiffness in connection with the system of equations with constant coef cients Ay yltogt yo 617 where A is an m gtlt m constant matrixi Equation 617 is the model problem for systems of ODE si In the same manner that the single equation model problem was helpful in analyzing numerical methods for a single ODE problem 617 is useful for analyzing numerical methods for systems From linear algebra we know that this system will have a bounded solution if all eigenvalues of A have negative real parts This is analogous to the one equation model problem y Ay where the real part of A was negative Applying the explicit Euler method to 617 leads to yn1 yn hAyn I M ym OI yn I 71A yo Buyoy where I is the identity matrix To have a stable numerical solution for large n the matrix B should be bounded For this to happen we know from linear algebra that the eigenvalues of B must have modulus less than or equal to unity The eigenvalues of B are a 1 hAi where the Ai s are the eigenvalues of the matrix A Thus for numerical stability we must have 1 hkil S 1 The eigenvalue having the largest magnitude places the most restriction on hi If the eigenvalues are real and negative then 2 7 Vilma If the range of the magnitudes of the eigenvalues is large 239e S E max Milmin gtgt 1 where S is called the stiffness ratio and the solution is desired over a large span of the independent variable t then the system of differential equations is called a stz systeml Stiffness arises in physical situations with many degrees of freedom but with widely different rates of response Examples include a system composed of two springs one very stiff and the other very exible a mixture of chemical species with very different reaction rates and a boundary layer with two disparate length scales Stiff systems are associated with numerical dif culties Problems arise if the system of equations are to be integrated to large values of the independent variable t Since the step size is limited by the part of the solution with fastest response time ie with largest eigenvalue magnitude the number of steps required can become enormousl In other words even if one is interested only in the long term behavior of the solution the time step must still be very small In practice to circumvent stiffness implicit methods are used With implicit methods there is no restriction on time step due to numerical stabilityl For high accuracy one can choose small time steps to resolve the rapidly varying portions of the solution fast parts and large time steps in the slowly varying portionsi There are stiff ODE solvers that have an adaptive time step selection mechanisml These are based on implicit methods and automatically reduce or increase the time step depending on the behavior of the solution Note that with explicit methods one cannot use large time steps in the slowly varying part of the solution Roundoff errors will trigger the numerical instability associated with the fast part of the solution Example Consider the uncoupled twodimensional linear system li lefillZZl 62 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS H l we want to obtain the solution by numerical integration To be speci c assume that Af 7106 and A5 71 A conservative estimate of the length of the transient is ten time constants for a rst order equation the time constant is Therefore yft m 0 for t gt tf m 10 s and y5t m 0 for t gt its m 10 Suppose that when either yf or ys is integrated alone using the forward Euler scheme accuracy considerations require ten time steps per time constant Thus the time step is hf 10 7 for the fast equation and hs 10 1 for the slow onei Note that the values that we have chosen individually satisfy the stability restrictions These values suggest that when the system is being integrated a step size of 10 7 should be used when 0 lt t lt tf and 10 1 during the interval tf lt t lt tsi Furthermore the integration is stopped when t gt its since then yf m 0 and ys m 0 Using these step sizes the total number of steps is with exact solution 6M ekst yfo yso Mt ys t 10 7 105 t5 7 tf 7 105 1071 7 1077 m 200 tf N hf h We can investigate the stability of the forward Euler scheme by calculating Ah for the four combina tions of eigenvalues and step sizesi During 0 lt t lt tf we have Afhf 710 1 and Ashf 710 7i Both these values lie within the circle of stability of the forward Euler schemei During tf lt t lt t5 we have A5115 710 1 and Afhs 7105 The value 7105 lies well outside the region of stability and thus the forward Euler scheme is unstable using these step sizes to solve the systemi Since both eigenvalues are real the region of stability for the forward Euler scheme limits the step sizes to h g ZWma xi Thus for our case hma x ZMfl 2 gtlt 10 6i Then the number of steps needed to integrate the sytem is tf t5 7 tf 105 10 7105 N m 10 i hf hm 1m7 1m6 Assuming that each integration step takes 10 3 s of CPU time the coupled integration takes 28 hrs while separate integrations would only take 02 s Not only the integration of the system takes an unreasonably long time but the cumulative effect of the roundoff error after 107 steps could be quite large Now suppose that the backward Euler scheme is used instead Then we have no stability limit on the step size Thus we can use hf and hs in the separate regions of time and thus integrate the system ef ciently The problem with implicit methods is that in general at each time step they require the solution of a nonlinear algebraic equation which often requires an iterative solution procedure such as the Newton Raphson method For nonlinear differential equations iteration can be avoided by the linean39zation technique We shall rst demonstrate the technique for the single ordinary differential equation y ft7yA Applying the trapezoidal method to this equation yields To solve for yn1 would require solving a nonlinear algebraic equationi yn1 yn 3 WWW fawn 0013 618 Consider the Taylor series expansion of ftn1yn1 about the solution yn 10 0f 2f flttnlvynl flttnlvyn Ty yquot yn1 yn 2 0 y2 yquot yn1 yn2 39 39 quot 6 19 613 BOUNDARY VALUE PROBLEMS 63 But from 618 we see that ynl yn 00quot therefore replacing ftn1 yn1 in equation 618 with the rst two terms in its Taylor series expansion does not alter the order of accuracy of 618 which for one step is Ohgi Making this substitution results in 0013 6 20 h 0 yn1 yn 5 flttnlvyn 0 yn1 7 yn ftmyn yn Rearranging and solving for yn1 yields ftn17yn ftmynl lih h yn1 yn 2 2 3y yquot Thus the solution can proceed without iteration and with retention of global second order accuracyi Linearization can also be applied to the implicit solution of a system of ODE si Consider the system dy f t m dt 71173127 y where bold letters again are used for vectorsi Applying the trapezoidal method results in h yltn1gt yltngt 5 f lttltn1gt7yltn1gtgt f lttltngt7yltngt Taylor series expansion of f with elements fi about the solution ylt gt yields m 1 tltn1gtyltn1gt fitltn1gtyltngt Z n1 W 2 ayj yj y gtOh z 12Himi n 11 y In matrix notation we have f lttltn1gt7yltn1gtgt f lttltn1gt7yltngtgt Am ltyltn1gt yltngtgt 002 where 2amp 2A 2L 3 Byz Bym AW am am 2m 3 Byz Bym n is the Jacobian matrix At each time step one must solve the following system of algebraic equations h h 1 yltn1gt yltngt 5 I Two f lttltn1gt7yltngtgt f lttltngt7yltngt Note that the matrix A is not constant in general and should be updated at every time step 613 Boundary Value Problems When data associated with a differential equation are prescribed at more than one value of the inde pendent variable then the problem is a boundary value problemi ln initial value problems all the data y0y0 i i are prescribed at one value of the independent variable in this case at t 0 To have a boundary value problem we must have at least a second order differential equation y fryy y0 yo yL my 621 where f is an arbitrary function There are two techniques for solving boundary value problems 64 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS H i Shooting Method Shooting is a technique which uses standard methods for initial value problems such as Runge Kutta metho sf 5 Direct Methods These methods are based on straight forward discretization of the derivatives in the differential equation and solving the resulting system of algebraic equations 6131 Shooting Method We shall begin with the discussion of the shooting methodi Let s reduce the second order ODE in 621 to two first order equations by letting u y and v y so that u U 1 f17uyv7 with end conditions 140 yo uL yLi To solve this system one needs one condition for each of the unknowns u and 1 rather than two for one and none for the otheri Therefore we use a guess for v0 and integrate the equations to z Li At this point uL is compared to yL if the agreement is not satisfactory most likely it will not be another guess is made for v0 and the iterative process is repeated For linear problems this iterative process is very systematic only two iterations are needed Consider the general linear equation WI AIy r MIMI HI 6322 y0 yo yL yL Let s denote two solutions of the equation as y1z and y2z which were obtained using two initial guesses for y 0i Since the differential equation is linear the exact solution can be formed as a linear combination of y1 and y2 MI C1y11 Czyzfx Since y10 y20 y0 it follows that cl Cg 1 Next we require that yL yL which in turn requires that C1y1L Czysz yL We have two equations for cl and Cg the solution is C yL y2L C y1LyL 1 y1ltLgtey2ltLgt 2 y1ltLgtey2ltLgt This approach also works for higher order linear systems In general if 71 conditions are specified at the final point 71 1 solutions must be generated to obtain the final solution by superpositioni For this reason if the number of conditions at the two endpoints are unequal it is more ef cient to start the integration from the side where the larger number of conditions are given Unfortunately when f is nonlinear we may have to perform several iterations to obtain the solution at L to within a prescribed accuracyi Any of the methods for solving a system of nonlinear equations can be applied see Chapter 2 One may encounter dif culty in obtaining a converged solution if yL is a very sensitive function of y 0i In this case it may be worthwhile to integrate from the opposite direction by guessing a value of yL and iterating until y0 is suf ciently close to yoi 613 BOUNDARY VALUE PROBLEMS 65 6132 Direct Methods With direct methods one simply approximates the derivatives in the differential equation with a nite difference or some other approximation of suitable accuracyl One obtains a system of algebraic equations for the dependent variables at the node points For linear differential equations the system is a linear system of algebraic equations for nonlinear equations it is a nonlinear system of algebraic equations Second order approximation to the linear equation 6 22 yields y39172y39ygt1 y39 liy39il Aj TJBjyjfj 112MN71 yj0 yo yjN yLy where a uniform grid 11 zj 7 zj1 is introduced between the boundary points 10 and IN Rearranging the terms yields ajyjl jyj7jyjilsz 11727 A 7Nily 1 Aj 2 1 Aj aFltW gtv BFltBFWgtV 7FltW l This is a tridiagonal system of linear algebraic equations The only special treatment comes at the points next to the boundaries j 1 and j N7 1 t j 1 we have where mm lyl f1 YlyoA Note that yo which is known appears on the right hand side Similarly yN appears on the right hand side of the corresponding equation for the right boundaryl Thus the unknowns y1y2 l l l yN1 are obtained from the solution of 51 a1 y1 f1 71y0 72 52 a2 yz f2 7N4 BN4 O Niij LyNiij L fNii j 39YNil 5N4 yNil fNil aNilyL Implementation of mixed boundary conditions such as ay0 by 0 g is also straight forward One simply approximates yO with a onesided nite difference approximation such as yO 3yo 4y1 yz and solves for yo in terms of y1 y2 and 9 The result is then substituted in the nite difference equation evaluated at j 1 Higher order nitedifference approximations can also be used The only dif culty with higher order methods is that near the boundaries they require data from points outside the domain The standard procedure is to use either lower order or onesided approximations for points near the boundaries Higher order nitedifferences lead to broad band matrices instead of a tridiagonal matrix For example a pentadiagonal system is obtained with a fourth order central difference approximation to equation 62 66 CHAPTER 6 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS 614 Eigenvalue Problems Eigenvalue problems are simply the homogeneous versions of boundary value problems with the presence of a parameter To illustrate the solution procedures consider the following prob em y kzy 0 y0 y1 0 whose eigenvalues are kn nTr 1 4 121 and eigenvectors sinn7rzi 6141 Shooting Method Letting 141 y and 142 y we can rewrite the equation as 141 7 0 1 141 142 7 7162 0 142 7 with boundary conditions 1410 1411 0 Since the system is second order it possesses 2 linearly independent solutions that we call 14Ejgtj 121 The general solution for 14239 will be 2 V 14 E cjug i F1 The essence of the method is to replace the boundary value problem by an initial value problemi The idea is to use the initial conditions for the 14Jgt which later simplify the problem considerably namely 14200 6 j121 Then the rst boundary condition gives 1410 cluglgt0 6214 0 clugl0 0 gt cl 0 and the second boundary condition gives 7411 W490 of If we want a nontrivial solution we must have 1491 01 Now the system is numerically integrated with the initial condition 1490 5427 i we take 1410 0 and 1420 l and a xed value of kl From the integration we nd that 1491 01 Generally a 0 in which case we change the value of k and repeat the process until 1491 m 0 New approximations for k are obtained by standard iteration methodsi Note that the eigenfunction if we um c2u 2gtltzgt 615 NONUNIFORM GRIDS 67 6142 Direct Method To illustrate the method we use secondorder central differences on a uniform grid to discretise the equation and obtain 4 7 2y39 y39 1 k2yi0 y0yN0 212HlN71l The above equations can be rewritten in the following form A 7 my o where 111 2 1 y2 1 72 1 y 7 A 7 yNii 1 2 1 yNil 1 2 I is the identity matrix and A iszzzl The above system is a standard eigenvalue problem in linear algebra and can be solved by standard methods for solving such problems 615 NonUniform Grids Often the solution of a differential equation varies rapidly in a part of the domain and it has mild variation elsewhere In such cases it is wasteful to use a ne grid capable of resolving the rapid variations everywhere in the domain One should use a nonuniform grid spacing In some problems such as boundary layers in uid ow problems the regions of rapid variations are known a pn39o m39 and grid points can be clustered where needed There are also adaptive techniques that estimate the grid requirements as the solution progresses and place additional grid points in the regions of rapid variations With nonuniform grids one can either use nitedifference formulas for nonuniform grids or use a coordinate transformation Typical nite difference formulas for rst and second derivatives are yj1 yjil 1 V 7 1 71177 y hjlhj 2 11 y 7 7 yj1 yj y 1 V 7 V m 7 I I y 2 15141011 hm hjhj1 MW WM 3 h l My 7 where hj zj 7 zj71l These formulas can be substituted in equation 62 and the resulting system of equations can be solved It should be noted that nite difference formulas for nonuniform meshes gener ally have a lower order of accuracy than their counterparts with the same stencil but de ned for uniform meshes The lower order accuracy arises because of reduced cancellations in Taylor series expansions due to lack of symmetry of the meshesl
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'