Report

ORDINARY DIFFERENTIAL EQUATIONS ENGR 351 Numerical Methods for Engineers Southern Illinois University Carbondale College of Engineering Dr. L.R. Chevalier Dr. B.A. DeVantier Ordinary Differential Equations… where to use them The dissolution (solubilization) of a contaminant into groundwater is governed by the equation: dC kl C s C dt where kl is a lumped mass transfer coefficient and Cs is the maximum solubility of the contaminant into the water (a constant). Given C(0)=2 mg/L, Cs = 500 mg/L and kl = 0.1 day-1, estimate C(0.5) and C(1.0) using a numerical method for ODE’s. Ordinary Differential Equations… where to use them A mass balance for a chemical in a completely mixed reactor can be written as: dc V F Qc kVc 2 dt where V is the volume (10 m3), c is concentration (g/m3), F is the feed rate (200 g/min), Q is the flow rate (1 m3/min), and k is reaction rate (0.1 m3/g/min). If c(0)=0, solve the ODE for c(0.5) and c(1.0) Ordinary Differential Equations… where to use them Before coming to an exam Friday afternoon, Mr. Bringer forgot to place 24 cans of a refreshing beverage in the refrigerator. His guest are arriving in 5 minutes. So, of course he puts the beverage in the refrigerator immediately. The cans are initially at 75, and the refrigerator is at a constant temperature of 40. Ordinary Differential Equations… where to use them The rate of cooling is proportional to the difference in the temperature between the beverage and the surrounding air, as expressed by the following equation with k = 0.1/min. dT k T Tair dt Use a numerical method to determine the temperature of the beverage after 5 minutes and 10 minutes. Ordinary Differential Equations • A differential equation defines a relationship between an unknown function and one or more of its derivatives • Physical problems using differential equations • electrical circuits • heat transfer • motion Ordinary Differential Equations • The derivatives are of the dependent variable with respect to the independent variable • First order differential equation with y as the dependent variable and x as the independent variable would be: • dy/dx = f(x,y) Ordinary Differential Equations • A second order differential equation would have the form: d2y dy f x, y, 2 dx dx } does not necessarily have to include all of these variables Ordinary Differential Equations • An ordinary differential equation is one with a single independent variable. • Thus, the previous two equations are ordinary differential equations • The following is not: dy f x1 , x2 , y dx1 Ordinary Differential Equations • The analytical solution of ordinary differential equation as well as partial differential equations is called the “closed form solution” • This solution requires that the constants of integration be evaluated using prescribed values of the independent variable(s). Ordinary Differential Equations • An ordinary differential equation of order n requires that n conditions be specified. • Boundary conditions • Initial conditions Ordinary Differential Equations • An ordinary differential equation of order n requires that n conditions be specified. • Boundary conditions • Initial conditions consider this beam where the deflection is zero at the boundaries x= 0 and x = L These are boundary conditions consider this beam where the deflection is zero at the boundaries x= 0 and x = L These are boundary conditions P a yo In some cases, the specific behavior of a system(s) is known at a particular time. Consider how the deflection of a beam at x = a is shown at time t =0 to be equal to yo. Being interested in the response for t > 0, this is called the initial condition. Ordinary Differential Equations • At best, only a few differential equations can be solved analytically in a closed form. • Solutions of most practical engineering problems involving differential equations require the use of numerical methods. Scope of Lectures on ODE • One Step Methods • • • • • Euler’s Method Heun’s Method Improved Polygon Runge Kutta Systems of ODE • Adaptive step size control Scope of Lectures on ODE • Boundary value problems • Case studies Specific Study Objectives • Understand the visual representation of Euler’s, Heun’s and the improved polygon methods. • Understand the difference between local and global truncation errors • Know the general form of the Runge-Kutta methods. • Understand the derivation of the secondorder RK method and how it relates to the Taylor series expansion. Specific Study Objectives • Realize that there are an infinite number of possible versions for second- and higherorder RK methods • Know how to apply any of the RK methods to systems of equations • Understand the difference between initial value and boundary value problems Review of Analytical Solution dy 4x2 dx 2 dy 4 x dx 3 4x y C 3 At this point lets consider initial conditions. y(0)=1 and y(0)=2 4x3 y C 3 for y0 1 40 1 C 3 then C 1 3 for y0 2 40 2 C 3 and C 2 3 What we see are different values of C for the two different initial conditions. The resulting equations are: 4x3 y 1 3 4x3 y 2 3 16 y (0 )=1 y (0 )=2 12 y (0 )=3 y y (0 )=4 8 4 0 0 0 .5 1 x 1 .5 2 2 .5 One Step Methods • Focus is on solving ODE in the form dy f x , y dx yi 1 yi h h y yi+1 yi slope = x This is the same as saying: new value = old value + slope x step size Euler’s Method • The first derivative provides a direct estimate of the slope at xi • The equation is applied iteratively, or one step at a time, over small distance in order to reduce the error • Hence this is often referred to as Euler’s One-Step Method Example For the initial condition y(1)=1, determine y for h = 0.1 analytically and using Euler’s method given: dy 2 4x dx Error Analysis of Euler’s Method • Truncation error - caused by the nature of the techniques employed to approximate values of y • local truncation error (from Taylor Series) • propagated truncation error • sum of the two = global truncation error • Round off error - caused by the limited number of significant digits that can be retained by a computer or calculator Ex a m ple 12 A n a ly t ica l Solu t ion 10 y 8 Nu m er ica l Solu t ion 6 4 2 0 0 0 .5 1 1 .5 2 2 .5 x ....end of example Higher Order Taylor Series Methods yi 1 f ' xi , yi 2 yi f xi , yi h h 2 • This is simple enough to implement with polynomials • Not so trivial with more complicated ODE • In particular, ODE that are functions of both dependent and independent variables require chain-rule differentiation • Alternative one-step methods are needed Modification of Euler’s Methods • A fundamental 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 will be demonstrated • These modification actually belong to a larger class of solution techniques called Runge-Kutta which we will explore later. Heun’s Method • Determine the derivative for the interval • the initial point • end point • Use the average to obtain an improved estimate of the slope for the entire interval y Use this “average” slope to predict yi+1 xi xi+1 { yi 1 yi f xi , yi f xi 1 , yi 1 h 2 y f xi , yi f xi 1 , yi 1 yi 1 yi h 2 y xi xi+1 xi xi+1 x f xi , yi f xi 1 , yi 1 yi 1 yi h 2 y yi 1 yi h xi xi+1 x Improved Polygon Method • Another modification of Euler’s Method • Uses Euler’s to predict a value of y at the midpoint of the interval • This predicted value is used to estimate the slope at the midpoint Improved Polygon Method • We then assume that this slope represents a valid approximation of the average slope for the entire interval • Use this slope to extrapolate linearly from xi to xi+1 using Euler’s algorithm Runge-Kutta Methods Both Heun’s and the Improved Polygon Method have been introduced graphically. However, the algorithms used are not as straight forward as they can be. Let’s review the Runge-Kutta Methods. Choices in values of variable will give us these methods and more. It is recommend that you use this algorithm on your homework and/or programming assignments. Runge-Kutta Methods • RK methods achieve the accuracy of a Taylor series approach without requiring the calculation of a higher derivative • Many variations exist but all can be cast in the generalized form: { yi 1 yi xi , yi , hh is called the incremental function , Incremental Function can be interpreted as a representative slope over the interval a1k1 a2 k 2 an k n where the a ' s are constant and the k ' s are: k1 f xi , yi k 2 f xi p1h , yi q11k1h k 3 f xi p2 h , yi q21k1h q22 k 2 h k n f xi pn h , yi qn 1,1k1h qn 1, 2 k 2 h qn 1, n 1k n 1h a1k1 a2 k 2 an k n where the a ' s are constant and the k ' s are: k1 f xi , yi k 2 f xi p1h , yi q11k1h k 3 f xi p2 h , yi q21k1h q22 k 2 h k n f xi pn h , yi qn 1,1k1h qn 1, 2 k 2 h qn 1, n 1k n 1h NOTE: k’s are recurrence relationships, that is k1 appears in the equation for k2 which appears in the equation for k3 This recurrence makes RK methods efficient for computer calculations Second Order RK Methods yi 1 yi a1k1 a2 k 2 h where k1 f xi , yi k 2 f xi p1h, yi q11k1h a1k1 a2 k 2 an k n where the a ' s are constant and the k ' s are: k1 f xi , yi k 2 f xi p1h , yi q11k1h k 3 f xi p2 h , yi q21k1h q22 k 2 h k n f xi pn h , yi qn 1,1k1h qn 1, 2 k 2 h qn 1, n 1k n 1h Second Order RK Methods • We have to determine values for the constants a1, a2, p1 and q11 • To do this consider the Taylor series in terms of yi+1 and f(xi,yi) yi 1 yi a1k1 a2 k 2 h h2 yi 1 yi f xi , yi h f ' xi , yi 2 Now, f’(xi , yi ) must be determined by the chain rule for differentiation f f dy f ' xi , yi x y dx substituting f f dy h 2 yi 1 yi f xi , yi h x y dx 2 The basic strategy underlying Runge-Kutta methods is to use algebraic manipulations to solve for values of a1, a2, p1 and q11 yi 1 yi a1k1 a2 k2 h f f dy h 2 yi 1 yi f xi , yi h x y dx 2 By setting these two equations equal to each other and recalling: k1 f xi , yi k2 f xi p1h, yi q11k1h we derive three equations to evaluate the four unknown constants a1 a2 1 1 2 1 a2 q11 2 a2 p1 Because we have three equations with four unknowns, we must assume a value of one of the unknowns. Suppose we specify a value for a2. What would the equations be? a1 1 a2 1 p1 q11 2a2 Because we can choose an infinite number of values for a2 there are an infinite number of second order RK methods. Every solution would yield exactly the same result if the solution to the ODE were quadratic, linear or a constant. Lets review three of the most commonly used and preferred versions. y i 1 y i a 1 k 1 a 2 k 2 h Consider the following: where k 1 f xi , yi k 2 f xi p1h, yi q11 k1h a1 a 2 1 1 a 2 p1 2 1 a 2 q11 2 Case 1: a2 = 1/2 Case 2: a2 = 1 These two methods have been previously studied. What are they? Case 1: a2 = 1/2 a1 1 a2 1 1 / 2 1 / 2 1 a2 p1 2 1 a2 q11 2 1 p1 q11 1 2a2 1 1 yi 1 yi k1 k 2 h 2 2 where k1 f xi , yi k 2 f xi h, yi k1h This is Heun’s Method with a single corrector. Note that k1 is the slope at the beginning of the interval and k2 is the slope at the end of the interval. yi 1 yi a1k1 a2 k 2 h where k1 f xi , yi k 2 f xi p1h, yi q11k1h a1 1 a2 1 1 0 Case 2: a2 = 1 1 2 1 a2 q11 2 This is the Improved Polygon Method. a2 p1 1 1 p1 q11 2 a2 2 yi 1 yi k 2 h where k1 f xi , yi 1 1 k 2 f xi h, yi k1h 2 2 yi 1 yi a1k1 a2 k 2 h where k1 f xi , yi k 2 f xi p1h, yi q11k1h Ralston’s Method Ralston (1962) and Ralston and Rabinowitiz (1978) determined that choosing a2 = 2/3 provides a minimum bound on the truncation error for the second order RK algorithms. This results in a1 = 1/3 and p1 = q11 = 3/4 1 2 yi 1 yi k1 k 2 h 3 3 where k1 f xi , yi 3 3 k 2 f xi h, yi k1h 4 4 Example dy 4x 2 y dx I .C. y 1 at x 1 i.e. y 1 1 step size h 0.1 Evaluate the following ODE using Heun’s Methods Third Order Runge-Kutta Methods • Derivation is similar to the one for the second-order • Results in six equations and eight unknowns. • One common version results in the following 1 yi 1 yi k 1 4 k 2 k 3 h 6 where k 1 f xi , yi Note the third term 1 1 k 2 f xi h, yi k1h 2 2 k 3 f xi h, yi hk1 2hk 2 NOTE: if the derivative is a function of x only, this reduces to Simpson’s 1/3 Rule Fourth Order Runge Kutta • The most popular • The following is sometimes called the classical fourth-order RK method 1 yi 1 yi k 1 2 k 2 2 k 3 k 4 h 6 where k 1 f xi , yi 1 1 k 2 f xi h, yi k1h 2 2 1 1 k 3 f xi h, yi hk 2 2 2 k 4 f xi h, yi hk 3 • Note that for ODE that are a function of x alone that this is also the equivalent of Simpson’s 1/3 Rule 1 yi 1 yi k1 2k 2 2k3 k 4 h 6 where k1 f xi , yi 1 1 k 2 f xi h, yi k1h 2 2 1 1 k3 f xi h, yi hk2 2 2 k 4 f xi h, yi hk3 Example Use 4th Order RK to solve the following differential equation: dy xy dx 1 x 2 I . C. y1 1 using an interval of h = 0.1 Higher Order RK Methods • When more accurate results are required, Bucher’s (1964) fifth order RK method is recommended • There is a similarity to Boole’s Rule • The gain in accuracy is offset by added computational effort and complexity Systems of Equations • Many practical problems in engineering and science require the solution of a system of simultaneous differential equations dy1 f 1 x , y1 , y2 , , yn dx dy2 f 2 x , y1 , y2 , , yn dx dyn f n x , y1 , y2 , , yn dx • Solution requires n initial conditions • All the methods for single equations can be used • The procedure involves applying the one-step technique for every equation at each step before proceeding to the next step dy1 f 1 x , y1 , y 2 , , y n dx dy2 f 2 x , y1 , y2 , , yn dx dyn f n x , y1 , y2 , , yn dx Boundary Value Problems • Recall that the solution to an nth order ODE requires n conditions • If all the conditions are specified at the same value of the independent variable, then we are dealing with an initial value problem • Problems so far have been devoted to this type of problem Boundary Value Problems • In contrast, we may also have conditions a different value of the independent variable. • These are often specified at the extreme point or boundaries of as system and customarily referred to as boundary value problems • To approaches to the solution • shooting method • finite difference approach General Methods for Boundary Value Problems The conservation of heat can be used to develop a heat balance for a long, thin rod. If the rod is not insulated along its length and the system is at steady state. The equation that results is: d 2T h ' Ta T 0 2 dx Ta T1 T2 Ta Ta T1 d 2T h ' Ta T 0 2 dx Clearly this second order ODE needs 2 conditions. This can be satisfied by knowing the temperature at the boundaries, i.e. T1 and T2 T2 Ta T(0) = T1 T(L) = T2 d 2T h ' Ta T 0 2 dx T(0) = T1 T(L) = T2 Use these conditions to solve the equation analytically. For a 10 m rod with Ta = 20 T(0) = 40 T(10) = 200 h’ = 0.01 T 73.45e0.1x 53.45e0.1x 20 Now that we have an analytical solution, lets evaluate our two proposed numerical methods. Shooting Method Given: d 2T h ' Ta T 0 2 dx dT z dx dz h ' Ta T dx We need an initial value of z. For the shooting method, guess an initial value. Guessing z(0) = 10 dz h ' Ta T dx Guessing z(0) = 10 Using a fourth-order RK method with a step size of 2, T(10) = 168.38 This differs from the BC T(10) = 200 Making another guess, z(0) = 20 T(10) = 285.90 Because the original ODE is linear, the estimates of z(0) are linearly related. Using a linear interpolation formula between the values of z(0), determine a new value of z(0) Recall: first estimate z(0) = 10 T(20) = 168.38 second estimate z(0)=20 T(20) = 285.90 What is z(0) that would give us T(20)=200? T(20) 300 250 200 150 0 5 10 15 z(0) 20 25 T(20) 300 250 200 150 0 5 10 15 20 25 z(0) 20 10 z 0 10 200 168.38 12.69 285.90 168.38 We can now use this to solve the first order ODE d 2T h ' Ta T 0 2 dx dT z dx dz h ' Ta T dx 250 150 Sh oot in g Met h od T 2 00 A n a ly t ica l Solu t ion 1 00 50 0 0 5 10 dist a n ce (m ) For nonlinear boundary value problems, linear interpolation will not necessarily result in an accurate estimation. One alternative is to apply three applications of the shooting method and use quadratic interpolation.. Finite Difference Methods The finite divided difference approximation for the 2nd derivative can be substituted into the governing equation. d 2 T Ti 1 2 Ti Ti 1 2 dx x 2 d 2T h ' Ta T 0 2 dx Ti 1 2 Ti Ti 1 h ' Ta Ti 0 2 x Collect terms Ti 1 2 Ti Ti 1 h ' Ta Ti 0 2 x Ti 1 2 h ' x 2 Ti Ti 1 h ' x 2 Ta We can now apply this equation to each interior node on the rod. Divide the rod into a grid, and consider a “node” to be at each division. i.e.. x = 2m x=2m T1 T2 L = 10 m Ti 1 2 h' x 2 Ti Ti 1 h' x 2Ta x=2m T(0) T(10) L = 10 m Consider the previous problem: L = 10 m We need to solve for the Ta = 20 temperature at the interior T(0) = 40 T(10) = 200 nodes (4 unknowns). h’ = 0.01 Apply the governing equation at these nodes (4 equations). What is the matrix? Ti 1 2 h' x 2 Ti Ti 1 h' x 2Ta x=0 2 4 6 8 10 T(0) T(10) i=0 1 2 3 4 5 Notice the labeling for numbering x and i Ti 1 2 h' x 2 Ti Ti 1 h' x 2Ta x=0 2 4 6 8 10 T(0) T(10) i=0 40 1 2 3 4 5 200 Note also that the dependent values are known at the boundaries (hence the term boundary value problem) Ti 1 2 h' x 2 Ti Ti 1 h' x 2Ta x=0 2 4 6 8 10 T(0) T(10) i=0 40 1 2 3 4 5 200 Apply the governing equation at node 1 T0 2 h' x 2 T1 T2 h' x 2Ta 40 2.04T1 T2 0.8 2.04T1 T2 40.8 Ti 1 2 h' x 2 Ti Ti 1 h' x 2Ta x=0 2 4 6 8 10 T(0) T(10) i=0 40 1 2 3 4 Apply the equation at node 2 T1 2 h' x 2 T2 T3 h' x 2Ta T1 2.04T2 T3 0.8 5 200 Ti 1 2 h' x 2 Ti Ti 1 h' x 2Ta x=0 2 4 6 8 10 T(0) T(10) i=0 1 2 3 4 5 40 200 We get a similar equation at node 3 T2 2 h' x 2 T3 T4 h' x 2Ta T2 2.04T3 T4 0.8 Ti 1 2 h' x 2 Ti Ti 1 h' x 2Ta x=0 2 4 6 8 10 T(0) T(10) i=0 40 1 2 3 4 At node 4, we consider the boundary at the right. T3 2 h' x 2 T4 T5 h' x 2Ta T3 2.04T4 200 0.8 T3 2.04T4 200.8 5 200 For the four interior nodes, we get the following 4 x 4 matrix 0 0 T1 40.8 2.04 1 1 2.04 1 0 T2 0.8 T 0 1 2 . 04 1 0 . 8 3 0 T 200.8 0 1 2 . 04 4 T T 65.97 93.78 124.54 159.48 250 A n a ly t ica l Solu t ion 2 00 Sh oot in g Met h od Fin it e Differ en ce T 150 1 00 50 0 0 5 dist a n ce (m ) 10 Example Consider the previous example, but with x=1. What is the matrix? Specific Study Objectives • Understand the visual representation of Euler’s, Heun’s and the improved polygon methods. • Understand the difference between local and global truncation errors • Know the general form of the Runge-Kutta methods. • Understand the derivation of the secondorder RK method and how it relates to the Taylor series expansion. Specific Study Objectives • Realize that there are an infinite number of possible versions for second- and higherorder RK methods • Know how to apply any of the RK methods to systems of equations • Understand the difference between initial value and boundary value problems … end of lecture on ODE