Transcript Document

INITIAL VALUE PROBLEMS
REMEMBER
NUMERICAL DIFFERENTIATION
The derivative of y(x) at x0 is:
y  x0   lim
y  x0  h  y  x0 
h 0
h
An approximation to this is:
y  x0  
y  x0  h  y  x0 
h
for small values of h.
Forward Difference
Formula
y  x0  
y  x0  h   y ( x 0 )
h
hy  x0   y  x0  h  y( x0 )
y( x0 )  hy  x0   y  x0  h
y  x0  h  y( x0 )  hy  x0 
Provided h is small.
y  x0  h  y( x0 )  hy  x0 
y  x0  h 
y   x0 
y  x0 
error
hy  x0 
h
x0
x1
y  x0  h  y( x0 )  hy  x0 
The value of a function at
x0  h
can be approximated if we know its value
and its slope at an earlier point,
provided h is small.
x0
A sequence can be written as:
y1  y0  hy0
y1  y0  hy0
y2  y1  hy1
yn1  yn  hyn
In an initial value problem,
and often written as:
Therefore,
yn
is given
f xn , yn 
yn1  yn  hf xn , yn 
Euler’s Method
yn1  yn  hf xn , yn 
yn  y( xn )
Value of y at
yn1  y( xn1 )
x  xn
Value of y at
x  xn1
f xn , yn   yxn 
First derivative of y at
x  xn
h is a small number by which x is incremented.
x1  x0  h
x2  x1  h  x0  2h
x3  x2  h  x0  3h
xn  x0  nh
So what are we doing in Euler’s method?
We are approximating the value of the
function y at x = xn+1 based on the
information that we obtained at the
previous point, x = xn .
The previous information are:
The value of y at x = xn
The value of the first derivative of y at x = xn
yn1  yn  hf xn , yn 
Consider the following initial value problem.
y  y  x
, y(0)  0.5
Find the value of y(0.8).
Let h  0.1
Given:
f xn , yn   y(at x  xn )  yn  xn
y( at x  x0 )  y0  0.5 Our starting point.
where
x0  0
x1  x0  h  0  0.1  0.1
x2  x0  2h  0  2  0.1  0.2
x3  x0  3h  0  3  0.1  0.3
Based on the sequence generated by Euler’s
method we can write:
y1  y0  hf x0 , y0 
f ( x0 , y0 )  y0  x0  0.5  0  0.5
y1  0.5  0.1  0.5  0.55
y2  y1  hf x1 , y1 
x1  x0  h  0  0.1  0.1
f ( x1 , y1 )  y1  x1  0.55  0.1  0.45
y2  0.55  0.1  0.45  0.595
y3  y2  hf x2 , y2 
x2  x0  2h  0  2  0.1  0.2
f ( x2 , y2 )  y2  x2  0.595 0.2  0.395
y3  0.595 0.1  0.395 0.6345
Approximate value of y at x = x3 is 0.6345.
Therefore, we can write that
y( 0.1)  y1  0.55
y( 0.2 )  y2  0.595
y( 0.3)  y3  0.6345
.
.
y( 0.8)  y8  0.728205
Modified Euler’s Method (Heun’s Method)
In the Modified Euler’s method an average value
of the slope is used.



1

*
yn1  yn  h f xn , yn   f xn1 , yn1 
2

Where
*
n1
y
*
n1 is calculated using the Euler’s method.
y
 yn  hf xn , yn 
*
n1
y
*
n1
y
 yn  hf xn , yn 
as calculated using the Euler’s
method shown above is used to
obtain the 2nd slope.
f xn1 , y
*
n1




1

*
yn1  yn  h f xn , yn   f xn1 , yn1 
2

Consider the previous initial value
problem.
y  y  x
, y(0)  0.5
Find the value of y(0.8).
Again,le t h  0.1
f ( x0 , y0 )  y0  x0  0.5  0  0.5
y  0.5  0.1  0.5  0.55
*
1
f ( x1, y )  y  x1  0.55  0.1  0.45
*
1
*
1



