# Basics of Ordinary Differential Equations

Document Sample

```					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 z2     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 differentiation
• Numerical differential equations
18
Interpolation
• 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
y 2  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
– 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
• 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! dx 2   xa
3! dx 3    xa

• Note: d0f/dx0 = f                                  
1 dn f
f ( x)              n
( x - a) n
and 0! = 1                                    n  0 n ! dx       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! dx 3   x  xi

• Compact derivative notation
df                   d2 f                                    dn f
fi 
'
fi  2
''
...   fi  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! dx 3   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 ''' h3
f i 1  f i  f i ' h                         .....
2!         3!                                    f i '' h 2 f i ''' h3
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 ''' h3
f i 1  f i  f i ' h                            .....
2!               3!                                  f i '' h 2 f i ''' h3
f i 1  f i  f i ' h                        .....
2!         3!
2 f i ''' h3 2 f i ''' h5
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 ''' h3
f i 1  f i  f i ' h                           .....
2!               3!                                    f i '' h 2 f i ''' h3
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 ''' h3
• k=2                    fi2  fi  2 fi 'h  4                      8                .....
2!                3!
f i '' h 2      f i ''' h3
• k = -2                 f i 2  f i  2 f i ' h  4                 8                .....
2!               3!
f i '' h 2            f i ''' h3
• k=3               f i 3  f i  3 f i ' h  9                   27                .....
2!                    3!
• k = -3                                             f i '' h 2            f i ''' h3
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      
                    h2         h3 
 4  f i  f i ' h  f i ''  f i '''  ..
                    2          6      
h 3    Second
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        
                        h2         h3 
 4  f i  f i ' h  f i ''  f i '''  ..
                         2         6      
h       Second    3
 f i 2  4 f i 1  3 f i  2hf i  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
e
log  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.718281828459  4.5 x109
2(0.0001)
2.71828210028724  2.71828155660388
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 E (h)
k h 
n     n 1
x  kh                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
y i01  y i  hi 1 f ( xi , y i )                    xi 1  xi  hi 1

y i 1    yi 
hi 1
2

f ( xi , y i )  f ( xi 1 , y i01 ) y i  y i01  hi 1 f ( xi 1 , y i0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  2k 2  2k 3  k 4
y i 1  y i                                 xi 1  xi  hi 1
6
k1  hi 1 f ( xi , y i )
 x  hi 1 , y  k1 
k 2  hi 1 f  i
             2 i          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

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 k 2      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
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: 5 posted: 11/27/2011 language: English pages: 86
How are you planning on using Docstoc?