Transcript Slide 1

Solutions to ODEs
• A general form for a first order ODE is
dy 

F  x, y ,   0
dx 

• Or alternatively
dy
 f  x, y 
dx
...(1)
• A general form for a first order ODE is
dy 

F  x, y ,   0
dx 

• Or alternatively
dy
 f  x, y 
dx
...(1)
• We desire a solution y(x) which satisfies
(1) and one specified boundary condition.
• To do this we divide the interval in the
independent variable x for the interval [a,b]
into subintervals or steps.
• To do this we divide the interval in the
independent variable x for the interval [a,b]
into subintervals or steps.
• Then the value of the true solution is
approximated at n+1 evenly spaced values
of x.
• Such that, xi  x0  ih i  0,1,...,n
ba
• where,h 
n
• We denote the approximation at the base
pts by yi so that yi  yxi 
• We denote the approximation at the base
pts by so that yi  yxi 
• The true derivation dy/dx at the base
points is approximated by f xi , yi  or just f i
• where fi  f xi , yi   f xi , yxi 
• We denote the approximation at the base
pts by so that yi  yxi 
• The true derivation dy/dx at the base
points is approximated by f xi , yi  or just f i
• where fi  f xi , yi   f xi , yxi 
• Assuming no roundoff error the difference
in the calculated and true value is the
truncation error,  i  yi  yxi 
• The only errors which will be examined are
those inherent to the numerical methods.
•
Numerical algorithms for solving 1st order
odes for an initial condition are based on
one of two approaches:
1. Direct or indirect use of the Taylor series
expansion for the solution function yx 
2. Use of open or closed integration
formulas.
•
The various procedures are classified
into:
1. One-step: calculation of yi 1 given the
differential equation and yi 1, xi
2. Multi-step: in addition the previous
information they require values of x and y
outside of the interval under
consideration
• The first method considered is the Taylor
series method.
• It forms the basis for some of the other
methods.
Taylor Expansion
• For this method we express the solution yx 
about some starting point x0 using a Taylor
expansion.
y  x0  h   y  x0   hf  x0 , y  x0 
h2

f  x0 , y  x0 
2
h3

f  x0 , y  x0   ...
6
• where f x, yx denotesd dx f x, yx
• Where f x, yx denotesd dx f x, yx
2
2






f
x
,
y
x
denotes
d
dx
f x, yx 
• And
• etc.


dy
2



f
x
,
y

x
where
dx
• Consider the example,
• And initial condition yx0   y0
dy
2



f
x
,
y

x
where
dx
• Consider the example,
• And initial condition yx0   y0
• Differentiating and then applying initial
conditions,
f x, y   2 x
dy
2



f
x
,
y

x
where
dx
• Consider the example,
• And initial condition yx0   y0
• Differentiating and then applying initial
conditions,
f x, y   2 x
f x0 , y0   2 x0
dy
2



f
x
,
y

x
where
dx
• Consider the example,
• And initial condition yx0   y0
• Differentiating and then applying initial
conditions,
f x, y   2 x
f x, y   2
f x0 , y0   2 x0
dy
2



f
x
,
y

x
where
dx
• Consider the example,
• And initial condition yx0   y0
• Differentiating and then applying initial
conditions,
f x, y   2 x
f x0 , y0   2 x0
f x, y   2
f x0 , y0   2
• Continuing,
f x, y   0
• Continuing,
f x, y   0
f x0 , y0   0
• Continuing,
f x, y   0
.
.
f n  x, y   0
f x0 , y0   0
.
.

n
x0 , y0   0 n  3
f
Here the third derivative and higher are
zero.
• Substituting these values into Taylor
expansion gives,
3
h
y x0  h   y x0   hx02  h 2 x0 
3
• Substituting these values into Taylor
expansion gives,
3
h
y x0  h   y x0   hx02  h 2 x0 
3
• To determine the error lets consider
compare this to the analytic solution.
• Substituting these values into Taylor
expansion gives,
3
h
y x0  h   y x0   hx02  h 2 x0 
3
• To determine the error lets consider
compare this to the analytic solution.
• To do this we will integrate the differential
equation.
dy
• Integrating,  x 2
dx
dy
• Integrating,  x 2
dx
y  x0  h 
 
