Docstoc

Basics of Ordinary Differential Equations

Document Sample
Basics of Ordinary Differential Equations Powered By Docstoc
					Introduction to Numerical Solutions
 of Ordinary Differential Equations

              Larry Caretto
      Mechanical Engineering 501AB
   Seminar in Engineering Analysis

          October 29, 2003
                  Outline
• Review last class and homework
• Midterm Exam November 4 covers
  material up to and including last week’s
  lecture
• Overivew of numerical solutions
  – Initial value problems in first-order
    equations
  – Systems of first order equations and initial
    value problems in higher order equations
  – Boundary value problems
  – Stiff systems and eigenvalues
                                                   2
        Review Last Class
• Systems of equations
  – Combine to higher order equation
  – Matrix approaches
  – Reduction of order
• Laplace transforms for ODEs
  – Transform ODE from y(t) to Y(s)
  – Solve for Y(s)
  – Rearrange solution
  – Transform back to y(t)

                                       3
      Simultaneous Solution
• Differentiate first (y1) equation to obtain
  derivative(s) of y2 that are present in
  second equation
• Solve first equation for y2
• Substitute equations for y2 and its deriv-
  atives into the result of the first step
• Result is a higher order equation that
  can be solved (if possible) by usual
  methods

                                            4
   Matrix Differential Equations II
                         dy
  • Matrix components in     Ay  r
                         dt
   y1                  dy1 
  y             y1   dt          a11   a12    a13   a1n 
                  y  dy2          a             a23   a2 n 
   2                                 21    a22                 
                  2       dt 
   y3    dy d  y3   dy3          a31   a32    a33   a3n 
y                    dt    A                           
        dt dt    
                                                         
                                                  
                                                          
                  yn  dyn 
                                    an1
                                             an 2   an 3   ann 
                                                                  
   yn 
                          dt 

                 r  r1 r2         r3   rn 
                                                       T

                                                              5
                    dy
       Solving          Ay  r
                    dt
• Assume that A(n x n) has n linearly
  independent eigenvectors
• Eigenvectors are columns of matrix, X
• Define a new vector s = X-1y (y = Xs)
• Substitute y = Xs into the matrix
  differential equation
• Obtain equation for s using L = X-1AX
  that gives independent equations
                                          6
                            ds
        Terms in Solution of  Ls  p
                            dt
          l1t
        e                                                         p1 
                    0        0          0            C1       l1 
                                                     C 
         0      e  l2 t    0          0                       p2 
                                                         2
                                                                   l 
         0                 e l3t                   C3 
E(t )  
                    0                     0
                                                  C    1      p 2
                                                       L p 3l 
                                                                   3
                                                    
        
                                 
                                                                   
         0                              ln t
                                         e                       
                   0        0                        Cn 
                                                                p    
                                                                   nl 
  • With these definitions, i            s = C elit + p /l            n
                              i                          i   i
    becomes s = EC + L-1p (E(0) = I)
                                                                   7
Apply Initial Conditions on y(0)
• Get y = Xs = XEC + XL-1p
• At t = 0, E = I, and y = y0 = initial y
  components, giving y0 = XIC + XL-1p
• Premultiply by X-1 to obtain X-1y0 =
  X-1XC + X-1XL-1p = C + L-1p
• Constant vector, C = X-1y0 - L-1p
• Result: y = XE [X-1y0 - L-1p] + XL-1p
• Homogenous(p = 0): y = XEX-1y0
                                       8
        Reduction of Order
• An nth order equation can be written as
  a system of n first order equations
• We can write a general nonlinear nth
  order equation as shown below
           dny            dy  d n1 y 
                 f  x, y, ,, n1 
                    
           dx n
                          dx  dx     
• Redefine y as z0 and define the
  following derivatives
         dy dz0         d 2 y d 2 z0 dz1
                z1        2
                                  2
                                          z2
         dx dx          dx      dx     dx        9
          Reduction of Order II
 • Continue this definition up to zn-1
  d n y d n z0 d n 1 z1 d n  2 z 2     dzn  2 dzn 1
      n
            n
                n 1        n2
                                          2
                                                
  dx      dx    dx       dx               dx      dx
 • z0 to zn-1 are n variables that satisfy
   simultaneous, first-order ODEs
dy n               d n1 y  dzn1
                             dx n  f x, z0 , z1 ,, zn 1 
                dy
   n
      f  x, y, ,, n1  
         
dx             dx  dx 

                     dzk
                          zk 1 k  0,1,, n  2
                     dx                                     10
      Reduction of Order III
• The main application of reduction of
  order is to numerical methods
• Numerical solutions of ODEs develop
  methods to solve a system of first order
  equations
• Higher order equations are solved by
  converting them to a system of first
  order equations

                                         11
   Simple Laplace Transforms
f(t)       F(s)           f(t)         F(s)
ctn        cn!/sn+1       ceatsin wt         cw
ctx        cG(x+1)/sx+1                ( s  a) 2  w 2
ceat       c/(s – a)      ceatcos wt       c( s  a)
csin wt    cw/(s2 + w2)                ( s  a) 2  w 2
ccos wt    cs/(s2 + w2)   Additional transforms
csinh wt   cw/(s2 - w2)   in Table 5.9, pp 297-
                          299 of Kreyszig
