# Numerical Solution of Ordinary Differential Equations

Document Sample

```					                                                             1

Numerical Solution of Ordinary Differential Equations

Consider the 1st order ordinary differential equation (ODE)

dy
 f ( x, y) .
dx

The initial condition can be taken as

y ( x0 )  y 0 .

Then we could use a Taylor series about x  x0 and obtain the complete solution, or

1                          1                                 1 (N)
y( x)  y( x0 )  y' ( x0 )(x  x0 )          y' ' ( x0 )(x  x0 ) 2  y' ' ' ( x0 )(x  x0 ) 3  ...     y ( x0 )(x  x0 ) N  ...
2!                         3!                                N!

Since

y ( x0 )  y 0
and

y ' ( x0 )  f ( x0 , y 0 )  y 0 '

then we can find the first two terms. For the second derivative

d2y                  f f dy 
y ''( x0 )  2                                    f, x ( x0 , y0 )  f, y ( x0 , y0 ) y0 ' .
dx        x  x0     x y dx  x x 0


Similarly the third derivative is

d3y
y '''( x0 )                    f, xx ( x0 , y0 )  2 f, xy ( x0 , y0 ) y0 ' f, yy ( x0 , y0 ) y0 '2 .
dx3   x  x0

If we truncate at the third derivative

1                          1
y( x)  y( x0 )  y' ( x0 )(x  x0 )       y' ' ( x0 )(x  x0 ) 2  y' ' ' ( x0 )(x  x0 ) 3  Error ,
2!                         3!

and

1 iv
Error            y ( )(x  x0 ) 4 , where x0    x .
4!
2

Euler’s Method

Take the Taylor series to 1st order, and let the interval h  x1  x0 , then

y ' ' ( ) 2
y1  y ( x0 )  f x0 , y0 h               h .
2
 
The error for a time step (the local error) is O h 2 . The global error, after many steps, is
Oh  . Then

y1  y 0  f ( x0 , y 0 )h where y0  y ( x0 ) ,
y 2  y1  f ( x1 , y1 )h where x1  x0  h ,
…

y N 1  y N  f ( x N , y N )h where x N 1  x N  h .

Example:
dy
 x  y, y (0)  1
dx

The exact solution can be found from
dy
 y  x.
dx

dy
Let y  yc  y p where c  y c  0 , or yc  Ce rx . Then rCe rx  Ce rx  0 for all x ,
dx

or r  1 , and yc  Ce x . Since the right hand side is linear in x try y p  Ax  B . Then

dy p             dy p
 A and           y p  x becomes A  Ax  B  x which must hold for all x. Hence
dx              dx

A  1, and B=-1, making y p  x  1 , and since y  yc  y p then

y  Ce x  x  1 .

But @ x  0 , y  1 or C  1  1 , and C  2 . Making the complete solution

y  2e x  x  1 .
3

Using Euler’s method and taking h  0.02

y 0  1, x0  0  x1  0.02 , y1  1  1  0.02  1.02 , since y 0 '  1 . In general,

y n 1  y n  y n ' h; xn 1  xn  h

n             xn                 yn                    yexact
0            0.00              1.0000                 1.0000
1            0.02              1.0200                 1.0204
2            0.04              1.0408                 1.0416
3            0.06              1.0624                 1.0637
4            0.08              1.0848                 1.0866
5            0.10              1.1081                 1.1103

For the error, y Euler  y5  y 0  0.1081 , y Exact  y5  y 0  0.1103 , can be defined
as

y Euler  y Exact          0.0022
Relative Error                                           2%
1
y Euler  y Exact        0.1092
2

The results plot as

It would be better to use the slope at the beginning and end of the increment (e.g., the
average at each end), and although we don’t know the slope at the end we can
approximate it.
4

Modified Euler Method

Let y n '  f ( xn , y n ) . Then an approximation for y at the end of the increment is
~
y n 1  y n  y n ' h
and an estimate for the slope at the end of the increment is
~'
y         f (x            ,~y     ).
n 1        n 1     n 1

We can now set
y n 1  y n 
1
2
            
y n ' ~n 1 h .
y'
The error can be found from
'
y n 1  y n  y n h 
1 '' 2
2
yn h  O h3 
and since
1  y n 1  y n        
 
'        '
'
y n 1  y n  y n h                     Oh  h 2  O h 3
2       h             

or
 y '  yn 
 
'
y n 1  y n   n 1     h  O h 3 .
     2    
          
Hence the local error is Oh 3  and the global error is Oh 2 . Another way to write our
results is

k1  hf ( xn , y n )
k 2  hf ( x n  h, y n  k1 )

y n 1  y n 
1
k1  k 2 
2

The previous example now can include modified Euler

n           xn                yeuler                    ymodified             yexact
0          0.00              1.0000                     1.0000               1.0000
1          0.02              1.0200                     1.0204               1.0204
2          0.04              1.0408                     1.0416               1.0416
3          0.06              1.0624                     1.0637               1.0637
4          0.08              1.0848                     1.0866               1.0866
5          0.10              1.1081                     1.1104               1.1103

which is much better.
5

Runge-Kutta Methods

The modified Euler method is actually a two step (second order) Runge-Kutta algorithm.
These methods can be readily extended to eight and even tenth order. The derivations
follow the same procedure. Assume for the second order method

y n 1  y n  ak1  bk 2
where
k1  hf ( xn , y n ) , and

k 2  hf ( xn  h, y n  k1 ) .

The parameters a, b,  and  are found by comparing to a Taylor series expansion.
Recall
1
2

y n 1  y n  hf n  h 2 f n, x  f n, y f n .
But
f ( xn  h, yn  k1 )  f n  hf n, x  k1 f n, y .
Since k1  hf n ,
f ( xn  h, yn  k1 )  f n  hf n, x  hf n f n, y
Using y n 1  y n  ak1  bk 2 ,


yn 1  yn  ahf n  bh f n  hf n, x  hf n f n, y ,  
or
y n 1  y n  a  bhf n  bh 2 f n, x  bh 2 f n f n, y

Comparing to the Taylor series

1            1
y n 1  y n  hf n  h 2 f n, x  h 2 f n f n, y .
2            2

1        1
Then a  b  1, b      , b  which is three equations in four unknowns. Hence we
2        2
can pick for the fourth equation any equation that is convenient.

For example we can take
a  2 / 3, b  1 / 3, and     3 / 2
or
a  0, b  1, and     1 / 2 .
Modified Euler is
a  b  1 / 2, and     1 .
6

Fourth Order Runge-Kutta

For fourth order Runge-Kutta the estimates for the changes are

k1  hf ( xn , y n )
1           1
k 2  hf ( x n     h, y n  k1 )
2           2
1           1
k3    hf ( xn  h, y n  k 2 )
2           2
k4    hf ( x n  h, y n  k 3 )

and the updated value for y is found from

y n 1  y n 
1
k1  2k 2  2k3  k 4 .
6

Note that all of the algorithms presented are for first order equations with only one
dependent variable. These can be readily extended to systems of higher order differential
equations. For those cases please see page 7.
7

Higher Order Differential Equations

Consider
dny
n
dx
                           
 y (n) ( x)  f x, y, y' , y ' ' ,..., y n 1  f x, y0 , y1 , y 2 ,..., y n 1 

and let the vector Y be defined as
 y0   y
'
          y0 
 '   1                               y 
 y1   y 2                            1 
      
  y'   y                              y 
Y '  2    3                         Y  2 
 ...              ...  where     ... 
 yn  2 
'              y n 1           yn  2 
 '              f               y       
 y n 1 
                                 n 1 

which is now a vector first order equation, and the first order rules can be applied to a
system of first order equations.

Systems of First Order Equations

 
      
Let Y   F x, Y
            
Euler: YN 1  YN  hFN

            
          1    
Modified Euler:                      YN 1  YN  K1  K 2
2
where
            
             
K1  hF x N , YN
                          
                    
K 2  hF x N  h, YN  K1

Runge-Kutta (4th order):             YN 1  YN 
1
6

K1  2K 2  2K3  K 4               
where
8

       
            
K1  hF x N , YN
               1      1  
K 2  hF  x N  h, YN  K1 
        2       2 
               1      1  
K 3  hF  x N  h, YN  K 2 
       2      2 
                

K 4  hF x N  h, YN  K 3
9

Problem

The equation for a pendulum with a mass, m, at the end of a rod of negligible mass with
length, L, is


mL2  mgL sin   0 .

For the initial conditions take

   0 @ t  0 , and   0 @ t  0 ,                            (1)

where
d( )              d 2( )
 is the angle from the vertical, ()            , and ()         .
dt                dt 2
g
Let       t . The equation becomes
L
d 2
 sin   0 .                             (2)
d 2
d
Now let p       then since
d
d 2         dp dp d    dp
           p    ,
d   2       d d d    d

The equation can be written
dp
p        sin   0 ,
d
and integrating
1 2
p  cos  E  constant .
2
d
This is the conservation of energy. The initial condition is p               0 @   0 ,
d
hence E   cos 0 or
2
1  d 
   cos   cos 0                                      (3)
2  d 

is also the governing differential equation.

Integrate the equation of motion (2) subject to the initial conditions (1),

   / 2 @   0 and   0 @   0 using Euler, modified Euler, and Runge-Kutta. Use
the expression for the energy to check the accuracy of your integration. Integrate for one-
quarter of a cycle and then determine the period for a complete cycle.

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 4 posted: 12/2/2011 language: English pages: 9