1
* 
y1  y0  h f x0 , y0   f x1 , y1 
2




1
* 
y1  y0  h f x0 , y0   f x1 , y1 
2

1

y1  0.5  0.1 0.5  0.45  0.5475
2

f ( x1 , y1 )  y1  x1  0.5475 0.1
 0.4475
y  0.5475 0.1  0.4475
*
2
 0.59225
f ( x2 , y )  y  x2  0.59225 0.2
*
2
*
2
 0.39225



1
* 
y2  y1  h f x1 , y1   f x2 , y2 
2

1

y2  0.5475 0.1 0.4475 0.39225
2

 0.589487
Taylor’s Method of Order Two
2
h
y( x  h)  y( x )  hy( x ) 
y( x )  R2
2!
R2 is called Taylor’s remainder.
If h is small then:
2
h
y( x  h)  y( x )  hy( x ) 
y( x )
2!
2
h
y( x  h)  y( x )  hy( x ) 
y( x )
2!
Taylor’s method of order two can be written
in the form of the following sequence.
2
h
yn1  yn  hyn 
yn
2
Solve the following IVP by using Taylor’s method of
order two.
y   y  x  1 , 0  x  1 , y(0)  1
Solution:
y   y  1  (  y  x  1)  1  y  x
2
h
yn1  yn  hyn 
yn
2
yn   yn  xn  1
yn  yn  xn
2
h
yn1  yn  hyn 
yn
2
2
h
yn1  yn  h(  yn  xn  1)  ( yn  xn )
2

h
  1  h 
2

2
yn1


h
 yn   h 
2


2

 xn  h

Let h = 0.1
yn1  0.905yn  0.095xn  0.1
Results with Taylor’s Method of Order Two
n
xn
yn+1
0
0.0
1.005
1
0.1
1.019025
2
0.2
1.041218
3
0.3
1.070802
4
0.4
1.107076
Interpolation produces a
function that matches the given
data exactly. The function then
should provide a good
approximation to the data
values at intermediate points.
Interpolation may also be used
to produce a smooth graph of a
function for which values are
known only at discrete points,
either from measurements or
calculations.
Given data points
Obtain a function, P(x)
P(x) goes through the data points
Use P(x)
To estimate values at intermediate
points
y
P(x)
8
?
3
2
4
5
x
Assume that a function goes through three points:
 x , y  x   ,  x , y  x   and  x , y  x   .
0
0
1
1
2
2
y( x )  P ( x )
P  x   L0  x  y  x0   L1  x  y  x1   L2  x  y  x2 
Lagrange Interpolating Polynomial
P  x   L0  x  y  x0   L1  x  y  x1   L2  x  y  x2 
( x  x1 )( x  x2 )
P  x 
y  x0 
( x0  x1 )( x0  x2 )
( x  x0 )( x  x2 )

y  x1 
( x1  x0 )( x1  x2 )
( x  x0 )( x  x1 )

y  x2 
( x2  x0 )( x2  x1 )
y( x )  P ( x )
2 x  x1  x2
P  x  
y  x0 
( x0  x1 )( x0  x2 )
2 x  x0  x 2

y  x1 
( x1  x0 )( x1  x2 )
2 x  x0  x1

y  x2 
( x2  x0 )( x2  x1 )
NUMERICAL INTEGRATION
b

a
y( x )dx  area under the curve f(x) between
x  a to x  b.
In many cases a mathematical expression for y(x) is
unknown and in some cases even if y(x) is known its
complex form makes it difficult to perform the integration.
Simpson’s Rule:

x2
x0
y  x  dx   P  x  dx
x2
x0
x1  x0  h and x2  x0  2h

x2
x0
P  x  dx  


x2
x0
x2
x0
x2
x0
( x  x1 )( x  x2 )
y  x0  dx
( x0  x1 )( x0  x2 )
( x  x0 )( x  x2 )
y  x1  dx
( x1  x0 )( x1  x2 )
( x  x0 )( x  x1 )
y  x2  dx
( x2  x0 )( x2  x1 )