ccosh wt   cs/(s2 - w2)
                                                     12
             Transforms of Derivatives
      • The main formulas for differential
        equations are the Laplace transforms of
        derivatives
      • These have transform of desired
        function, L [f(t)] = F(s) and the initial
        conditions, f(0), f’(0), etc.
                               [ f ' (t )]  sF (s)  f (0)

                       [ f ' ' (t )]  s 2 F ( s)  sf (0)  f ' (0)

[ f ( n ) (t )]  s n F ( s)  s n 1 f (0)  s n  2 f ' (0)    sf ( n  2) (0)  f ( n 1) (0)
                                                                                               13
Solving Differential Equations
• Transform all terms in the differential
  equation to get an algebraic equation
  – For a differential equation in y(t) we get the
    transforms Y(s) = L [y(t)]
  – Similar notation for other transformed
    functions in the equation R(s) = L [r(t)]
• Solve the algebraic equation for Y(s)
• Obtain the inverse transform for Y(s)
  from tables to get y(t)
                                                 14
          Partial Fractions
• Method to convert fraction with several
  factors in denominator into sum of
  individual factors (in denominator)
• Necessary to modify complicated
  solution for Y(s) into simpler functions
  whose transforms are known
• Special rules for repeated factors and
  for complex conjugates

                                             15
         Shifting Theorems
• Sometimes required for transform table
  – First theorem is defined for function f(t)
    with transform F(s)
  – F(s - a) has inverse of e-atf(t)
  – Example Y(s) = [2(s+1)+4]/[(s + 1)2 + 1]
  – For F(s) = s/(s2 + w2); f(t) = cos wt and f(t) =
    cos wt for F(s) = w/(s2 + w2)
  – Here Y(s+1) = 2(s+1)+4]/[(s + 1)2 + 1] +
    4/[(s + 1)2 + 1] so y(t) = e-t times Y(s)
    inverse
  – y(t) = e-t[2 cos t + 4 sin t]

                                                  16
       Shifting Theorems II
• Second shifting theorem
  – Applies to e-asF(s), where F(s) is known
    transform of a function f(t)
  – Inverse transform is f(t – a) u(t – a) where
    u(t – a) is the unit step function
  – If we have e-ass/(s2 + w2), we recognize
    s/(s2 + w2) as F(s) for f(t) = cos wt
  – Thus e-ass/(s2 + w2) is the Laplace
    transform for cos[wt(t – a)] u(t – a)

                                                   17
Numerical Analysis Problems
• Numerical solution of algebraic
  equations and eigenvalue problems
• Solution of one or more nonlinear
  algebraic equations f(x) = 0
• Linear and nonlinear optimization
• Constructing interpolating polynomials
• Numerical quadrature
• Numerical differentiation
• Numerical differential equations
                                           18
            Interpolation
• Start with N data pairs xi, yi
• Find a function (polynomial) that can be
  used for interpolation
• Basic rule: the interpolation polynomial
  must fit all points exactly
• Denote the polynomial as p(x)
• The basic rule is that p(xi) = yi
• Many different forms
                                         19
       Newton Polynomials
• p(x) = a0 + a1(x – x0) + a2(x – x0)(x – x1)
  + a3(x – x0)(x – x1)(x – x2) + … + an-1(x –
  x0)(x – x1)(x – x2) … (x – xn-2)
• Terms with factors of x – xi are zero
  when x = xi
  – Use this and rule that p(xi) = yi to find ai
• a0 = y0, a1 = (y1 – y0) / (x1 – x0)
• y2 = a0 + a1(x2 – x0) + a2(x2 – x0)(x2 – x1)
  – Solve for a2 using results for a0 and a1

                                                   20
        Newton Polynomials II
• y2 = a0 + a1(x2 – x0) + a2(x2 – x0)(x2 – x1)
                                           y1  y0
                                y2  y0           ( x2  x0 )
     y2  a0  a1 ( x2  x0 )              x1  x0
a2                           
      ( x2  x0 )(x2  x1 )         ( x2  x0 )(x2  x1 )

• y2 = a0 + a1(x2 – x0) + a2(x2 – x0)(x2 – x1)
• Data determine coefficients
• Develop scheme known as divided
  difference table to compute ak
                                                                 21
     Divided Difference Table
x0   y0   a0
               y1  y0
          F0             a1
               x1  x0
                               F1  F0   a2
x1   y1                   S0 
                               x2  x0
               y2  y1                        S1  S 0
          F1                            T0 
               x2  x1                        x3  x0
                               F2  F1
x2   y2                   S1                  a3
                               x3  x1
               y3  y 2
          F2 
               x3  x2
x3   y3
                                                         22
    Divided Difference Example
0     0     a0
                10  0
           F0          1 a1
                10  0
                                 3 1
10    10                   S0          .1 a2
                                20  0
               40  10                         .2  .15   1
          F1          3                 T0           
               20  10                          30  0 600
20    40                       63                  a3
                         S1           .15
                              30  10
              100  40
         F2           6
               30  20
30    100
                                                         23