dy 
y x0
x0  h

x0
x 2 dx
dy
• Integrating,  x 2
dx
y  x0  h 
 
dy 
y x0
x0  h

x0
3
x
x dx 
3
2
x0  h


 x0
dy
• Integrating,  x 2
dx
y  x0  h 
 
dy 
y x0
•
x0  h

x0
3
x
x dx 
3
2
x0  h


 x0
3
h
2
2




y
x

h

y
x

hx

h
x0 
Simplifying, 0
0
0
3
• In this case the two expressions are
identical, therefore there is no truncation
error.
Homework (Due Thursday)
dy
Consider the function  f x, y   2 y
dx
Initial condition y x0   y0
•
•
1. Write the Taylor expansion and show that
it can be written in the form,
2
3
4



2h
2h 2h
yx0  h  yx0 1  2h 



2!
3!
4! 

n 1
f
2. Show if terms up to n 1
are retained,

2h 
the error is,  y  
,  in x0 , x0  h 
n  1!
• Stepping from x0 to x0  h follows from the
Taylor expansion of y x  about x0.
• Stepping from x0 to x0  h follows from the
Taylor expansion of y x  about x0.
• Equivalently stepping from x i to xi  h can
be accomplished from the Taylor
expansion of y x  about x i .
• Stepping from x0 to x0  h follows from the
Taylor expansion of y x  about x0.
• Equivalently stepping from x i to xi  h can
be accomplished from the Taylor
expansion of y x  about x i .
y xi 1   y x  h   y xi   hf xi , y xi 
h 2 f xi , y xi  h 3 f xi , y xi 

2!

3!
 ...
h n f n 1 xi , y xi  h n 1 f n   , y  


n  1!
n!
 in xi , xi 1 
• Algorithms for which the last term in
expansion is dropped are of order hn.
• The error is or order hn+1.
• The local truncation
error
is
bounded
as
h n 1
M
follows et 
n  1!
• where
M  f n   , y  
max
 inxi , xi 1 
• However differentiation of f x, y  can be
complicated.
• However differentiation of f x, y  can be
complicated.
• Direct Taylor expansion is not used other
than the simplest case,
 