x2
x0
y  x  dx   P  x  dx
x2
x0
h
  y  x0   4 y  x1   y  x2  
3
Runge-Kutta Method of Order
Four
Runge-Kutta method uses a
sampling of slopes through an
interval and takes a weighted average
slope to calculate the end point. By
using fundamental theorem as shown
in Figure 1 we can write:
y  x0  h  y( x0 )  hy  x0 
y  x0  h 
y   x0 
y ( x0 )
error
hy  x0 
h
x0
x0  h
y  x0  h 
y   x0 
y ( x0 )
error
hy  x0 
h
x0
x0  h
dy  y x dx
 (1)
y  x n  h
y  xn 
dy
dx
xn
xn  h
Integrating both sides of Eqn. (1) we get

dy

y
   x dx
 (2)
Applying appropriate limits to Eqn. (2)
we get
y xn1   y xn   
xn  1
xn
y xn  h  y xn   
y x dx
xn  h
xn
y x dx
 (3)
 (4)
Let us now concentrate on the right-hand
side of Eqn. (4).

xn  h
xn
y x  dx  
xn  h
xn
P  x  dx  (5)
P(x) can be generated by utilizing Lagrange
Interpolating Polynomial. Assume that the
only information we have about a function,
f(x) is that it goes through three points:
 x , y  x   ,  x , y  x   and  x , y  x   .
0
0
1
1
2
2
P ( x )  y( x )
P  x   L0  x  y  x0   L1  x  y  x1   L2  x  y  x2 
( x  x1 )( x  x2 )
P  x 
y  x0 
( x0  x1 )( x0  x2 )
( x  x0 )( x  x2 )

y  x1 
( x1  x0 )( x1  x2 )
( x  x0 )( x  x1 )

y  x2 
( x2  x0 )( x2  x1 )
Using Simpson’s Integration,

x2
x0
y  x  dx   P  x  dx
x2
x0
x1  x0  h and
x2  x0  2h

x2
x0
P  x  dx  



x2
x0
x2
x0
x2
x0
x2
x0
( x  x1 )( x  x2 )
y  x0  dx
( x0  x1 )( x0  x2 )
( x  x0 )( x  x2 )
y  x1  dx
( x1  x0 )( x1  x2 )
( x  x0 )( x  x1 )
y  x2  dx
( x2  x0 )( x2  x1 )
y  x  dx  
x2
x0
h
P  x  dx   y  x0   4 y  x1   y  x2  
3
...………………….(6)
We can apply Simpson’s Integration to Eqn.
(6) with the following substitutions:
x0  x n
h  h / 2

xn  h
xn
;
x2  xn  h
;
and the midpoint,
h
x1  xn 
2

h
h

y x dx   y xn   4 y xn    y xn  h
6
2


….(7)
Compared to Euler’s Formula, an average of
six slopes is used in Eqn. (7) instead of just
one slope. Actually, the slope at the midpoint
has a weight of 4. The slope at the midpoint
can be estimated in two ways.

xn  h
xn
y x dx 

h
h
h






y  xn   2 y  x n    2 y  x n    y  x n  h

6
2
2



y xn  h  y xn  

h
h
h










y
x

2
y
x


2
y
x


y
x

h
 n

 n

n
n


6
2
2





y xn  h  y xn  

h
h
h










y
x

2
y
x


2
y
x


y
x

h
 n

 n

n
n


6
2
2






h
h
h


y xn1   y xn    y xn   2 y xn    2 y xn    y xn1 
6
2
2



........... (8)
The slopes can be estimated in the
following manner:
y xn   k1  f  xn , yn 
h
h
h 


y xn    k2  f  xn  , yn  k1 
2
2
2 


h
h
h 


y xn    k3  f  xn  , yn  k2 
2
2
2 


y xn  h  k4  f  xn  h, yn  hk3 
Substituting the estimated slopes into Eqn.
(8) gives the formula for Runge-Kutta
Method of Fourth Order:
yn1
h
 yn  k1  2k2  2k3  k4 
6