Divided Difference Example II
• Divided difference table gives a0 = 0, a1
  = 1, a2 = .1, and a3 = 1/600
• Polynomial p(x) = a0 + a1(x – x0) + a2(x
  – x0)(x – x1) + a3(x – x0)(x – x1)(x – x2) =
  0 + 1(x – 0) + 0.1(x – 0)(x – 10) +
  (1/600)(x – 0)(x – 10)(x – 20) = x +
  0.1x(x – 10) + (1/600)x(x – 10)(x – 20)
• Check p(30) = 30 + .1(30)(20) + (1/600)
  (30)(20)(10) = 30 + 60 + 10 = 100 (correct)
                                             24
        Constant Step Size
• Divided differences work for equal or
  unequal step size in x
• If Dx = h is a constant we have simpler
  results
  – Fk = Dyk/h = (yk+1 – yk)/h
  – Sk = D2yk/h2 = (yk+2 – 2yk-1 + yk)/h2
  – Tk = D3yk/h3 = (yk+3 – 3yk+2 + 3yk+1 – yk)/h3
  – Dnyk is called the nth forward difference
  – Can also define backwards and central
    differences
                                                    25
   Interpolation Approaches
• When we have N data points how do we
  interpolate among them?
  – Order N-1 polynomial not good choice
  – Use piecewise polynomials of lower order
    (linear or quadratic)
  – Can match first and or higher derivatives
    where piecewise polynomials join
  – Cubic splines are piecewise cubic
    polynomials that match first and second
    derivatives (as well as values)

                                                26
                         Cubic Spline Interpolation
           0.7

           0.6

           0.5
y values




                                                               Known f'
           0.4
                                                               Natural
                                                               No Knot
           0.3                                                 Data


           0.2

           0.1

            0
                 0   1       2      3       4     5   6

                                 x values

                                                          27
                    Newton Interpolating Polynomial
           5


           4


           3                 Polynomial
Y Values




                             Data
           2


           1


           0


           -1
                0      1      2       3      4   5         6

                                  X Values

                                                      28
     Polynomial Applications
• Data interpolation
• Approximation functions in numerical
  quadrature and solution of ODEs
• Basis functions for finite element
  methods
• Can obtain equations for numerical
  differentiation
• Statistical curve fitting (not discussed
  here) usually used in practice
                                             29
              Derivative Expressions
     • Obtain from differentiating interpolation
       polynomials or from Taylor series
     • Series expansion for f(x) about x = a
                  df                  1 d2 f                        1 d3 f
f ( x)  f (a )           ( x  a)                 ( x - a) 2                     ( x - a)3  ....
                  dx x  a            2! dx2   xa
                                                                    3! dx3     xa



    • Note: d0f/dx0 = f                                 
                                                        1 dn f
                                          f ( x)                          ( x - a) n
      and 0! = 1                                   n 0 n ! dx n     x a


    • What is error from truncating series?
                                                                                         30
                    Truncation Error
 • If we truncate series after m terms
                                                  
                m
                    1 dn f                               1 dn f
      f ( x)                      ( x - a) n                          ( x - a) n
               n 0 n ! dx n   x a             n  m 1 n ! dx n   xa

        Terms used      Truncation error, em
 • Can write truncation error as single term
   at unknown location (derivation based
   on the theorem of the mean)
        
             1 dn f                          1 d m1 f
em                         ( x - a) n                                     ( x - a) m1
    n  m 1 n ! dx n    xa
                                          (m  1) ! dx m1            x 
                                                                                       31
                Derivative Expressions
    • Look at finite-difference grid with equal
      spacing: h = Dx so xi = x0 + ih
    • Taylor series about x = xi gives f(xi + kh)
      = f[x0 + (i+k)h] = fi+k in terms of f(xi) = fi
                          df             1 d2 f                         1 d3 f
f ( xi  kh)  f ( xi )            kh                      (kh) 2                      (kh)3  .....
                          dx x  xi      2! dx 2    x  xi
                                                                        3! dx3   x  xi


  • Compact derivative notation
               df                 d2 f                                   dn f
          fi 
            '
                              fi  2
                                ''
                                                                ... f i  n
                                                                        n

               dx   x  xi        dx       x  xi
                                                                         dx           x  xi
                                                                                               32
           Derivative Expressions II
    • Combine all definitions for compact
      series notation
                          df             1 d2 f                            1 d3 f
f ( xi  kh)  f ( xi )            kh                       (kh)   2
                                                                                            (kh)3  .....
                          dx x  xi      2! dx 2     x  xi
                                                                           3! dx3   x  xi


                                      f i '' (kh) 2 f i ''' (kh)3
            fik    f i  f i ' kh                              .....
                                              2!            3!
  • Use this formula to get expansions for
    various grid locations about x = xi and
    use results to get derivative expressions
                                                                                               33
               Derivative Expressions III
     • Apply general                                                     f i '' (kh) 2 f i ''' (kh)3
                                               fik    f i  f i ' kh                              ..
       equation for k                                                            2!            3!
       = 1 and k = –1
                           f i '' h 2 f i ''' h 3
 f i 1  f i  f i ' h                          .....
                              2!         3!                                      f i '' h 2 f i ''' h 3
                                                        f i 1  f i  f i ' h                         .....
                                                                                    2!         3!
         f i 1  f i f i '' h f i ''' h                 f i 1  f i
fi 
   '
                                           .....                   Ah           Forward
               h          2!          3!                       h
     f i  f i 1 f i '' h f i ''' h           f i  f i 1
