Docstoc

Ordinary Differential Equations

Document Sample
Ordinary Differential Equations Powered By Docstoc
					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:

                   kl Cs  C 
               dC
               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  kVc2
             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.


            k T  Tair 
        dT
        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, 
 dx 2           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 second-
  order 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 higher-
   order 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                  At this point lets consider
      4x2
  dx                  initial conditions.
   dy   4 x 2 dx
                      y(0)=1
          3
     4x
  y    C            and
      3
                      y(0)=2
    4x3
y      C     What we see are different
     3
for y0  1   values of C for the two
               different initial conditions.
   40
       3

1       C
     3
               The resulting equations
then C  1
               are:
for y0  2
   40                  4x3
       3

2      C            y     1
     3                    3
and C  2             y
                         4x3
                             2
                          3
    16

                    y (0 )=1
                    y (0 )=2
    12
                    y (0 )=3
                    y (0 )=4

     8
y




     4



     0
         0   0 .5         1        1 .5   2   2 .5
                               x
One Step Methods

• Focus is on solving ODE in the form

   dy                             h
         f  x , y
                       y

   dx
                                        yi+1
   yi 1  yi  h
                             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
                  4x 2

              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
    10
                    Solu t ion
     8              Nu m er ica l
                    Solu t ion
     6
y




     4

     2

     0
         0   0 .5         1            1 .5   2         2 .5
                                 x

                                                  ....end of example
Higher Order Taylor Series
Methods
                                  f '  xi , yi  2
yi  1    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




                                                            x
                       xi             xi+1
             f  xi , yi   f xi 1 , yi 1 
yi 1  yi                                     h
                            2



                                               y


     y i 1  y i  h




                                                                x
                                                    xi   xi+1
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 , h h
                          {
             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 k 2 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 
k 2  f  xi  p1h, yi  q11k1h

we derive three equations to evaluate the four unknown
constants
           a1  a2  1
                    1
           a2 p1 
                    2
                    1
           a2 q11 
                    2

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                       Case 1: a2 = 1/2
k 2  f  xi  p1h, yi  q11 k1h        Case 2: a2 = 1
a1  a 2  1
          1                              These two methods
a 2 p1                                  have been previously
          2
          1                              studied.
a 2 q11 
          2
                                         What are they?
                                  Case 1: a2 = 1/2
a1  1  a2  1  1 / 2  1 / 2
                                  This is Heun’s Method with
         1
a2 p1                            a single corrector.
         2
         1                        Note that k1 is the slope at
a2 q11 
         2                        the beginning of the interval
            1                     and k2 is the slope at the
p1  q11       1
           2a2                    end of the interval.
              1    1 
yi  1  yi   k1  k 2  h       yi  1  yi  a1k1  a2 k 2 h
              2    2 
                                   where
where
                                   k1  f  xi , yi 
k1  f  xi , yi 
                                   k 2  f  xi  p1h, yi  q11k1h
k 2  f  xi  h, yi  k1h
a1  1  a2  1  1  0       Case 2: a2 = 1
         1
a2 p1                        This is the Improved
         2
                              Polygon Method.
         1
a2 q11 
         2
               1      1
p1  q11           
             2a2 2            yi  1  yi  a1k1  a2 k 2 h
yi  1  yi  k 2 h           where
where                         k1  f  xi , yi 
k1  f  xi , yi 
                              k 2  f  xi  p1h, yi  q11k1h
             1       1 
k 2  f  xi  h, yi  k1h
             2       2 
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 k  2 k h
                             yi 1  yi   1       2
                                          3     3 
                            where
                            k1  f  xi , yi 
                                    x  3 h, y  3 k h
                            k2  f  i               1 
                                        4
                                               i
                                                  4 
Example
dy                                      Evaluate the following
     4x 2 y
dx                                      ODE using Heun’s
I .C. y  1 at x  1 i.e. y 1  1   Methods
step size h  0.1
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 

             k3  f  xi  h, yi  hk2 
                            1         1
                                         
                           2         2   
             k 4  f  xi  h, yi  hk3 
Example

Use 4th Order RK to solve the following differential equation:


     dy   xy
                  I . C. y1  1
     dx 1  x 2


 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 , y2 , , yn 
        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
   2
      h '  Ta  T   0
dx
                                  Ta
                     T1                            T2
                                  Ta
                            Ta
                T1                      T2

d 2T                        Ta
   2
      h '  Ta  T   0
dx
Clearly this second order   T(0) = T1
ODE needs 2 conditions.
This can be satisfied by    T(L) = T2
knowing the temperature
at the boundaries,
i.e. T1 and T2
d 2T
   2
      h '  Ta  T   0          Use these conditions to solve
dx
                                   the equation analytically.
  T(0) = T1
                                   For a 10 m rod with
  T(L) = T2                        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
    2
       h '  Ta  T   0
 dx
 dT
     z                      We need an initial value
 dx                          of z.
 dz
     h '  Ta  T 
 dx                          For the shooting method, guess
                             an initial value.

                             Guessing z(0) = 10
dz
    h '  Ta  T     Guessing z(0) = 10
dx

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?

                  300

                  250
          T(20)




                  200

                  150
                        0   5   10          15   20   25
                                     z(0)
          300

          250
  T(20)



          200

          150
                0   5   10          15   20   25
                             z(0)


                   20  10
  z 0  10                   200  168.38  12.69
               285.90  168.38
                                                   d 2T
                                                      2
                                                         h '  Ta  T   0
                                                   dx
We can now use                                     dT
this to solve the first                                z
                                                   dx
order ODE                                          dz
                                                       h '  Ta  T 
                                                   dx
     250

                A n a ly t ica l
     2 00       Solu t ion


                Sh oot in g
     150        Met h od
 T




     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
             
      dx  2
                       x 2
      d 2T
          2
              h ' Ta  T   0
      dx
      Ti 1  2 Ti  Ti 1
                            h ' Ta  Ti   0
              x 2
Collect terms
                Ti 1  2 Ti  Ti 1
                                      h ' Ta  Ti   0
                        x  2


                 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
Ta = 20            We need to solve for the
T(0) = 40          temperature at the interior
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      1       2         3        4   5




40                                                           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      1       2         3        4   5




40                                             200
          Apply the governing equation at node 1

           T0  2  h' x 2 T1  T2  h' x 2Ta
           40  2.04T1  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      1       2         3        4   5




40                                                           200
              Apply the equation at node 2


           T1  2  h' x 2 T2  T3  h' x 2Ta
           T1  2.04T2  T3  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      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      1       2         3        4   5



              At node 4, we consider the
40            boundary at the right.                         200

           T3  2  h' x 2 T4  T5  h' x 2Ta
           T3  2.04T4  200  0.8
           T3  2.04T4  200.8
For the four interior nodes, we get the following
4 x 4 matrix


      2.04  1   0    0  T1   40.8 
        1 2.04  1   0  T2   0.8 
                            
                              
                                        
                                        
       0     1 2.04  1  T3   0.8 
       0          1 2.04 T4  200.8
            0                      



     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
    150        Fin it e
               Differ en ce
T




    1 00


     50


      0
           0                 5          10
                     dist a n ce (m )
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 second-
  order 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 higher-
   order 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

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:9
posted:9/28/2012
language:English
pages:81