Lecture notes introduction to Ordinary differential Equations

lecture notes on numerical solution of ordinary differential equations and lecture notes on theory of ordinary differential equations pdf free download
GregDeamons Profile Pic
GregDeamons,New Zealand,Professional
Published Date:03-08-2017
Your Website URL(Optional)
cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 547 PART SIX Ordinary Differential Equations 6.1 OVERVIEW The fundamental laws of physics, mechanics, electricity, and thermodynamics are usually based on empirical observations that explain variations in physical properties and states of systems. Rather than describing the state of physical systems directly, the laws are usually couched in terms of spatial and temporal changes. These laws define mechanisms of change. When combined with continuity laws for energy, mass, or momentum, differential equations result. Subsequent integration of these differential equations results in mathe- matical functions that describe the spatial and temporal state of a system in terms of energy, mass, or velocity variations. As in Fig. PT6.1, the integration can be implemented analyti- cally with calculus or numerically with the computer. The free-falling bungee jumper problem introduced in Chap. 1 is an example of the de- rivation of a differential equation from a fundamental law. Recall that Newton’s second law was used to develop an ODE describing the rate of change of velocity of a falling bungee jumper: dv c d 2 = g − v (PT6.1) dt m where g is the gravitational constant, m is the mass, and c is a drag coefficient. d Such equations, which are composed of an unknown function and its derivatives, are called differential equations. They are sometimes referred to as rate equa- tions because they express the rate of change of a variable as a function of vari- ables and parameters. In Eq. (PT6.1), the quantity being differentiated v is called the dependent variable. The quantity with respect to which v is differentiated t is called the independent variable. When the function 547cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 548 548 PART 6 ORDINARY DIFFERENTIAL EQUATIONS Physical law F = ma dv c d 2 ODE = g − v dt m Analytical Numerical (calculus) (computer)       gm gc c d d 2 Solution v = tanh t v = v + g − v t i+1 i c m m d FIGURE PT6.1 The sequence of events in the development and solution of ODEs for engineering and science. The example shown is for the velocity of the free-falling bungee jumper. involves one independent variable, the equation is called an ordinary differential equation (or ODE). This is in contrast to a partial differential equation (or PDE) that involves two or more independent variables. Differential equations are also classified as to their order. For example, Eq. (PT6.1) is called a first-order equation because the highest derivative is a first derivative. A second- order equation would include a second derivative. For example, the equation describing the position x of an unforced mass-spring system with damping is the second-order equation: 2 d x dx m + c + kx = 0 (PT6.2) 2 dt dt where m is mass, c is a damping coefficient, and k is a spring constant. Similarly, an nth- order equation would include an nth derivative. Higher-order differential equations can be reduced to a system of first-order equations. This is accomplished by defining the first derivative of the dependent variable as a new variable. For Eq. (PT6.2), this is done by creating a new variable v as the first derivative of displacement dx v = (PT6.3) dt where v is velocity. This equation can itself be differentiated to yield 2 dv d x = (PT6.4) 2 dt dt Equations (PT6.3) and (PT6.4) can be substituted into Eq. (PT6.2) to convert it into a first- order equation: dv m + cv + kx = 0 (PT6.5) dtcha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 549 6.1 OVERVIEW 549 As a final step, we can express Eqs. (PT6.3) and (PT6.5) as rate equations: dx = v (PT6.6) dt dv c k =− v − x (PT6.7) dt m m Thus, Eqs. (PT6.6) and (PT6.7) are a pair of first-order equations that are equivalent to the original second-order equation (Eq. PT6.2). Because other nth-order differential equations can be similarly reduced, this part of our book focuses on the solution of first- order equations. A solution of an ordinary differential equation is a specific function of the independent variable and parameters that satisfies the original differential equation. To illustrate this concept, let us start with a simple fourth-order polynomial, 4 3 2 y =−0.5x + 4x − 10x + 8.5x + 1 (PT6.8) Now, if we differentiate Eq. (PT6.8), we obtain an ODE: dy 3 2 =−2x + 12x − 20x + 8.5 (PT6.9) dx This equation also describes the behavior of the polynomial, but in a manner different from Eq. (PT6.8). Rather than explicitly representing the values of y for each value of x, Eq. (PT6.9) gives the rate of change of y with respect to x (i.e., the slope) at every value of x. Figure PT6.2 shows both the function and the derivative plotted versus x. Notice how the zero values of the derivatives correspond to the point at which the original function is flat—that is, where it has a zero slope. Also, the maximum absolute values of the deriva- tives are at the ends of the interval where the slopes of the function are greatest. Although, as just demonstrated, we can determine a differential equation given the original function, the object here is to determine the original function given the differential equation. The original function then represents the solution. Without computers, ODEs are usually solved analytically with calculus. For example, Eq. (PT6.9) could be multiplied by dx and integrated to yield  3 2 y = (−2x + 12x − 20x + 8.5) dx (PT6.10) The right-hand side of this equation is called an indefinite integral because the limits of in- tegration are unspecified. This is in contrast to the definite integrals discussed previously in Part Five compare Eq. (PT6.10) with Eq. (19.5). An analytical solution for Eq. (PT6.10) is obtained if the indefinite integral can be eval- uated exactly in equation form. For this simple case, it is possible to do this with the result: 4 3 2 y =−0.5x + 4x − 10x + 8.5x + C (PT6.11) which is identical to the original function with one notable exception. In the course of dif- ferentiating and then integrating, we lost the constant value of 1 in the original equation and gained the value C. This C is called a constant of integration. The fact that such an arbitrary constant appears indicates that the solution is not unique. In fact, it is but one ofcha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 550 550 PART 6 ORDINARY DIFFERENTIAL EQUATIONS y 4 x 3 (a) dy/dx 8 3 x –8 (b) FIGURE PT6.2 4 3 2 Plots of (a) y versus x and (b) dy/dx versus x for the function y=−0.5x + 4x − 10x + 8.5x + 1. an infinite number of possible functions (corresponding to an infinite number of possible values of C) that satisfy the differential equation. For example, Fig. PT6.3 shows six pos- sible functions that satisfy Eq. (PT6.11). Therefore, to specify the solution completely, a differential equation is usually accom- panied by auxiliary conditions. For first-order ODEs, a type of auxiliary condition called an initial value is required to determine the constant and obtain a unique solution. For example, the original differential equation could be accompanied by the initial condition that at x = 0, y = 1. These values could be substituted into Eq. (PT6.11) to determine C = 1. Therefore, the unique solution that satisfies both the differential equation and the specified initial condition is 4 3 2 y =−0.5x + 4x − 10x + 8.5x + 1 Thus, we have “pinned down” Eq. (PT6.11) by forcing it to pass through the initial condi- tion, and in so doing, we have developed a unique solution to the ODE and have come full circle to the original function Eq. (PT6.8). Initial conditions usually have very tangible interpretations for differential equations derived from physical problem settings. For example, in the bungee jumper problem, thecha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 551 6.2 PART ORGANIZATION 551 y C =3 C =2 C =1 C =0 x C =–1 C =–2 FIGURE PT6.3 3 2 Six possible solutions for the integral of −2x + 12x − 20x + 8.5. Each conforms to a different value of the constant of integration C. initial condition was reflective of the physical fact that at time zero the vertical velocity was zero. If the bungee jumper had already been in vertical motion at time zero, the solu- tion would have been modified to account for this initial velocity. When dealing with an nth-order differential equation, n conditions are required to ob- tain a unique solution. If all conditions are specified at the same value of the independent variable (e.g., at x or t = 0), then the problem is called an initial-value problem. This is in contrast to boundary-value problems where specification of conditions occurs at different values of the independent variable. Chapters 22 and 23 will focus on initial-value prob- lems. Boundary-value problems are covered in Chap. 24. 6.2 PART ORGANIZATION Chapter 22 is devoted to one-step methods for solving initial-value ODEs. As the name suggests, one-step methods compute a future prediction y , based only on information at i+1 a single point y and no other previous information. This is in contrast to multistep ap- i proaches that use information from several previous points as the basis for extrapolating to a new value. With all but a minor exception, the one-step methods presented in Chap. 22 belong to what are called Runge-Kutta techniques. Although the chapter might have been organized around this theoretical notion, we have opted for a more graphical, intuitive approach to in- troduce the methods. Thus, we begin the chapter with Euler’s method, which has a very straightforward graphical interpretation. In addition, because we have already introduced Euler’s method in Chap. 1, our emphasis here is on quantifying its truncation error and de- scribing its stability. cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 552 552 PART 6 ORDINARY DIFFERENTIAL EQUATIONS Next, we use visually oriented arguments to develop two improved versions of Euler’s method—the Heun and the midpoint techniques. After this introduction, we formally de- velop the concept of Runge-Kutta (or RK) approaches and demonstrate how the foregoing techniques are actually first- and second-order RK methods. This is followed by a discus- sion of the higher-order RK formulations that are frequently used for engineering and scientific problem solving. In addition, we cover the application of one-step methods to systems of ODEs. Note that all the applications in Chap. 22 are limited to cases with a fixed step size. In Chap. 23, we cover more advanced approaches for solving initial-value problems. First, we describe adaptive RK methods that automatically adjust the step size in response to the truncation error of the computation. These methods are especially pertinent as they are employed by MATLAB to solve ODEs. Next, we discuss multistep methods. As mentioned above, these algorithms retain in- formation of previous steps to more effectively capture the trajectory of the solution. They also yield the truncation error estimates that can be used to implement step-size control. We describe a simple method—the non-self-starting Heun method—to introduce the essential features of the multistep approaches. Finally, the chapter ends with a description of stiff ODEs. These are both individual and systems of ODEs that have both fast and slow components to their solution. As a con- sequence, they require special solution approaches. We introduce the idea of an implicit solution technique as one commonly used remedy. We also describe MATLAB’s built-in functions for solving stiff ODEs. In Chap. 24, we focus on two approaches for obtaining solutions to boundary-value problems: the shooting and finite-difference methods. Aside from demonstrating how these techniques are implemented, we illustrate how they handle derivative boundary conditions and nonlinear ODEs.cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 553 22 Initial-Value Problems CHAPTER OBJECTIVES The primary objective of this chapter is to introduce you to solving initial-value problems for ODEs (ordinary differential equations). Specific objectives and topics covered are Understanding the meaning of local and global truncation errors and their • relationship to step size for one-step methods for solving ODEs. Knowing how to implement the following Runge-Kutta (RK) methods for • a single ODE: Euler Heun Midpoint Fourth-order RK Knowing how to iterate the corrector of Heun’s method. • Knowing how to implement the following Runge-Kutta methods for systems • of ODEs: Euler Fourth-order RK YOU’VE GOT A PROBLEM e started this book with the problem of simulating the velocity of a free-falling bungee jumper. This problem amounted to formulating and solving an ordinary W differential equation, the topic of this chapter. Now let’s return to this problem and make it more interesting by computing what happens when the jumper reaches the end of the bungee cord. 553cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 554 554 INITIAL-VALUE PROBLEMS To do this, we should recognize that the jumper will experience different forces de- pending on whether the cord is slack or stretched. If it is slack, the situation is that of free fall where the only forces are gravity and drag. However, because the jumper can now move up as well as down, the sign of the drag force must be modified so that it always tends to retard velocity, dv c d 2 = g − sign(v) v (22.1a) dt m 2 where v is velocity (m/s), t is time (s), g is the acceleration due to gravity (9.81 m/s ), c is d 1 the drag coefficient (kg/m), and m is mass (kg). The signum function, sign, returns a −1 or a 1 depending on whether its argument is negative or positive, respectively. Thus, when the jumper is falling downward (positive velocity, sign = 1), the drag force will be negative and hence will act to reduce velocity. In contrast, when the jumper is moving upward (negative velocity, sign =−1), the drag force will be positive so that it again reduces the velocity. Once the cord begins to stretch, it obviously exerts an upward force on the jumper. As done previously in Chap. 8, Hooke’s law can be used as a first approximation of this force. In addition, a dampening force should also be included to account for frictional effects as the cord stretches and contracts. These factors can be incorporated along with gravity and drag into a second force balance that applies when the cord is stretched. The result is the following differential equation: dv c k γ d 2 = g − sign(v) v − (x − L) − v (22.1b) dt m m m where k is the cord’s spring constant (N/m), x is vertical distance measured downward from the bungee jump platform (m), L is the length of the unstretched cord (m), and γ is a damp- ening coefficient (N · s/m). Because Eq. (22.1b) only holds when the cord is stretched (x L), the spring force will always be negative. That is, it will always act to pull the jumper back up. The damp- ening force increases in magnitude as the jumper’s velocity increases and always acts to slow the jumper down. If we want to simulate the jumper’s velocity, we would initially solve Eq. (22.1a) until the cord was fully extended. Then, we could switch to Eq. (22.1b) for periods that the cord is stretched. Although this is fairly straightforward, it means that knowledge of the jumper’s position is required. This can be done by formulating another differential equa- tion for distance: dx = v (22.2) dt Thus, solving for the bungee jumper’s velocity amounts to solving two ordinary dif- ferential equations where one of the equations takes different forms depending on the value 1 Some computer languages represent the signum function as sgn(x). As represented here, MATLAB uses the nomenclature sign(x).cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 555 22.2 EULER’S METHOD 555 of one of the dependent variables. Chapters 22 and 23 explore methods for solving this and similar problems involving ODEs. 22.1 OVERVIEW This chapter is devoted to solving ordinary differential equations of the form dy = f (t, y) (22.3) dt In Chap. 1, we developed a numerical method to solve such an equation for the velocity of the free-falling bungee jumper. Recall that the method was of the general form New value = old value + slope × step size or, in mathematical terms, y = y + φh (22.4) i+1 i where the slope φ is called an increment function. According to this equation, the slope es- timate of φ is used to extrapolate from an old value y to a new value y over a distance i i+1 h. This formula can be applied step by step to trace out the trajectory of the solution into the future. Such approaches are called one-step methods because the value of the increment function is based on information at a single point i. They are also referred to as Runge- Kutta methods after the two applied mathematicians who first discussed them in the early 1900s. Another class of methods called multistep methods use information from several previous points as the basis for extrapolating to a new value. We will describe multistep methods briefly in Chap. 23. All one-step methods can be expressed in the general form of Eq. (22.4), with the only difference being the manner in which the slope is estimated. The simplest approach is to use the differential equation to estimate the slope in the form of the first derivative at t . In i other words, the slope at the beginning of the interval is taken as an approximation of the average slope over the whole interval. This approach, called Euler’s method, is discussed next. This is followed by other one-step methods that employ alternative slope estimates that result in more accurate predictions. 22.2 EULER’S METHOD The first derivative provides a direct estimate of the slope at t (Fig. 22.1): i φ = f (t ,y ) i i where f (t , y ) is the differential equation evaluated at t and y . This estimate can be sub- i i i i stituted into Eq. (22.1): y = y + f (t ,y )h (22.5) i+1 i i i This formula is referred to as Euler’s method (or the Euler-Cauchy or point-slope method). A new value of y is predicted using the slope (equal to the first derivative at the original value of t) to extrapolate linearly over the step size h (Fig. 22.1).cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 556 556 INITIAL-VALUE PROBLEMS y Predicted error True h t t t i i1 FIGURE 22.1 Euler’s method. EXAMPLE 22.1 Euler’s Method  0.8t Problem Statement. Use Euler’s method to integrate y = 4e − 0.5y from t = 0 to 4 with a step size of 1. The initial condition at t = 0 is y = 2. Note that the exact solution can be determined analytically as 4 0.8t −0.5t −0.5t y = (e − e ) + 2e 1.3 Solution. Equation (22.5) can be used to implement Euler’s method: y(1) = y(0) + f (0, 2)(1) where y(0) = 2 and the slope estimate at t = 0 is 0 f (0, 2) = 4e − 0.5(2) = 3 Therefore, y(1) = 2 + 3(1) = 5 The true solution at t = 1 is  4 0.8(1) −0.5(1) −0.5(1) y = e − e + 2e = 6.19463 1.3 Thus, the percent relative error is 6.19463 − 5 ε = × 100% = 19.28% t 6.19463 For the second step: y(2) = y(1) + f (1, 5)(1) 0.8(1) = 5 + 4e − 0.5(5) (1) = 11.40216cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 557 22.2 EULER’S METHOD 557 TABLE 22.1 Comparison of true and numerical values of the integral 0.8t of y’ = 4e − 0.5y, with the initial condition that y = 2 at t = 0. The numerical values were computed using Euler’s method with a step size of 1. ty y ε (%) true Euler t 0 2.00000 2.00000 1 6.19463 5.00000 19.28 2 14.84392 11.40216 23.19 3 33.67717 25.51321 24.24 4 75.33896 56.84931 24.54 y 60 40 True solution 20 Euler solution 0 t 01234 FIGURE 22.2 Comparison of the true solution with a numerical solution using Euler’s method for the integral of 0.8t y’ = 4e − 0.5y from t = 0 to 4 with a step size of 1.0. The initial condition at t = 0 is y = 2. The true solution at t = 2.0 is 14.84392 and, therefore, the true percent relative error is 23.19%. The computation is repeated, and the results compiled in Table 22.1 and Fig. 22.2. Note that although the computation captures the general trend of the true solution, the error is considerable. As discussed in the next section, this error can be reduced by using a smaller step size. 22.2.1 Error Analysis for Euler’s Method The numerical solution of ODEs involves two types of error (recall Chap. 4): 1. Truncation, or discretization, errors caused by the nature of the techniques employed to approximate values of y. 2. Roundoff errors caused by the limited numbers of significant digits that can be retained by a computer. The truncation errors are composed of two parts. The first is a local truncation error that results from an application of the method in question over a single step. The second is a propagated truncation error that results from the approximations produced during thecha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 558 558 INITIAL-VALUE PROBLEMS previous steps. The sum of the two is the total error. It is referred to as the global trunca- tion error. Insight into the magnitude and properties of the truncation error can be gained by de- riving Euler’s method directly from the Taylor series expansion. To do this, realize that the differential equation being integrated will be of the general form of Eq. (22.3), where  dy/dt = y , and t and y are the independent and the dependent variables, respectively. If the solution—that is, the function describing the behavior of y—has continuous derivatives, it can be represented by a Taylor series expansion about a starting value (t , y ), as in i i recall Eq. (4.13): (n)  y y  i 2 i n y = y + y h + h +···+ h + R (22.6) i+1 i n i 2 n where h = t − t and R = the remainder term, defined as i+1 i n (n+1) y (ξ) n+1 R = h (22.7) n (n + 1) where ξ lies somewhere in the interval from t to t . An alternative form can be devel- i i+1 oped by substituting Eq. (22.3) into Eqs. (22.6) and (22.7) to yield  (n−1) f (t , y ) f (t , y ) i i i i 2 n n+1 y = y + f (t , y )h + h +···+ h + O(h ) (22.8) i+1 i i i 2 n n+1 where O(h ) specifies that the local truncation error is proportional to the step size raised to the (n + 1)th power. By comparing Eqs. (22.5) and (22.8), it can be seen that Euler’s method corresponds to the Taylor series up to and including the term f (t , y )h. Additionally, the comparison i i indicates that a truncation error occurs because we approximate the true solution using a fi- nite number of terms from the Taylor series. We thus truncate, or leave out, a part of the true solution. For example, the truncation error in Euler’s method is attributable to the remain- ing terms in the Taylor series expansion that were not included in Eq. (22.5). Subtracting Eq. (22.5) from Eq. (22.8) yields  f (t , y ) i i 2 n+1 E = h +···+ O(h ) (22.9) t 2 where E = the true local truncation error. For sufficiently small h, the higher-order terms t in Eq. (22.9) are usually negligible, and the result is often represented as  f (t , y ) i i 2 E = h (22.10) a 2 or 2 E = O(h ) (22.11) a where E = the approximate local truncation error. a According to Eq. (22.11), we see that the local error is proportional to the square of the step size and the first derivative of the differential equation. It can also be demon- strated that the global truncation error is O(h)—that is, it is proportional to the step sizecha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 559 22.2 EULER’S METHOD 559 (Carnahan et al., 1969). These observations lead to some useful conclusions: 1. The global error can be reduced by decreasing the step size. 2. The method will provide error-free predictions if the underlying function (i.e., the solution of the differential equation) is linear, because for a straight line the second derivative would be zero. This latter conclusion makes intuitive sense because Euler’s method uses straight-line seg- ments to approximate the solution. Hence, Euler’s method is referred to as a first-order method. It should also be noted that this general pattern holds for the higher-order one-step methods described in the following pages. That is, an nth-order method will yield perfect results if the underlying solution is an nth-order polynomial. Further, the local truncation n+1 n error will be O(h ) and the global error O(h ). 22.2.2 Stability of Euler’s Method In the preceding section, we learned that the truncation error of Euler’s method depends on the step size in a predictable way based on the Taylor series. This is an accuracy issue. The stability of a solution method is another important consideration that must be con- sidered when solving ODEs. A numerical solution is said to be unstable if errors grow exponentially for a problem for which there is a bounded solution. The stability of a par- ticular application can depend on three factors: the differential equation, the numerical method, and the step size. Insight into the step size required for stability can be examined by studying a very simple ODE: dy =−ay (22.12) dt If y(0) = y , calculus can be used to determine the solution as 0 −at y = y e 0 Thus, the solution starts at y and asymptotically approaches zero. 0 Now suppose that we use Euler’s method to solve the same problem numerically: dy i y = y + h i+1 i dt Substituting Eq. (22.12) gives y = y − ay h i+1 i i or y = y (1 − ah) (22.13) i+1 i The parenthetical quantity 1 − ah is called an amplification factor. If its absolute value is greater than unity, the solution will grow in an unbounded fashion. So clearly, the stability depends on the step size h. That is, if h 2/a, y→∞ as i →∞. Based on this analy- i sis, Euler’s method is said to be conditionally stable. Note that there are certain ODEs where errors always grow regardless of the method. Such ODEs are called ill-conditioned.cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 560 560 INITIAL-VALUE PROBLEMS Inaccuracy and instability are often confused. This is probably because (a) both repre- sent situations where the numerical solution breaks down and (b) both are affected by step size. However, they are distinct problems. For example, an inaccurate method can be very stable. We will return to the topic when we discuss stiff systems in Chap. 23. 22.2.3 MATLAB M-file Function: eulode We have already developed a simple M-file to implement Euler’s method for the falling bungee jumper problem in Chap. 3. Recall from Section 3.6, that this function used Euler’s method to compute the velocity after a given time of free fall. Now, let’s develop a more general, all-purpose algorithm. Figure 22.3 shows an M-file that uses Euler’s method to compute values of the dependent variable y over a range of values of the independent variable t. The name of the function holding the right-hand side of the differential equation is passed into the function as the FIGURE 22.3 An M-file to implement Euler’s method. function t,y = eulode(dydt,tspan,y0,h,varargin) % eulode: Euler ODE solver % t,y = eulode(dydt,tspan,y0,h,p1,p2,...): % uses Euler's method to integrate an ODE % input: % dydt = name of the M-file that evaluates the ODE % tspan = ti, tf where ti and tf = initial and % final values of independent variable % y0 = initial value of dependent variable % h = step size % p1,p2,... = additional parameters used by dydt % output: % t = vector of independent variable % y = vector of solution for dependent variable if nargin4,error('at least 4 input arguments required'),end ti = tspan(1);tf = tspan(2); if (tfti),error('upper limit must be greater than lower'),end t = (ti:h:tf)'; n = length(t); % if necessary, add an additional value of t % so that range goes from t = ti to tf if t(n)tf t(n+1) = tf; n = n+1; end y = y0ones(n,1); %preallocate y to improve efficiency for i = 1:n-1 %implement Euler's method y(i+1) = y(i) + dydt(t(i),y(i),varargin:)(t(i+1)-t(i)); endcha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 561 22.3 IMPROVEMENTS OF EULER’S METHOD 561 variable dydt. The initial and final values of the desired range of the independent variable is passed as a vector tspan. The initial value and the desired step size are passed as y0 and h, respectively. The function first generates a vector t over the desired range of the dependent variable using an increment of h. In the event that the step size is not evenly divisible into the range, the last value will fall short of the final value of the range. If this occurs, the final value is added to t so that the series spans the complete range. The length of the t vector is deter- mined as n. In addition, a vector of the dependent variable y is preallocated with n values of the initial condition to improve efficiency. At this point, Euler’s method (Eq. 22.5) is implemented by a simple loop: for i = 1:n-1 y(i+1) = y(i) + dydt(t(i),y(i),varargin:)(t(i+1)- t(i)); end Notice how a function is used to generate a value for the derivative at the appropriate val- ues of the independent and dependent variables. Also notice how the time step is automat- ically calculated based on the difference between adjacent values in the vector t. The ODE being solved can be set up in several ways. First, the differential equation can be defined as an anonymous function object. For example, for the ODE from Example 22.1: dydt=(t,y) 4exp(0.8t) - 0.5y; The solution can then be generated as t,y = eulode(dydt,0 4,2,1); disp(t,y) with the result (compare with Table 22.1): 0 2.0000 1.0000 5.0000 2.0000 11.4022 3.0000 25.5132 4.0000 56.8493 Although using an anonymous function is feasible for the present case, there will be more complex problems where the definition of the ODE requires several lines of code. In such instances, creating a separate M-file is the only option. 22.3 IMPROVEMENTS OF EULER’S METHOD A fundamental source of error in Euler’s method is that the derivative at the beginning of the interval is assumed to apply across the entire interval. Two simple modifications are available to help circumvent this shortcoming. As will be demonstrated in Section 22.4, both modifications (as well as Euler’s method itself) actually belong to a larger class of so- lution techniques called Runge-Kutta methods. However, because they have very straight- forward graphical interpretations, we will present them prior to their formal derivation as Runge-Kutta methods.cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 562 562 INITIAL-VALUE PROBLEMS 22.3.1 Heun’s Method One method to improve the estimate of the slope involves the determination of two deriv- atives for the interval—one at the beginning and another at the end. The two derivatives are then averaged to obtain an improved estimate of the slope for the entire interval. This ap- proach, called Heun’s method, is depicted graphically in Fig. 22.4. Recall that in Euler’s method, the slope at the beginning of an interval  y = f (t ,y ) (22.14) i i i is used to extrapolate linearly to y : i+1 0 y = y + f (t ,y )h (22.15) i i i i+1 For the standard Euler method we would stop at this point. However, in Heun’s method the 0 y calculated in Eq. (22.15) is not the final answer, but an intermediate prediction. This is i+1 why we have distinguished it with a superscript 0. Equation (22.15) is called a predictor equa- tion. It provides an estimate that allows the calculation of a slope at the end of the interval:   0 y = f t , y (22.16) i+1 i+1 i+1 Thus, the two slopes Eqs. (22.14) and (22.16) can be combined to obtain an average slope for the interval:  0 f (t , y ) + f t , y i i i+1  i+1 y ¯ = 2 This average slope is then used to extrapolate linearly from y to y using Euler’s i i+1 method:  0 f (t , y ) + f t , y i i i+1 i+1 y = y + h (22.17) i+1 i 2 which is called a corrector equation. FIGURE 22.4 Graphical depiction of Heun’s method. (a) Predictor and (b) corrector. y 0 y Slope  f(t , y ) i1 i1 0 f(t , y )  f(t , y ) i i i1 i1 Slope  Slope  f(t , y ) 2 i i t t t t t t i i1 i i1 (a) (b)cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 563 22.3 IMPROVEMENTS OF EULER’S METHOD 563 j1 m f (t , y )  f(t , y ) i i i1 i1 m j y  h y i i1 2 FIGURE 22.5 Graphical representation of iterating the corrector of Heun’s method to obtain an improved estimate. The Heun method is a predictor-corrector approach. As just derived, it can be ex- pressed concisely as 0 m Predictor (Fig. 22.4a): y = y + f (t ,y )h (22.18) i i i+1 i    j−1 m f t , y + f t , y i i+1 i i+1 j m Corrector (Fig. 22.4b): y = y + h (22.19) i+1 i 2 (for j = 1,2,..., m) Note that because Eq. (22.19) has y on both sides of the equal sign, it can be applied in i+1 an iterative fashion as indicated. That is, an old estimate can be used repeatedly to provide an improved estimate of y . The process is depicted in Fig. 22.5. i+1 As with similar iterative methods discussed in previous sections of the book, a termi- nation criterion for convergence of the corrector is provided by j j−1 y − y i+1 i+1 ε = × 100% a j y i+1 j−1 j where y and y are the result from the prior and the present iteration of the corrector, i+1 i+1 respectively. It should be understood that the iterative process does not necessarily con- verge on the true answer but will converge on an estimate with a finite truncation error, as demonstrated in the following example. EXAMPLE 22.2 Heun’s Method  0.8t Problem Statement. Use Heun’s method with iteration to integrate y = 4e − 0.5y from t = 0 to 4 with a step size of 1. The initial condition at t = 0 is y = 2. Employ a stop- ping criterion of 0.00001% to terminate the corrector iterations. Solution. First, the slope at (t , y ) is calculated as 0 0  0 y = 4e − 0.5(2) = 3 0 Then, the predictor is used to compute a value at 1.0: 0 y = 2 + 3(1) = 5 1cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 564 564 INITIAL-VALUE PROBLEMS 0.8t TABLE 22.2 Comparison of true and numerical values of the integral of y’ = 4e − 0.5y, with the initial condition that y = 2 at t = 0. The numerical values were computed using the Euler and Heun methods with a step size of 1. The Heun method was implemented both without and with iteration of the corrector. Without Iteration With Iteration ty y ε (%) y ε (%) y ε (%) true Euler t Heun t Heun t 0 2.00000 2.00000 2.00000 2.00000 1 6.19463 5.00000 19.28 6.70108 8.18 6.36087 2.68 2 14.84392 11.40216 23.19 16.31978 9.94 15.30224 3.09 3 33.67717 25.51321 24.24 37.19925 10.46 34.74328 3.17 4 75.33896 56.84931 24.54 83.33777 10.62 77.73510 3.18 Note that this is the result that would be obtained by the standard Euler method. The true value in Table 22.2 shows that it corresponds to a percent relative error of 19.28%. 0 Now, to improve the estimate for y , we use the value y to predict the slope at the i+1 1 end of the interval   0 0.8(1) y = f x , y = 4e − 0.5(5) = 6.402164 1 1 1 which can be combined with the initial slope to yield an average slope over the interval from t = 0 to 1: 3 + 6.402164  y ¯ = = 4.701082 2 This result can then be substituted into the corrector Eq. (22.19) to give the prediction at t = 1: 1 y = 2 + 4.701082(1) = 6.701082 1 which represents a true percent relative error of −8.18%. Thus, the Heun method without iteration of the corrector reduces the absolute value of the error by a factor of about 2.4 as compared with Euler’s method. At this point, we can also compute an approximate error as 6.701082 − 5 ε = × 100% = 25.39% a 6.701082 Now the estimate of y can be refined by substituting the new result back into the 1 right-hand side of Eq. (22.19) to give 0.8(1) 3 + 4e − 0.5(6.701082) 2 y = 2 + 1 = 6.275811 1 2 which represents a true percent relative error of 1.31 percent and an approximate error of 6.275811 − 6.701082 ε = × 100% = 6.776% a 6.275811cha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 565 22.3 IMPROVEMENTS OF EULER’S METHOD 565 The next iteration gives 0.8(1) 3 + 4e − 0.5(6.275811) 2 y = 2 + 1 = 6.382129 1 2 which represents a true error of 3.03% and an approximate error of 1.666%. The approximate error will keep dropping as the iterative process converges on a sta- ble final result. In this example, after 12 iterations the approximate error falls below the stopping criterion. At this point, the result at t = 1 is 6.36087, which represents a true rel- ative error of 2.68%. Table 22.2 shows results for the remainder of the computation along with results for Euler’s method and for the Heun method without iteration of the corrector. Insight into the local error of the Heun method can be gained by recognizing that it is related to the trapezoidal rule. In the previous example, the derivative is a function of both the dependent variable y and the independent variable t. For cases such as polynomials, where the ODE is solely a function of the independent variable, the predictor step Eq. (22.18) is not required and the corrector is applied only once for each iteration. For such cases, the technique is expressed concisely as f (t ) + f (t ) i i+1 y = y + h (22.20) i+1 i 2 Notice the similarity between the second term on the right-hand side of Eq. (22.20) and the trapezoidal rule Eq. (19.11). The connection between the two methods can be formally demonstrated by starting with the ordinary differential equation dy = f (t) (22.21) dt This equation can be solved for y by integration:   y t i+1 i+1 dy = f (t) dt (22.22) y t i i which yields  t i+1 y − y = f (t) dt (22.23) i+1 i t i or  t i+1 y = y + f (t) dt (22.24) i+1 i t i Now, recall that the trapezoidal rule Eq. (19.11) is defined as  t i+1 f (t ) + f (t ) i i+1 f (t) dt = h (22.25) 2 t icha01102_ch22_547-587.qxd 12/17/10 8:23 AM Page 566 566 INITIAL-VALUE PROBLEMS where h = t − t . Substituting Eq. (22.25) into Eq. (22.24) yields i+1 i f (t ) + f (t ) i i+1 y = y + h (22.26) i+1 i 2 which is equivalent to Eq. (22.20). For this reason, Heun’s method is sometimes referred to as the trapezoidal rule. Because Eq. (22.26) is a direct expression of the trapezoidal rule, the local truncation error is given by recall Eq. (19.14)  f (ξ) 3 E =− h (22.27) t 12 where ξ is between t and t . Thus, the method is second order because the second deriv- i i+1 ative of the ODE is zero when the true solution is a quadratic. In addition, the local and 3 2 global errors are O(h ) and O(h ), respectively. Therefore, decreasing the step size decreases the error at a faster rate than for Euler’s method. 22.3.2 The Midpoint Method Figure 22.6 illustrates another simple modification of Euler’s method. Called the midpoint method, this technique uses Euler’s method to predict a value of y at the midpoint of the interval (Fig. 22.6a): h y = y + f (t ,y ) (22.28) i+1/2 i i i 2 Then, this predicted value is used to calculate a slope at the midpoint:  y = f (t ,y ) (22.29) i+1/2 i+1/2 i+1/2 FIGURE 22.6 Graphical depiction of midpoint method. (a) Predictor and (b) corrector. y y Slope  f(t , y ) i12 i12 Slope  f(t , y ) i12 i12 t t t t t t t i i12 i1 i i1 (b) (a)

Advise: Why You Wasting Money in Costly SEO Tools, Use World's Best Free SEO Tool Ubersuggest.