fi 
  '
                                    .....                Ah
          h        2!       3!                      h                               Backward
                                                                                                   34
              Derivative Expressions IV
   • Subtract fi+1 and fi-1 expressions
                            f i '' h 2 f i ''' h 3
f i 1  f i  f i ' h                              .....
                               2!                 3!                                   f i '' h 2 f i ''' h 3
                                                              f i 1  f i  f i ' h                         .....
                                                                                          2!         3!
                                    2 f i '' ' h 3 2 f i ' '' h 5
 f i 1  f i 1  2 f i ' h                                       .....
                                       3!             5!
            f i 1  f i 1 f i ''' h 2 f i ''''' h 4           f i 1  f i 1
       fi 
         '
                                                     .....                   Ah 2
                   h           3!            5!                       2h
  • Result called central difference expression
                                                                                                         35
         Order of the Error
• Forward and backward derivative have
  error term that is proportional to h
• Central difference error is proportional
  to h2
• Error proportional to hn called nth order
• Reducing step size by a factor of a
  reduces nth order error by an
                                     n
                           h2   
                 e 2  e1 
                          h     
                                 
                           1                36
    Order of the Error Notation
 • Write the error term for nth error term as
   O(hn)
    – Big oh notation, O, denotes order
    – Recognizes that factor multiplying hn may
      change slightly with h
 First order forward            First order backward
        f i 1  f i                      f i  f i 1
   fi 
    '
                      O ( h)        fi 
                                       '
                                                        O ( h)
              h                                h

Second order central                 f i 1  f i 1
                                fi 
                                 '
                                                      O( h 2 )
                                           2h
                                                                  37
                  Higher Order Derivatives
      • Add fi+1 and fi-1 expressions
                              f i '' h 2 f i ''' h 3
   f i 1  f i  f i ' h                            .....
                                 2!                3!                                       f i '' h 2 f i ''' h 3
                                                                 f i 1  f i  f i ' h                           .....
                                                                                               2!         3!
                                        f i ' ' h 2 2 f i '' '' h 4 2 f i '' '''' h 6
      f i 1  f i 1  2 f i  2                                                    .....
                                           2!           4!                   6!
     f i 1  f i 1  2 f i f i ''' h 2 f i ''''' h 4           f i 1  f i 1  2 f i
fi 
  ''

              h2
                            
                                3!
                                        
                                              5!
                                                        ..... 
                                                                          h2
                                                                                          O h2                    
     • f’’ is second-order, central difference
                                                                                                             38
    Higher Order Directional
• We can get higher order derivative
  expressions at the expense of more
  computations
• Get second order forward and backward
  derivative expressions from fi+2 and fi-2,
  respectively
• Combine with previous expressions for
  fi+1 and fi-1 to eliminate first order error
  term
                                            39
           Specific Taylor Series
• General                             f i '' (kh) 2 f i ''' (kh)3
           f i  k  f i  f i ' kh                              .....
  equation                                    2!            3!
                                                         f i '' h 2       f i '' ' h 3
                         fi 2  fi  2 fi 'h  4                    8                   .....
• k=2                                                        2!               3!
                                                          f i '' h 2       f i ''' h 3
                         f i 2  f i  2 f i 'h  4                  8                  .....
• k = -2                                                      2!               3!
                                                    f i '' h 2            f i ' '' h 3
 • k=3               f i 3  f i  3 f i ' h  9                  27                   .....
                                                       2!                    3!
 • k = -3                                              f i '' h 2           f i '' ' h 3
                        f i 3  f i  3 f i ' h  9                  27                  .....
                                                           2!                     3!
                                                                                               40
       Second Order Forward
• Subtract 4fi+1 from fi+2 to eliminate h2
  term                  '' h
                              2
                                 ''' h
                                       3
                                         
   f i  2  4 f i 1   f i  2 f i h  4 f i
                    '
                                                        8 fi    ...
                                                   2         6      
                            '' h
                                  2
                                         ''' h
                                               3
                                                     
    4 f i  f i h  f i
                     '
                                     fi          ..
                               2            6       
                                                    h     3Second
   f i  2  4 f i 1  3 f i  2hf i '  4 f i '''  ... order
                                                    6
                                                                         error
                        f i  2  4 f i 1  3 f i      '' ' h
                                                                    2
                  fi 
                    '
                                                     fi         ...
                                   2h                         3
                                                                                 41
    Second Order Backwards