yxi 1   yxi   hf xi , yxi   O h 2
• Big Oh is the order of the algorithm.
• Usually only the value of y x0  is the only
value of yxi  that is known, therefore yxi 
must be replaced by y i .
• Thus,
y1  yx0   hf x0 , yx0 
• In general, yi 1  yi  hf xi , yxi   yi  hfi i  1
which is known as Euler’s method.
Euler’s Method
• The problem with the simple Euler method
is the inherent inaccuracy in the formula.
• The geometric interpretation is shown the
diagram.
y
y1
y(x)
y(x0)
y(x1
h
x0
x1
x
• The solution across the interval[x0,x1] is
assumed to follow the line tangent to y(x)
at x0.
y
y1
y(x)
y(x0)
y(x1
h
x0
x1
x
• When the method is applied repeatedly
across several intervals in sequence, the
numerical solution traces out a polygon
segment with sides of slope fi,
i=0,1,2,…,(n-1).
• The simple Euler method is a linear
approximation, and only works if the
function is linear (or at least linear in the
interval).
• The simple Euler method is a linear
approximation, and only works if the
function is linear (or at least linear in the
interval).
• This is inherently inaccurate.
• Because of this inaccuracy small step
sizes are required when using the
algorithm.
y
y1
i
y(x0)
y(x1
• From the graph, as
h -> 0, the
approximation
y(x) becomes better
because the curve in
the region becomes
approx. linear.
h
x0
x1
x
• The modified Euler method improves on
the simple Euler method.
Modified Euler
• One modification is to determine the
average slope in the region.
• The average slope may be approximated
by the mean of the slopes at the beginning
and the end of the interval.
• One modification is to determine the
average slope in the region.
• The average slope may be approximated
by the mean of the slopes at the beginning
and the end of the interval.
• This modified Euler may be written as the
arithmetic average, yi 1  yi  h  yi  yi1
or yi 1
f i  f i 1
 yi  h 
2
2
dy
where
 f  x, y 
dx
• However it is not possible to implement
this directly.
• Instead we first “predict” a value for yi 1
using the simple Euler method.
• Then use this value to determine the
gradient of the slope at yi 1. This gives a
“corrected” value for yi 1 .
• However since the predicted value is not
usually accurate, yi1 is also inaccurate.
Corrector curve using
values of predictor
Predictor curve using
the Simple Euler
y1
y0
True value
x0
x1
Runge-Kutta Methods
• The solution of a differential equation
using higher order derivatives of the Taylor
expansion is not practical.
• The solution of a differential equation
using higher order derivatives of the Taylor
expansion is not practical.
• Since for only the simplest functions, these
higher orders are complicated. Also there
is no simple algorithm which can be
developed.
• This is because each series expansion is
unique.
• However we have methods which use only
1st order derivates while simulating higher
order(producing equivalent results).
• These one step methods are called
Runge-Kutta methods.
• However we have methods which use only
1st order derivates while simulating higher
order(producing equivalent results).
• These one step methods are called
Runge-Kutta methods.
• Approximation of the second, third and fourth
order (retaining h2, h3, h4 respectively in the
Taylor expansion) require estimation at 2, 3 ,
4 pts respectively in the interval (xi,xi+1).
• Method of order m (where m >4) require
derivate evaluation at more than m nearby
points.
• The Runge-Kutta methods have
algorithms of the form, yi 1  yi  h xi , yi , h
• where is the increment function.
• The increment function is a suitably
chosen approximation to f x, y  on the
interval xi  x  xi 1 .
• The Runge-Kutta methods have
algorithms of the form, yi 1  yi  h xi , yi , h
• where is the increment function.
• The increment function is a suitably
chosen approximation to f x, y  on the
interval xi  x  xi 1 .
• Because of the amount of algebra
involved, on the simplest case(2nd order)
will be derived in detail.
Derivation of the 2nd order
Runge-Kutta
• Let  be the weighted average of two
derivative evaluations k1and k 2 on the
interval xi  x  xi 1,
• Let  be the weighted average of two
derivative evaluations k1and k 2 on the
interval xi  x  xi 1,   ak1  bk2
• Let  be the weighted average of two
derivative evaluations k1and k 2 on the
interval xi  x  xi 1,   ak1  bk2
• Then, yi 1  yi  hak1  bk2 
• Let  be the weighted average of two
derivative evaluations k1and k 2 on the
interval xi  x  xi 1,   ak1  bk2
• Then, yi 1  yi  hak1  bk2 
• Let, k1  f xi , yi 
k 2  f xi  ph, yi  qhf xi , yi 
 f xi  ph, yi  qhk1 
• The constants p and q will be established.
• Expanding k 2in a Taylor series for a
function of two variables and dropping
terms higher than h,
k2  f xi  ph, yi  qhf xi , yi 
 f xi , yi   phfx xi , yi 
 
 qhf xi , yi  f y xi , yi   O h 2
• NB: 1st few terms in two-variable Taylor series,
f x  r , y  s   f x, y   rf x x, y   sf y x, y   r 2 f xx x, y  / 2

 rsf xy x, y   s 2 f yy x, y  / 2  O  r  s 
3

 
2








k

f
x
,
y

phf
x
,
y

qhf
x
,
y
f
x
,
y

O
h
2
i
i
x
i
i
i
i
y
i
i
Recall:
•
• Substituting for k 2,
yi 1  yi  haf xi , yi   bf xi , yi 

  
 h 2 bpfx xi , yi   bqf xi , yi  f y xi , yi   O h3
 
2








k

f
x
,
y

phf
x
,
y

qhf
x
,
y
f
x
,
y

O
h
2
i
i
x
i
i
i
i
y
i
i
Recall:
•
• Substituting for k 2,
yi 1  yi  haf xi , yi   bf xi , yi 
 h 2 bpfx xi , yi   bqf xi , yi  f y xi , yi   O h3 ...1

  
• Expanding the function yx  about x i
y xi  h   y xi 1   y xi   hf xi , y xi 
h2
h3

f xi , y xi  
f  , y  
2
6
 in xi , xi 1  ...2
• Using the chain rule f xi , yxi  can be written
as, f xi , yxi  f x xi , yxi   f xi , yxi  f y xi , yxi 
• Using the chain rule f xi , yxi  can be written
as, f xi , yxi  f x xi , yxi   f xi , yxi  f y xi , yxi 
• Equating like terms in equation 1 and 2,
Power
of h
Expansion of y(x)
Runge-Kutta
y  xi 
0
yi
1
f xi , yxi 
2
1
f x  xi , y xi   f  xi , y xi  f y xi , y xi i 
2

a  b f xi , yi 

bpfx xi , yi   bqf xi , yi  f y xi , yi 
• Assuming that yi  yxi  we compare
coefficients. So that, a  b  1, bp  1 2 and bq  1 2
• Assuming that yi  yxi  we compare
coefficients. So that, a  b  1, bp  1 2 and bq  1 2
1
• Therefore, a  1  b and p  q 
2b
• Assuming that yi  yxi  we compare
coefficients. So that, a  b  1, bp  1 2 and bq  1 2
1
• Therefore, a  1  b and p  q 
2b
• However we have four unknowns and only
three equations.
• We have the variable b which can be
chosen arbitrarily.
• Two common choices are b  12 and b  1
• Forb  12 . a  12 , p  1, q  1 and we get,
yi 1  yi 
h
 f xi , yi   f xi  h, yi  hf xi , yi 
2
• which can be written as,
yi 1  yi 
h
 f xi , yi   f xi 1 , yi 1 
2
yi 1  yi  hf xi , yi 
• where
• This improved Euler method or Heun’s
Method.
• Geometric Interpretation
f xi 1 , yi 1 
yi 1
yi 1
yi
h
2
h
xi
xi 1
• For b  1 .a  0, p  12 , q  12 and we get,
h


yi 1  yi  hf  xi  , yi  1 
2
2


yi  1  yi 
h
f  xi , yi 
2
• where
• This improved polygon method or the
modified Euler’s Method.
2
• Geometric Interpretation

f xi  h2 , yi  i
2

yi 1
yi  i
2
yi
h
2
xi
h
2
xi 1
Higher Orders
• The higher orders are developed in the
way as the 2nd order Runge-Kutta.
• Consider the 3rd order Runge-Kutta,
  ak1  bk2  ck3
• k1 , k2 and k3 approximate the derivative at
the points on the interval.
• The higher orders are developed in the
way as the 2nd order Runge-Kutta.
• Consider the 3rd order Runge-Kutta,
  ak1  bk2  ck3
• k1 , k2 and k3 approximate the derivative at
the points on the interval.
• For this case,k1  f xi , yi 
k2  f xi  ph, yi  phk1 
k3  f xi  rh, yi  shk2  r  s hk1 
• Third order Runge-Kutta algorithms are of
the form, yi 1  yi  hak1  bk2  ck3 
• Third order Runge-Kutta algorithms are of
the form, yi 1  yi  hak1  bk2  ck3 
• To determine a,b,c,p,r and s,
First expand k2 and k3 about xi , yi 
• Third order Runge-Kutta algorithms are of
the form, yi 1  yi  hak1  bk2  ck3 
• To determine a,b,c,p,r and s,
First expand k2 and k3 about xi , yi 
Expand yx  as a Taylor series.
• Third order Runge-Kutta algorithms are of
the form, yi 1  yi  hak1  bk2  ck3 
• To determine a,b,c,p,r and s,
First expand k2 and k3 about xi , yi 
Expand yx  as a Taylor series.
Compare coefficients of the Runge-Kutta
and Taylor series.
• The equations formed are,
a b  c 1
bp  cr  1 2
bp  cr  13
2
2
cps 
1
6
• Two of the constants a,b,c,p,r,s are
arbitrary.
• For one set of constant, selected by Kutta
the third order algorithm is,
h
yi 1  yi  k1  4k 2  k3 
6
k1  f xi , yi 
k2  f xi  12 h, yi  12 hk1 
k3  f xi  h, yi  2hk2  hk1 
• For one set of constant, selected by Kutta
the third order algorithm is,
h
yi 1  yi  k1  4k 2  k3 
6
k1  f xi , yi 
k2  f xi  12 h, yi  12 hk1 
k3  f xi  h, yi  2hk2  hk1 
• If f x, y  is a function of x only then this
reduces to the Simpson’s rule.
• Runge-Kutta formulas of higher orders can
be developed using a similar procedure.