• Add 4fi-1 to –fi-2 to eliminate h2 term
                                                 '' h
                                                       2
                                                                ''' h
                                                                      3
                                                                             
   f i 2  4 f i 1    f i  2 f i h  4 f i
                                        '
                                                          8 fi          ...
                                                    2              6        
                        '' h
                              2
                                     ''' h
                                           3
                                                 
   4 f i  f i h  f i
                '
                                 fi          ..
                            2            6      
                                               h       Second     3
   f i  2  4 f i 1  3 f i  2hfi  4 f i  '
                                                  ...     ' ''

                                               6       order
                                                                            error
                           f i  2  4 f i 1  3 f i      ''' h
                                                                        2
                      fi 
                        '
                                                       fi        ...
                                      2h                       3
                                                                                    42
Other Derivative Expressions
• Can continue in this fashion
  – Write Taylor series for fi+1, fi-1, fi+2, fi-2, fi+3,
    fi-3, etc.
  – Create linear combinations with factors that
    eliminate desired terms
  – Eliminate fi term to obtain central difference
  – Keep only terms in fk with k  i for forward
    difference expressions
  – Keep only terms in fk with k  i for forward
    difference expressions
  – See results on page 271 of Hoffman
                                                       43
    Order of Error Examples
• Table 2-1 in notes shows error in first
  derivative for ex around x = 1
  – Using first- and second-order forward and
    second-order central differences
  – Step h = 0.4, 0.2, and 0.1
  – Error ratio for doubling step size
     • 4.01 to 4.02 for central differences
     • 2.07 to 2.15 for first-order forward differences
     • 4.32 to 4.69 for second-order forward

               log  e 2 
                    e  log(e )  log(e )
            n          1
                              2           1

               log  2  log( h2 )  log( h1 )
                     h
                    h                                   44
                    1
                   Roundoff Error
    • Possible in derivative expressions from
      subtracting close differences
    • Example f(x) = ex: f’(x)  (ex+h – ex-h)/(2h)
      and error at x = 1 is (e1+h – e1-h)/(2h) - e
          3.004166  2.722815
       E                      2.718282  4.5 x10 3
                 2(0.1)
                         Second order error
    2.7185536702  2.7180100139
E                                 2.7182818284 59  4.5 x10 9
              2(0.0001 )
   2.7182821002 8724  2.7182815566 0388
E                                        2.718281828  5.9 x10 9
                2(0.0000001 )
                                                            45
                        Figure 2-1. Effect of Step Size on Error

        1.E+01
        1.E+00
        1.E-01
        1.E-02
        1.E-03
        1.E-04
Error




        1.E-05
        1.E-06
        1.E-07
        1.E-08
        1.E-09
        1.E-10
        1.E-11
             1.E-17   1.E-15   1.E-13   1.E-11   1.E-09   1.E-07   1.E-05   1.E-03   1.E-01


                                                 Step Size
                                                                                        46
   Numerical ODE Solutions
• The initial value problem
• Euler method as a prototype for the
  general algorithm
• Local and global errors
• More accurate methods
• Step-size control for error control
• Applications to systems of equations
• Higher-order equations
                                         47
   The Initial Value Problem
• dy/dx = f(x,y) (known) with y(x0) = y0
• Basic numerical approach
  – Use a finite difference grid: xi+1 – xi = h
  – Replace derivative by finite-difference
    approximation: dy/dx  (yi+1 – yi) / (xi+1 – xi)
    = (yi+1 – yi) / h
  – Derive a formula to compute favg the
    average value of f(x,y) between xi and xi+1
  – Replace dy/dx = f(x,y) by (yi+1 – yi) / h = favg
  – Repeatedly compute yi+1 = yi + h favg
                                                   48
                 Notation
• xi is the value of the independent
  variable at point i on the grid
  – Determined from the user-selected value of
    step size (or the series of hi values)
  – Can always specify exactly the
    independent variable’s value, xi
• yi is the value of the numerical solution
  at the point where x = xi
• fi is derivative value found from xi and
  the numerical value, yi. I.e., fi = f(xi, yi).
                                                   49
            More Notation
• y(xi) is the exact value of y at x = xi
  – Usually not known but notation is used in
    error analysis of algorithms
• f(xi,y(xi)) is the exact value of the
  derivative at x = xi
• e1 = y(x1) – y1 is the local truncation
  error
  – This is error for one step of algorithm
    starting from known initial condition
                                                50
   Local versus Global Error
• At the initial condition we know the
  solution, y, exactly
• First step introduces some error
• Remaining steps have single step error
  plus previous accumulated error
• Ej = y(xj) – yj is global truncation error
  – Difference between numerical and exact
    solution after several steps
  – This is the error we want to control

                                               51
           Euler’s Method
• Simplest algorithm, example used for
  error analysis, not for practical use
• Define favg = f(xi, yi)
• Euler’s method algorithm is yi+1 = yi + hifi
  = yi + hi f(xi, yi)
• Example dy/dx = x + y, y = 0 at x = 0
• Choose h = 0.1
• We have x0 = 0, y0 = 0, f0 = x0 + y0 = 0
  x1 = x0 + h = 0.1, y1 = y0 + hf0 = 0
                                            52
    Euler Example Continued
• Next step is from x1 = 0.1 to x2 = 0.2
• f1 = x1 + y1 = 0.1 + 0 = 0.1
• y2 = y1 + hf1 = 0 + (.1)(.1) = .01
•  For dy/dx = x + y, we know the exact
  solution is (x0 + y0 + 1)ex-x0 – x – 1
• For x0 = y0 = 0, y = ex – x – 1
• Look at application of Euler algorithm
  for a few steps and compute the error
                                           53
            Euler Example
xi   yi     fi     y(xi)     E(xi)
0    0      0      0         0
.1   0      .1     .0052     .0052
.2   .01    .21    .0214     .0114
.3   .031   .331   .04986    .01886
.4   .0641 .4641   .091825   .027725
                                       54
                    Euler Example Plot
    0.8

    0.7

    0.6

    0.5

    0.4
y




                                      Exact
                                      Numerical
    0.3

    0.2

    0.1

     0
          0   0.1   0.2   0.3   0.4   0.5    0.6   0.7   0.8   0.9        1
                                      x

                                                                     55
         Error Propagation
• Behavior of Euler algorithm is typical of
  all algorithms for numerical solutions
• Error grows at each step
• We usually do not know this global
  error, but we would like to control it
• Look at local error for Euler algorithm
• Then discuss general relationship
  between local and global error

                                              56
           Taylor Series to Get Error
     • Expand y(x) in Taylor series about x = a
                  dy                   1 d2y                       1 d3y
y ( x)  y (a )            ( x  a)                 ( x - a) 2                  ( x - a) 3  ....
                  dx   x a            2! dx 2   x a
                                                                   3! dx 3   xa

   • Look at one step from known initial
     condition, a = x0, to x0 + h so x – a = h
                               dy      1 d2y 2 1 d3y 3
    y ( x 0  h)  y ( x 0 )       h      2
                                               h       3
                                                           h  ....
                               dx 0    2! dx 0     3! dx 0
   • In ODE notation, dy/dx|0 = f(x0, y(x0))
                                                                                         57
                Local Euler Error
• Result of Taylor series on last chart
                                             1 d2y 2 1 d3y 3
y( x0  h)  y( x0 )  hf ( x0 , y ( x0 ))       2
                                                     h       3
                                                                 h  ....
                                             2! dx 0     3! dx 0

 Euler Algorithm          Truncation Error
• This is only the Euler algorithm for the
  first step when we know f(x0,y(x0))
• This gives the local truncation error
• Local truncation error for Euler algorithm
  is second order
                                                                    58
              Global Error
• We will show that a local error of order
  n, has a global error of order n-1
• To show this consider the global error at
  x = x0 + kh after k algorithm steps
  – Is approximately k times the local error
  – If local error is O(hn), approximate global
    error after k steps is k O(hn)  kAhn
  – A new step size, h/r, takes kr steps to get
    to the same x value
                                                  59
     Global Error Concluded
• Compare error for same x = kh with
  step sizes h and h/r
• Ex=kh(h)  kAhn                   
                          Ex kh h
                                   r 
                                       kr h
                                               n

                                            r  1
• Ex=kh(h/r) = krA(h/r) n
                                         k h 
                                             n   n 1
                            Ex kh (h)                 r

• When we reduce the step size by a
  factor of 1/r we reduce the error by a
  factor of 1/rn-1; this is the behavior of
  an algorithm whose error is order n-1
                                                  60
 Euler Local and Global Error
• Previously showed Euler algorithm to
  have second order local error
• Should have first order global error
• Results for previous Euler example at x
  = 1 with different step sizes
    Step size      First step   Final error
      h = 0.1     5.17x10-3     1.25 x10-1
     h = 0.01     5.02 x10-5    1.35 x10-2
    h = 0.001     5.00 x10-7    1.36 x10-3
                                              61
          Better Algorithms
• Seek high accuracy with low
  computational work
• Could improve Euler accuracy by
  cutting step size, but this is not efficient
• Use other algorithms that have higher
  order errors
• Runge-Kutta methods typically used
  – This is a class of methods that use several
    function evaluation methods per step

                                                 62
            Second-order Runge Kutta
        • Huen’s method
yi01  yi  hi 1 f ( xi , yi )                      xi 1  xi  hi 1

yi 1    yi 
               hi 1
                2
                       
                     f ( xi , yi )  f ( xi 1 , yi01 )   
                                                           yi  yi01  hi 1 f ( xi 1 , yi01 )
                                                                          2

        • Modified Euler method
                         hi 1                                                      hi 1
        yi 1  yi              f ( xi , y i )                 xi  1        xi 
              2          2                                               2           2
        y i 1  y i  hi 1 f ( xi  1 , y i  1 )
                                       2        2

                                                                                              63
      Fourth-order Runge Kutta
• Uses four derivative evaluations per step
               k1  2 k 2  2 k 3  k 4
y i 1  y i                                  xi 1  xi  hi 1
                          6
k1  hi 1 f ( xi , y i )

k 2  hi 1 f  xi  i 1 , y i  1 
                     h                  k
                                         2
                           2              
               x  hi 1 , y  k 2 
k 3  hi 1 f  i
                           2 i           2
                                           
k 4  hi 1 f ( xi  hi 1 , y i  k 3 )
                                                            64
    Comparison of Methods
• Look at Euler, Heun, Modified Euler and
  fourth-order Runge-Kutta
• Solve dy/dx = e-y-x with y(0) = 1
• Compare numerical values to exact
  solution y = ln( ey0 + e-x0 – e-x)
• Look at errors in the methods at x = 1
  as a function of step size
• Compare error propagation (increase in
  error as x increases)
                                        65
                      Error versus Step Size for Simple ODE Solvers
        1.E+00
                  dy
                      e x y        with y  y0  1 at x  x0  0
        1.E-01
        1.E-02
        1.E-03
                  dx
        1.E-04
        1.E-05
        1.E-06
        1.E-07                                                                       Euler
Error




        1.E-08                                                                       Huen
        1.E-09                                                                       Modified Euler
                                                                                     Runge-Kutta
        1.E-10
        1.E-11
        1.E-12
        1.E-13
        1.E-14
        1.E-15
        1.E-16
             1.E-07      1.E-06   1.E-05      1.E-04      1.E-03   1.E-02   1.E-01
                                           Step size, h
                                                                                          66
                     Error Propagation in Solutions of ODEs

        1.E+00
        1.E-01

        1.E-02
        1.E-03
        1.E-04
                                                              Euler (N = 10)
        1.E-05                                                Euler (N = 100)
        1.E-06                                                Euler N = 1,000)
Error




                                                              Huen (N =10)
        1.E-07
                                                              Huen (N = 100)
        1.E-08                                                Huen (N = 1,000)
        1.E-09                                                RK4 (H = 10)
                                                              RK4 (N = 100)
        1.E-10
        1.E-11
        1.E-12
        1.E-13
        1.E-14
                 0    0.2      0.4       0.6    0.8     1           67
                                     x
     Some Basic Concepts
• A numerical method is convergent with
  the solution of the ODE if the numerical
  solution approaches the actual solution
  as h  0 (with increase in numerical
  precision at smaller h)
  – Mainly a theoretical concept; algorithms in
    use are convergent
• Stability refers to the ability of a
  numerical algorithm to damp any errors
  introduced during the solution
                                                  68
         More on Stability
• Finite difference equations in numerical
  algorithms, when iterated, may be
  numerically increase without bound
• Stability usually is obtained by keeping
  step size h small, sometimes smaller
  than the h required for accuracy
• For most ODEs stability is not a
  problem, but it is for stiff systems of
  ODEs and for partial differential
  equations
                                         69
              Euler Solution of Decaying Exponential
     5

    4.5

     4                                                         y(0) =   1
                                                               y(0) =   2
    3.5                Euler algorithm errors move
                                                               y(0) =   3
                      numerical solution of dy/dx = -y
     3                                                         y(0) =   4
                       to solution for anohter initial
                       condition. Error bounded by             y(0) =   5
y




    2.5
                          decreasing exponential               Euler
     2

    1.5

     1

    0.5

     0
          0   0.5      1           1.5          2        2.5                3
                                   x


                                                                            70
                    Euler Solution for Increasing Exponential
    100
                                                       y(0) =   1
     90
                                                       y(0) =   2
     80        Euler algorithm errors move             y(0) =   3
              numerical solution of dy/dx = -y
     70                                                y(0) =   4
               to solution for anohter initial
     60                                                y(0) =   5
              condition. Error becomes large
                for increasing exponential             Euler
y




     50

     40

     30

     20

     10

      0
          0           0.5          1             1.5   2            2.5   3

                                                 x


                                                                              71
            Error Control
• How do we choose h?
• Want to obtain a result with some
  desired small global error
• Can just repeat calculations with smaller
  h until two results are sufficiently close
• Can use algorithms that estimate error
  and adjust step size during the
  calculation based on the error

                                          72
     Runge-Kutta-Fehlberg
• Uses two equations to compute yn+1,
  one has O(h5), the other O(h6) error
• Requires six derivative evaluations per
  step (same evaluations used for both
  equations)
• The error estimate can be used for step
  size control based on an overall 5th
  order error
• Cask-Karp version different coefficients

                                         73
      Runge-Kutta-Fehlberg II
•   See page 377 in Hoffman for algorithm
•   Typical formula components below
•   yn+1= yn + (16k1/135 + 6656k2/12825 …
•   k3 = hf(xn + 3h/8, yn + 3k1/32 + 9k2/32)
•   Error = k1/360 - 128k3/4275 …
•   hnew = hold|EDesired/Error|1/5
•   EDesired is set by user
•   RKF45 code by Watts and Shampine
                                               74
       Other Error Controls
• Do fourth-order Runge-Kutta two times,
  with step sizes different by a factor of
  two and compare results
• Runge-Kutta-Verner is similar to Runge-
  Kutta-Fehlberg, but uses higher order
  expressions




                                         75
      Systems of Equations
• Any problem with one or more ODEs of
  any order can be reduced to a system
  for first order ODEs
• Last week we reduced a system of two
  second order equations to a system of
  four first order equations
• In this process the sum of the orders is
  constant
• Example from last week

                                         76
    Reduction of Order Example
  d 2 y1 k1  k2      k2              d 2 y2 k2       k3  k 2
      2
                y1     y2  0           2
                                                y1           y2  0
   dt      m1         m1               dt     m2        m2
• Define variables y3 = dy1/dt and y4 = dy2/dt
• Then dy3/dt = d2y1/dt2 and dy4/dt = d2y2/dt2
• Have four simultaneous first-order ODEs
                dy1                       dy2
                     y3  f1 ( x, y )          y 4  f 2 ( x, y )
                dt                         dt
                dy3     k1  k 2       k2
                               y1      y 2  f 3 ( x, y )
                dt        m1           m1
                dy4 k 2      k3  k 2
                       y1           y 2  f 4 ( x, y )
                 dt m2         m2                                     77
 Solving Simultaneous ODEs
• Apply same algorithms used for single
  ODEs
  – Must apply each step and substep to all
    equations in system
  – Key is consistent determination of fi(x,y)
  – All yi values in y must be available at the
    same x point when computing the fi
  – E.g., in Runge-Kutta we must evaluate k1
    for all equations before finding k2

                                                  78
Runge-Kutta for ODE System
 – y(i) is vector of dependent variables at x = xi
 – k(1), k(2), k(3), and k(4), are vectors containing
   intermediate Runge-Kutta results
 – f is a vector containing the derivatives
 – k(1) = hf(xi, y(i))
 – k(2) = hf(xi + h/2, y(i) + k(1)/2)
 – k(3) = hf(xi + h/2, y(i) + k(2)/2)
 – k(4) = hf(xi + h, y(i) + k(3))
 – y(i+1) = (k(1) + 2k(2) + 2k(3) + k(4))/6
                                                   79
       ODE System by RK4
• dy/dx = -y + z and dz/dx = y – z with
  y(0) = 1 and z(0) = -1 with h = .1
• k(1)y = h[-y + z] = 0.1[-1 + (-1)] = -.2
• k(1)z = h[y - z] = 0.1[1 - (-1)] = .2
• k(2)y = h[-(y+ k(1)y/2) + z + k(1)z/2] = 0.1[
  -(1 + -0.2/2) + (-1 + .2/2)] = -.18
• k(2)z = h[(y+ k(1)y/2) –(z + k(1)z/2)] = 0.1[(1
  + -0.2)/2 - (-1 + .2/2)] = .18

                                                80
      ODE System by RK4 II
• k(3)y = h[-(y+ k(2)y/2) + z + k(2)z/2] = 0.1[
  -(1 + -0.18/2) + (-1 + .18/2)] = -.182
• k(3)z = h[(y+ k(2)y/2) –(z + k(2)z/2)] = 0.1[(1
  + -0.18)/2 - (-1 + .18/2)] = .182
• k(4)y = h[-(y+ k(3)y) + z + k(3)z] = 0.1[ -(1 +
  -0.182) + (-1 + .182)] = -.1636
• k(4)z = h[(y+ k(3)y) –(z + k(3)z)] = 0.1[ (1 +
  -0.182) - (-1 + .182)] = .1636

                                               81
     ODE System by RK4 III
• yi+1 = yi + (k(1)y + 2k(2)y + 2k(3)y + k(4)y )/6
  = 1 + [ (–.2) + 2(–.18) + 2(–.182) +
  (–.1636)]/6 = .8187
• zi+1 = zi + (k(1)z + 2k(2)z + 2k(3)z + k(4)z )/6
  = –1 + [(.2) + 2(.18) + 2(.182) +
  (.1636)]/6 = –.8187
• Continue in this fashion until desired
  final x value is reached
• No x dependence for f in this example
                                                 82
Numerical Software for ODEs
• Usually written to solve a system of N
  equations, but will work for N = 1
• User has to code a subroutine or
  function to compute the f array
  – Input variables are x and y; f is output
  – Some codes have one dimensional
    parameter array to pass additional
    information from main program into the
    function that computes derivatives

                                               83
  Derivative Subroutine Example
  • Visual Basic code   dy1
                              y1  y2  y3e 2 x
    for system of       dx
    ODEs at right is    dy2           dy3
                             2 y12
                                           3 y1 y2
    shown below         dx            dx
Sub fsub( x as Double, y() as Double, _
  f() as Double )
  f(1) = -y(1) + Sqr(y(2)) + y(3)*Exp(2*x)
 f(2) = -2 * y(1)^2
  f(3) = -3 * y(1) * y(2)
End Sub                               84
  Derivative Subroutine Example
  • Fortran code for    dy1
                              y1  y2  y3e 2 x
    system of ODEs      dx
    at right is shown   dy2           dy3
                             2 y12
                                           3 y1 y2
    below               dx            dx
subroutine fsub( x, y, f )
  real(KIND=8) x, y(*), f(*)
  f(1) = -y(1) + sqrt(y(2)) +y(3)*exp(2*x)
 f(2) = -2 * y(1)**2
  f(3) = -3 * y(1) * y(2)
end subroutine fsub                   85
   Derivative Subroutine Example
   • C++ code for        dy1
                               y1  y2  y3e 2 x
     system of ODEs      dx
     at right is shown   dy2           dy3
                              2 y12
                                            3 y1 y2
     below               dx            dx
void fsub(double x, double y[], double f[])
{
 f[1] = -y[1] + sqrt(y[2]) + y[3]*exp(2*x);
  f[2] = -2 * y[1] * y[1];
  f[3] = -3 * y[1] * y[2];
}                                       86

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:6
posted:3/13/2012
language:English
pages:86