A short note on Numerical Analysis II meant for self-taught

students until face to face class will resume

Complied by

Hailu Muleta (Assistant Professor of Mathematics)

April, 2020
Jimma, Ethiopia
 Chapter One

Revision on Numerical Integration

This chapter is highly devoted to find the numerical integration for a

given set of data points. The topics covered here are

- Numerical Integration

- Newton-Cote’s quadrature formula

- Trapezoidal Rule

- Simpson’s one-third Rule

- Simpson’s three-eighths Rule

- Weddle’s rule

Consider the definite integral a f  x dx .

b

This integral represents the area between y  f (x) , the x –axis and the lines x  a & x  b .

This integration is possible as far as f (x ) is explicitly given and the function is integrable.

Now suppose set of (n+ 1) paired values are given. First as we did in the case of numerical

 p x dx
b
differentiation, we find f (x ) by an interpolating polynomial Pn (x) and obtain a
n which

 p x dx .
b
can approximate the value for a
n

A General Quadrature Formula for Equidistant Spacing (Newton-Cote’s Formula)

For equally spaced intervals, we have Newton’s forward difference formula as

u u  1 2 u u  1u  2  3
y  x   y0  uy0   y0   y0  ...1
2! 3!

2
 x  x0
Here, u , h is the interval of differencing.
h

Now, instead of f (x ) , we will replace it by this interpolating polynomial y (x ) of Newton.

x  x0 x  x0
Since xn  x0  nh and u  un
h h
x 0  nh
 f x dx   f  x dx
xn
Thus, x0 x0

x 0  nh
 Pn x dx . Where Pn (x) is interpolating polynomial of degree n.
x0

n uu  1 2 uu  1u  2 3 

   y0  uy0   y0   y0  ... hdu
0
 2! 3! 

 n u u  1 2 u 3  3u 2  2u 3 
 h   y0  uy0 
  y0   y0  ... du
0
 2! 3! 

n
 u3 u 2 u4 
 2   u3  u 2 
u
 h uy0  y0  3 2 2 y0  4 3 y0  ...
 2 2! 3! 
  0

 n2 1  n3 n 2  1  n4 
 h ny0  y0    2 y0    n3  n 2 3 y0  ...(2)
 2 2 3 2 6 4 

This equation is called Newton-cote’s quadrature formula and is a general quadrature

formula.

By taking different values for n we get a number of special formulas. Here we try to look for

some values of n for which their practical application is very important in different disciplines of

science. Detail information and results are explained below.

3
 Trapezoidal Rule

Put n = 1, in the quadrature formula

x0  nh x0  h  
f x dx   f x dx  h 1. y0  y0 
1
x0 x0
 2 

Since other differences do not exist if n = 1.

 
 h y0   y1  y0 
1
 2 


h
 y0  y1 
2

x 0  nh
 f x dx   f x dx
xn

x0 x0

x 0  nh x0  2 h x 0  nh
 f x dx   f x dx  ...   f x dx
x0 x0  h x 0   n 1h


h
 y0  y1   h  y1  y2   ...  h  yn 1  yn 
2 2 2


h
 y0  yn   2 y1  y2  ...  yn 1 
2

This is known as a trapezoidal rule.

Even though this method is very simple for calculation, the error in this case is significant.

Truncation Error in Trapezoidal Rule

In the neighborhood of x  x0 , we can expand y  f (x) by Taylor series in powers of x  x0 .

That is,

yx   y0  x  x0  y 
 x  x0  ''
2
y0  ...
'
1
0
2!

4
 x1 
ydx    y0   x  x0  y0 
x1  x  x0  ''
2

 y0  ...dx
'

 2! 
x0 x0

 x  x0  ' x  x0  ''
x1
 2 3

  y0 x  y0  y0  ...
 2! 3!  x0

 y0 x1  x0  
x1  x0 2 y '  x1  x0 3 y ''  ...
0 0
2! 3!

h 2 ' h3 ''
 hy0  y0  y0  ... (2)
2! 3!

If h is the equidistant length, then also

h
 y0  y1   Area of the first trapezium = A0 say
x1

x0
ydx 
2

Putting x  x1 in (1), we get

yx1   y1  y0  x1  x0  y 
 x1  x0  ''
2
y0  ...
'
0
2!

h 2 ''
 y0  hy0'  y0  ...
2!

h h 2 '' 
A0   y0  y1    y0  y0  hy0 
h '
y0  ... 
2 2 2! 

h2 ' h3 ''
 hy0  y0  y0  ...
2 2.2!

Subtracting A0 from (2), we obtain

3 ''  1 1 
x1
x0 ydx  A0  h y0     ...
 3! 2.2! 

5
 1
  h3 y0''  ...
2

1 3 ''
 The error made in the first interval ( x0 , x1 ) is  h y0  ...
2

1
Similarly the error in the ith interval   h yi 1
3 ''

Hence, the total cumulative error E is

1

E   h3 y0''  y1''  y2''  ...  yn'' 1
2


 E 
nh3
12
M 
where M  max y0'' , y1'' , y 2'' ,... 

b  a h 2 M ba
if the interval is (a, b) and h 
12 n

 The error in the trapezoidal rule is of order h 2 .

The accuracy of the result can be improved by increasing the number of intervals and

decreasing the value of h .

Simpson’s One-Third Rule

Setting n = 2 in Newton-cotes quadrature formula, we have

 18 4 
 f xdx  h2 y
x2
0  4y0    2 y0  since the other terms vanish (become zero).
x 23 2 

 
 h2 y0  2 y1  y0   E  1 y0 
1 2

 3 

 
 h2 y0  2 y1  2 y0   y2  2 y1  y0 
1
 3 

6
 1 4 1 
 h y2  y1  y0 
3 3 3 


h
 y2  4 y1  y0 
3

 f x dx  3  y  4 y3  y4  and f  x dx   yi  4 yi 1  yi  2 

x4 h xi  2 h
Similarly,
x2
2 
xi 3

If n is an even integer, then the last integral will be

f  x dx 
h
 y n2  4 y n1  y n 
xn
 xn  2 3

Adding all these integrals, if n is even positive integer, then y 0 , y1 , y 2 , y3 ,..., y n are odd in

number; we have

x0 f ( x)dx  x0 f ( x)dx  x24 f x dx  ...  xn2 f x dx

xn x2 x xn


h
 y0  4 y1  y2    y2  4 y3  y4   ...   yn  2  4 yn 1  yn 
3


h
 y0  yn   2 y2  y4  ...  yn  2   4 y1  y3  ...  yn 1 
3

Simpson’s Three-Eighths Rule

Putting n = 3 in Newton-cotes formula, we get

x  9 19 1  81  
x03 f ( x)dx  h3 y 0  y0   2 y0    27  9 3 y0 
2 2 2 6 4  

 
 h3 y0   y1  y0   E  1 y0  E  1 y0 
9 9 2 3 3

 2 4 8 

 
 h3 y0  y1  y0   y2  2 y1  y0    y3  3 y2  3 y1  y0 
9 9 9 3
 2 2 4 8 


3h
 y3  3 y2  3 y1  y0 
8

7
If n is a multiple of 3,

x  nh x 3h x 6h x  nh
x00 f ( x)dx   0 f ( x)dx   0 f x dx  ...   0 f x dx
x0 x 3h x0 ( n3) h
0


3h
 y0  3 y1  3 y2  y3    y3  3 y4  3 y5  y6   ...   yn 3  3 yn  2  3 yn 1  yn 
8


3h
 y0  yn   3 y1  y2  y4  y5  ...  yn  2  yn 1   2 y3  y6  yn 
8

This is Simpson’s three-eighths rule and is applicable only when n is a multiple of 3.

Weddle’s Rule

Putting n = 6 in Newton-cotes formula

x0 6h  
f ( x)dx  h 6 y0  18y0  72  182 y 0  324  216  363 y0  ...
1 1
x0  2 6 

 123 4 33 41 6 
 h6 y0  18y0  272 y0  243 y0   y0  5 y0   y0 
 10 10 140 

41 6 42 6 h 6
Now replace the term  y0 by  y0 , by doing this, the error introduced is only  y0
140 140 140

which is negligible when h and 6 y0 are small.

Using   E  1 and replacing all differences in terms of y ’s,

we get

x0 6h

x0
f x dx 
3h
 y0  5 y1  y 2  6 y3  y 4  5 y5  y6 
10

x 12h
x006h f x dx  10 y  5 y 7  y8  6 y 9  y10  5 y11  y12 
3h
Similarly, 6

x  nh
x00(n6)h f x dx  10 y  5 y n 5  y n  4  6 y n 3  y n  2  5 y n 1  y n 
3h
and n 6

Adding all these integrals, we get

8
 x  nh
x00 f x dx 
3h
 y0  5 y1  y 2  6 y3  y 4  5 y5   2 y6  5 y7  y8  6 y9  y0  5 y11 
10
 ...  2 y n 6  5 y n 5  y n 4  6 y n 3  y n 2  5 y n 1  y n 

This equation is called Weddle’s rule.

Truncation Error in Simpson’s Formula

By Taylor expansion of y  f (x) in the neighborhood of x  x0 , we obtain

y  y0  x  x0  y0' 
x  x0 2 y ''  x  x0 3 y '''  ... ...1
0 0
2! 3!

x x   x  x0  ' x  x0  '' x  x0  '''

2 3

 2 ydx  2 y0
x0 x0    y  0 y  y  ... dx
0 0 
 1! 2! 3! 

 x  x0  ' x  x0  '' x  x0  ''' 

x2
 2 3 4
  y0 x  y0  y0  y0  ..
 21! 3! 4!  x0

 y0 x2  x0  
x2  x0 2 y '  x2  x0 3 y ''  x2  x0 4 y '''  ...
0 0 0
2! 3! 4!

4h 2 ' 8h3 '' 16 2 '''

 2hy0  y0  y0  h y0  ...
2! 3! 4!

4 2 h 4 ' ' ' 4 h 5 4 

 2hy0  2h 2 y0'  h3 y0''  y0  y0  ...
3 3 15

h
 y0  4 y1  y 2  , by Simpson’s rule
x2
Now let A1  area  x ydx 
0 3

Putting x  x1 in (1), to get

y1  y0   x1  x0  y0' 
x1  x2  y ''  ...
0
2!

h 2 '' h3 '''
 y0  hy0'  y0  y0  ...
2! 3!

Putting x  x2 in (1), we have

9
 4h 2 '' 8h3 '''
y2  y0  2hy  '
0 y0  y0  ...
2! 3!

h   
y0  y0  ...    y0  2hy0'  2h 2 y0''  ...
h 2 '' h3 '''
 A1   0
y  4
 y0  hy0
'

3  2! 3!  

4 3 ' ' 2 4 ' ' ' 5 5 4 

 2hy0  2h 2 y0'  h y0  h y0  h y0  ...
3 3 18

4 5
  ydx  A1    h5 y04  ...
x2

x0
 15 18 

 h 5 4 
 y0  ... Omitting the remaining terms involving h 6 and higher powers of h .
90

 h 5 4 
This means that the error made in x0 , x2  is  y0  ...
90

 h 5 4 
Similarly, the error made in x2 , x4  is y2 and so on.
90

Hence the total error E is x0 , xn  is  

9

h 5 4 
y0  y24   ... 

nh5
E  M , where M is the maximum value of y04, y24  ..., y24n 2
90

Since x2 n , y2 n  is the last paired value because we require odd number of ordinates to apply

Simpson’s one-third rule.

If the interval is (a, b) , then b  a  h(2n) , using this

E 
b  a h 4
M
180

Hence, the error in Simpson’s one-third rule is of the order h 4 .

Examples

10
 7
1. Evaluate 
1
x 2 dx by using

a) Trapezoidal rule

b) Simpson’s rule and verify your results by actual integration

Solution

Here y ( x)  x 2 and the interval length is 7-1 = 6 so we divide this interval into 6 equal parts

6
with h  1
6

x: 1 2 3 4 5 6 7

y: 1 4 9 16 25 36 49

a) By Trapezoidal rule

h
 y1  y7   2 y2  y3  y4  y5  y6 
7
1
x 2 dx 
2

1
1  49  24  9  16  25  36  115
2

b) By Simpson’s one-third rule

1
 y1  y7   2 y3  y5 )  4( y2  y4  y6 
7
1
x 2 dx 
3


1
1  49  29  25  44  16  36
3

= 114

Since n = 6, we can also use Simpson’s three-eighth rule.

3
 y1  y7   3 y2  y3  y5  y6   2 y4 
7
So 
1
x 2 dx 
8


3
1  49  34  9  25  36  216
8

11
 = 114

7
7 x3 1 1 342
 x dx   73    114
2
1 3 1 3 3 3

dx
h  0.2 , and hence
1
2. Evaluate  1 x
0 2
, using Trapezoidal rule with

obtain an approximate value of  .

Solution

yx  
1
Let
1  x2

x: 0 0.2 0.4 0.6 0.8 1

y: 1 0.961538461 0.862068965 0.735294117 0.609756097 0.5

dx h
 y0  y5   2 y1  y2  y3  y4 
1
 1 x
0 2

2


0.2
1  0.5  20.961538461 0.862068965 0.735294117 0.6097566097
2

= 0.783731528

But by actual integration

1 dx 
 1 x  tan 1 x 0 
1
2
0 4


  0.783731528
4

  3.134926112

To compare this approximated value of  with its actual value using calculator 3.141592654 ,

the error is 0.00666654159 which is 0.6% .

12
3. From the table below, find the area bounded by the curve and

the x -axis from x  7.47 to x  7.53

x: 7.47 7.48 7.49 7.50 7.51 7.52 7.53

y: 1.93 1.95 1.98 2.01 2.03 2.06 2.08

Solution

Let us compare the results obtained by different methods

i) By Trapezoidal rule

f  x dx 
0.01
1.93  2.08  21.95  1.98  2.01  2.03  2.06
7.53
7.47 2

= 0.12035

ii) By Simpson’s one-third rule

f  x dx 
0.01
1.93  2.08  21.98  2.03  4  41.95  2.01  2.06
7.53
7.47 3

 0.120366667

iii) By Simpson’s three-eighths rule

f  x dx 
3(0.01)
1.93  2.08  31.95  1.98  2.03  2.06  22.01
7.53
7.47 8

 0.1203375

iv) By Weddle’s rule

f  x dx 
3(0.01)
1.93  51.951.98  62.01  2.03  5(2.06)  2.08
7.53
 7.47 10

 0.12039

As we can see from these rules the area is 0.1203 (correct to four decimal places)

13
 5.2
4. Evaluate the integral I  4 nxdx using the rules so far developed.

Solution

1 .2
Since b  a  5.2  4  1.2 let us divide the interval into 6 equal parts, i.e.; h   0 .2
6

x ln x x ln x

4 1.386294361 4.6 1.526056303

4.2 1.435084525 4.8 1.568615918
4.4 1.481604541 5.0 1.609437912
5.2 1.648658626

i) By Trapezoidal rule


0.2
1.386294361 1.648658626  2(1.435084525 1.481604541
5.2
nxdx 
4 2
1.526056303 1.568615918 1.609437912)]

= 1.827655139

ii) By Simpson’s one-third rule


0.2
1.386294361 1.648658626  2(1.481604541 1.568615918)  4
5.2
nxdx 
4 3
(1.435084525 1.526056303 1.609437912)]

 1.827847258

iii) By Simpson’s three-eighths rule


0.2
1.386294361 1.648658626  3(1.435084525 1.481604541)  3
5.2
nxdx 
4 3
(1.568615918 1.609437912)  2(1.526056303)

 1.827847258

iv) By Weddle’s Rule

14
 3(0.2)
1.386294361 51.435084525  1.481604541 61.526056303 
5.2

4
nxdx 
10
6 (1.526056303)  1.1.568615918 51.609437912  1.648658626

 1.827847407

By actual integration,

nxdx  xnx  1
5.2
  1.827847409
5.2
4 4

Here Weddle’s rule best approximates the exact value.

1
 e dx by Simpson’s one-third rule correct to five decimal places.
x
5. Evaluate 0

Solution

The interval b  a  1

Since error E 
b  a  h 4 M , where M = Max (e x ) in the range
180

1 4
 h e
180
6
Now we require E  10

h 4e
  10 6
180

1
 180  10 6  4
 h     0.090207886 0.1
 e 

Hence we take h  0.1

15
 1
  e x dx 
0
0.1
3
   
1  e  2 e0.2  e0.4  e0.6  e0.8  4 e0.1  e0.3  e0.5  e0.7  e0.9 
 1.718282782

By the actual integration,

1

1
e x dx  e x  e  1  1.718281828
0 0

1
 e dx  1.71828 which is the same as the exact value.
x
Correct to five decimal places 0

6. A curve passes through the points (1,2), (1.5,2.4),(2.0,2.7),(2.5,2.8), (3,3), (3.5,2.6)

and (4,2.1) . Obtain the area bounded by the curve, the x -axis, x  1 and x  4 .

b 4
Area =  ydx  
a 1
ydx


0.5
2  2.1  22.7  3  4(2.4  2.8  2.6)
3


1
4.1  11.4  31.2  7.7833 sq. units
6

b 4
Volume =  a y dx   1 y dx
2 2


0.5 2
3

2 2
  2 2
 
2  2.1  2 2.7   3 2  4 2.4   2.8  2.6 
2


 8.41  32.58  81.44  20.405 = 64.1041981 cubic units.
6

7. A river is 80 meters wide. The depth ‘d’ in meters at a distance x

meters from one bank is given by the following table.

x: 0 10 20 30 40 50 60 70 80

d: 0 4 7 9 12 25 14 8 3

Calculate the area of cross-section of the river using Simpson’s rule.

16
 Solution

80
Area of cross-section = 
0
ydx


10
0  3  2(7  12  14)  4(4  9  15  8) = 710 sq. meters
3

8. The table below gives the velocity v of a moving particle at time t second. Find the

distance (S) covered by the particle in 12 seconds and also the acceleration at t  2

seconds.

t: 0 2 4 6 8 10 12

v: 4 6 16 34 60 94 136

Solution

ds dv
We know that v  and a 
dt dt

To get S

2
4  136  216  60  46  34  94
12
S  
0
vdt 
3

= 552 meters

 dv 
To find a, a    first form the difference table
 dt  t 2

t v v 2 v 3v

0 4
2
8
2 6
10 0
8
4 16
18 0
8
6 34
26
0
8 60 8
34
0
17
10 94
42 8
12 136

 dv  1 1 2 1 3 
   v0   v0   v0 
 dt t  2 h  2 3 

Taking v 0 = 6

1 1 
 10  8  3m / sec2
2 2 

Exercise

2
dx
1. Evaluate 1 x
1
2
taking h  0.2 , using Trapezoidal rule. Can you use

Simpson’s rule? Justify your reasons.

1
2. Compute the value of 
0
sin x  cos x dx correct to four decimal places

taking h  0.8 .

1 1
x2
3. Find the value of log 2 from 3
0 1  x 3 dx using Simpson’s one-third rule

with h  0.25 .

4. When a train is moving at 30m / sec. steam is shut off and brakes are applied. The speed of the

train per second after t seconds is given by

Time t : 0 5 10 15 20 25 30 35 40

Speed v : 30 24 19.5 16 13.6 11.7 10 8.5 7.0

Using Simpson’s rule, determine the distance moved by the train in 40 seconds.

1
5. Evaluate  e  x dx a) dividing the range into four equal parts
2

18
 b) dividing the range into ten equal parts by

i) Trapezoidal rule and

ii) Simpson’s one-third rule


2

6. Evaluate  e dx taking h 
sin x
.
0 6



7. Calculate  sin 3 xdx taking h  .
0
6

x
2
8. Evaluate log xdx taking four strips.
3

0 .7

e
x
9. Calculate x dx taking 5 ordinates by Simpson’s rule.
0 .5

0 .5
dx
10. Evaluate 
0 1 x2
dx by Weddle’s rule, dividing the range into six

parts.


sin x
11. Evaluate 
0
x
dx dividing into six equal parts using Simpson’s rule,

Weddle’s rule and Trapezoidal rule.

19
 Chapter Two

Curve Fitting

Fitting of curves to a set of numerical data is of considerable importance- theoretical as well as

practical. Theoretically it is useful in the study of correlation and regression. In practice it
enables us to represent the relationship between two variables by simple algebraic expressions
(polynomials, exponential or logarithmic functions or any). Besides, it may be used to estimate
the values of one variable which would correspond to the specified values of the other
variable(s).

This chapter covers how to fit a curve for a given set of data points using different methods and
it focuses on the following points:

Regression

- Linear regression

- quadratic regression

- polynomial regression

- multiple regression

- fitting an exponential curve

- curve fitting with Sinusoidal Functions

In most of the fields of engineering and science, we come across experiments which

involve many variables, and most of the time data is collected or given for discrete values

along a continuum; the relation between these variables can be discussed so easily and for

many of these variables it is very difficult to identify the relation unless we can model the

system mathematically. When the system is explained in terms of mathematical models we

have the following relationships about the variables:

1. The relationship between these variables is given in terms of

20
 mathematical rules, formulae if any, to determine the quantities of these variables. Actually

it is simple to use these rules for application.

2. The quantities/ variables are given so that we will be interested in finding the

relationships between these variables. This process is a little bit difficult because to write

one variable in terms of the other variables (called empirical equation). Most of the time

we may not be able to get an exact relation between these variables and we may get only

an approximate relation or curve.

This approximating curve is an empirical equation and the method of finding such an

approximating curve is called curve fitting.

Suppose xi , yi  , i  1,2,3,..., n be n sets of observations and the law relating x and

y can be determined by different mathematical systems that clearly explains the relationship

between these sets of n observations xi , yi  . Actually, here we may have different approaches to

fit the given data, and one system may approximate better than the other system on the same given

set of data points.

Now we will see some of these different approaches:

REGRESSION

1. LINEAR REGRESSION

Suppose that the relationship is given

y i  a  bxi , i  1,2,3,..., n . (1)

Equation (1) represents a family of straight lines for different values of the arbitrary constants ' a '

and ' b '. The problem now is to determine ' a ' and ' b ' so that the line (1) is the line of "best fit".

The term best fit is interpreted in accordance with the Legendre's principle of least squares which

consists of the deviations of the actual values as given by the line of best fit. As a matter of

21
chance all the points may lie on a straight line and in this case the line is a 'perfect fit' and the

sum of the squares of the deviations is zero.

Let Pi ( xi , y i ) be any point in the scatter diagram. Draw Pi M perpendicular to the x-axis meeting

the line y  a  bx in H i . The coordinates of H i are xi , a  bxi  .

Pi H i  Pi M  H i M
 yi  (a  bxi )

PiHi is called the error of estimates or the residual for y i .

According to the principle of least squares, we have to determine a and b so that

n n
E   Pi H i   ( yi  a  bxi ) 2 is the minimum.
2

i 1 i 1

Using the principle of maxima and minima what we have studied in calculus, the partial

derivatives of E with respect to a and b should vanish separately.

E n
That is,  2 ( yi  a  bxi )  0
a i 1

n n n
  y i   a  b xi
i 1 i 1 i 1
n
 na  b xi i 
i 1

22
 E n
 2 xi ( y i  a  bxi )  0
b i 1
. n n n
  xi y i  a  xi  b  xi
2
ii 
i 1 i 1 i 1

Equations (i) and (ii) are called normal equations.

Solving for a and b from (i) and (ii), we get the values of a and b , and with these values of a

and b so obtained, equation (1) is the line of best fit to the given set of points  xi , y i  ,
i  1,2,3,..., n .

Now let us see some examples to illustrate the above discussion.

Example 1. By the method of least squares find the best fitting

straight line to the data given below:

x: 5 10 15 20 25

y: 15 19 23 26 30

Solution

Let the line of best be y  a  bx

The normal equations are

n n

 yi  5a  b xi
i 1 i 1

n n n

 xi yi  a xi  b xi
2

i 1 i 1 i 1

We calculate  x,  y,  x ,  xy and form the table below

2

23
x y x2 xy
5 15 25 75
10 19 100 190
15 23 225 345
20 26 400 520
25 30 625 750
75 114 1375 1885
Using these values in the normal equations, we get

5a  75b  114
75a  1375b  1885

Solving for a and b we get a=12.3 and b=0.7 and thus the line of best fit is y  12.3  0.7 x

Example 2. Find the best fitting straight line to the data given below

by the method of least squares and also estimate y when

x is 70.

x : 71 68 73 69 67 65 66 67

y : 69 72 70 70 68 67 68 64

Solution

First transform the values of x and y to X  x  68 and Y  y  70

and the normal equations are

b X  8a   Y
b X 2  a  X   XY

Calculations:

x y X Y X2 XY

71 69 3 -1 9 -3

68 72 0 2 0 0

24
73 70 5 0 25 0

69 70 1 0 1 0

67 68 -1 -2 1 2

65 67 -3 -3 9 9

66 68 -2 -2 4 4

67 64 -1 -6 1 6

2 -12 50 18

Substituting these values in the normal equations, we get

2b  8a  12
50b  2a  18

35 16
Solving for a and b , we get b  and a
99 99

35 16
Thus the line of best fit is of the form Y  X
99 99

35 16
This implies y  70  ( x  68) 
99 99

35 4566
y x
99 99

When x  70  y  70.87

1. Quadratic Regression (Fitting Of Second Degree Parabola)

Let y  a  bx  cx 2 be the second degree parabola of best fit to set of n points xi , yi  ,

i  1,2,3,..., n .

Using the principle of least squares, we have to determine a,b, and c so that

 
n
E   yi  a  bxi  cx 2
i is minimum.
i 1

25
 Equating to zero the partial derivatives of E with respect to a,b, and c separately, we get the

normal equations for estimating a,b, and c as

E
 2 yi  a  bxi  cxi2   0
n

a i 1

n n n
  yi  na  b xi  c xi2 …………………. (1)
i 1 i 1 i 1

E
 2 yi  a  bxi  cxi2 xi  0
n

b i 1

n n n n
  xi yi  a xi  b xi2  c xi3.......... .......... ....( 2)
i 1 i 1 i 1 i 1

E
 2 yi  a  bxi  cxi2 xi2  0
n

c i 1

n n n n
  xi2 yi  a xi2  b xi3  c xi4 .......... .......... ...(3)
i 1 i 1 i 1 i 1

Solving for a,b, and c from (1), (2) and (3), we get with these values of a,b, and c the

parabola of best fit.

Example 1 Fit a parabola of second degree to the following data

X: 0 1 2 3 4

Y: 1 1.8 1.3 2.5 6.3

Solution

X Y X2 X3 X4 XY X2Y

0 1 0 0 0 0 0

1 1.8 1 1 1 1.8 1.8

2 1.3 4 8 16 2.6 5.2

26
 3 2.5 9 27 81 7.5 22.5

4 6.3 16 64 256 25.2 100.8

10 12.9 30 100 354 37.1 130.3

Substituting these values in the normal equations, we get

12.9  5a  10b  30c

37.1  10a  30b  100c
130.3  30a  100b  354c

Solving for a,b, and c , we get

a  1.42, b  1.07 , and c  0.55

 y  1.42  1  07 x  0.55x 2 is the best fit.

Exercise

1. Find the best fitting parabola to the data given below by the

method of least squares and also estimate y when x is 70.

x : 71 68 73 69 67 65 66 67

y : 69 72 70 70 68 67 68 64

2. Polynomial Regression Fitting Of A Polynomial Of kth Degree

If y  a0  a1 x  a2 x 2  ...  ak x k is the kth degree polynomial of best fit to the set of points

xi , yi  ; i  1,2,3,..., n the constants a 0 , a1 , a 2 , ... , a n are to be obtained so that

E   yi  a0  a1 xi  a2 xi2 ...ak xik  is minimum.

n
2

i 1

Thus the normal equations for estimating a 0 , a1 , a 2 , ... , a n are obtained on equating to zero the

partial derivatives of E with respect to a 0 , a1 , a 2 , ... , a n separately.

27
E
 0   yi  na0  a1  xi  a2  xi2  ...  ak  xik
a0

E n
 0   xi yi  a0  xi  a1  xi2  ...  ak  x k 1
a1 i 1

E n
 0   xik yi a0  xik  a1  xik 1  ...  ak  xi2 k
ak i 1

Exercise

1. Find the polynomial of degree three that best fits the data given below by the method of

least squares and also estimate y when x is 70.

x : 71 68 73 69 67 65 66 67

y : 69 72 70 70 68 67 68 64

4. Multiple Regressions

There are different multiple regression forms. For the sake of discussion let us see the following

regression type.

Suppose we want to fit the set of data points by the relation

Z  ax  bxy  cy , here we need the points to be of the form xi , y i , z i  .

We determine the values of a,b, and c , from

n
E   zi  axi  bxi yi  cyi 
2

i 1

E n
 2 zi  axi  bxi yi  cyi xi  0
a i 1

n n
  xi zi  a xi2  b xi2 yi  c xi yi .......... .......... .......... (1)
i 1 i 1

28
 E n
 2 zi  axi  bxi yi  cyi xi yi
b i 1

n n n n
  zi xi yi  a xi2 yi  b xi2 yi2  c xi yi2 .......... .......... .......... (2)
i 1 i 1 i 1 i 1

E n
 2 zi  axi  bxi yi  cyi yi
b i 1

n n n n
  zi yi  a xi yi  b xi yi2  c yi2 .......... .......... .......... (3)
i 1 i 1 i 1 i 1

By solving (1), (2) and (3) for a,b, and c , we get the best approximation.

5. Fitting an Exponential Curve

Let xi , yi  , i  1,2,3,..., n be the n sets of observations of related data and let y  ab be the
x

best fit for the data.

Then taking logarithm on both sides,

y a b
log  log  x log (*)
10 10 10

y a b
Let Y  log10 , A  log10 , and B  log , then (*) reduces to
10

Y  A  Bx which is linear in x and Y , we can find A, B since x and Y are known, and from

A, B we can get a, b and hence y  ab is found out.

x

Fitting a curve of the form y  ax b

y  ax b

y a x
 log  log  b log
10 10 10

y a x
 Y  A  bx, letting y  log , A  log , x  log
10 10 10

Using this linear fit, we find A, b .

29
 a, b are known and thus y  ax b is found out.

Example 1. From the table given below, find the best values of

a and b in the law y  aebx by the method of least squares.

x: 0 5 8 12 20

y : 3.0 1.5 1.0 0.55 0.18

Solution

Let y  aebx be the approximating curve.

e e
 Y  A  bx log  Y  A  Bx where B = b log
10 10

So the normal equations are

B x  5 A   Y
B x 2  A x   xY

x y Y x2 xY

0 3.0 0.4771 0 0

5 1.5 0.1761 25 0.8805

8 1.0 0 64 0

12 0.55 -0.2596 144 -3.1152

20 0.18 -0.7447 400 -14.894

45 -0.3511 633 -17.1287

Substituting these values, we get

5 A  45B  0.3511
45A  633B  17.1287

Solving for A and B , we get

30
 A  0.4815

B  0.0613

So a  10 A  100.4815  3.0304

e
B  b log  0.0613
10

10
 b  (0.0613) log  0.1411
e

0.1411x
Hence the curve is y  3.0304e

Exercise

1. Fit a straight line to the following data and hence find y ( x  25)

x: 0 5 10 15 20

y: 7 11 16 20 26

2. Fit a straight line to the data

x : 0.5 1.0 1.5 2.0 2.5 3.0

y : 0.31 0.82 1.29 1.85 2.51 3.02

3. Fit a parabola to the data

x: 1 2 3 4 5

y: 2 3 5 8 10

4. Fit a curve of the form y  aebx to the data given below:

x: 1 2 3 4 5 6 7 8

y: 15.3 20.5 27.4 36.6 49.1 65.6 87.8 117.6

5. Fit a curve of the form z  ax  by  cxy to the data given below

31
 (0,0,1), (0,1,2), (1,0,4), (1,1,1), (2,0,4), (1,2,5)
a
6. It is given that x , and y are related by y   bx to the data below and obtain the best
x

values of a and b.

x: 1 2 4 6 8

y: 5.43 6.28 10.32 14.86 19.5

32
 Chapter Three

Numerical Solution of Ordinary Differential Equations

In this chapter we are highly interested in finding the solution of ordinary differential

equation numerically using different methods.

The topics covered here is

- Taylor series method

- Taylor series method for simultaneous first order DE

- Taylor series method for second order differential equation

- Picard’s method of successive approximations
- Euler’s Method
- Runge-Kutta method

- Predictor-corrector method

In the fields of Engineering and Science, we come across through the natural phenomena that can

be represented by mathematical models which happen to be in the form of differential equations;

for instance, the equation of motion, the equation of deflection of a beam, etc. The solution of

these differential equations is very essential in the studies of such phenomena.

While finding the solution of these differential equations there are number of differential

equations that we cannot solve analytically; however, in such situations, depending on the nature

of the model, we go for numerical solutions of these differential equations. In many researches,

especially after the advent of modern computers, the numerical solutions of the differential

equations have become very easy for manipulation.

33
Thus, in this part we try to look some of the methods of numerical solutions that are

approximate solutions and in many cases these solutions are in the required (desired) degree of

accuracy and are quite sufficient.

dy
Suppose we want to solve  f ( x, y ) with the initial condition y ( x0 )  y 0 .
dx

Let y ( x0 )  y 0 , y( x1 )  y1 , y( x2 )  y 2 , … be the solution of y at x  x0 , x1 , x2 ,...

Let y  y (x) be the exact solution. If we plot and draw the graph of

y  y (x ) , (the exact curve) and also draw the approximate curve by plotting

( x0 , y 0 ), ( x1 , y1 ), ( x 2 , y 2 ),... (the approximated solution graph) we get two curves.

Approximate solution
Q

P Exact solution
(x0,y0)
QM = approximate value

PM = exact value at x = xi

Then QP = QM-PM
M
= yi – y(xi) = is called
 f  x, y ........(1)
dy
Suppose y ' 
dx

The equation y  f ( x, y ) subject to the initial condition

y ( x0 )  y0 is called an initial-value problem.

Using Taylor series we can expand y (x ) in the neighborhood of x 0 as a power

series of x  x0 . That is, if x is close to x 0 , then by Taylor’s series, we have

34
 y ' ' x0 x  x0  2 y ' ' ' x0 x  x0  3
y x   y x0   y ' x0 x  x0     .... 
2! 3!

 dy   d2y 
where y ' x0     , y ' ' x0    2  , etc
 dx  x  x0  dx  x  x0

If x  x1 is close to x 0 , substitute x  x1 in (*) and get y1  y( x1 )

Again starting from x1 , express y (x ) in a power series of x  x1 and then substitute x  x2 to

get y2  y( x2 ) . In this way we can get the sequence of y values y 0 , y1 , y 2 ,...

If x  x0  0 , we get the Maclaurin’s series expansion,

x2
y x   y 0  xy' 0  y ' ' 0  ...
2!

dy
Example 1. Evaluate the solution of the differential equation  y2  1
dx

by taking four terms of its Maclaurin’s series for x  0 , x  0.2 ,

x  0.6 given y (0)  0 and compare this result with its exact

solution.

Solution

y'  y 2  1 y' 0  1

y' '  2 y y' y' ' 0  2 y0y' 0  0

y ' ' '  2 y '2 2 yy' ' y' ' ' 0  2

y 4   2 yy' ' '6 y ' y ' ' y 4  0   0

By Maclaurin’s series, we have

x2 x3
y  x   y 0  xy' 0  y ' ' 0  y ' ' ' 0  ...
2! 3!

35
 x 3 16x 5
 0 x2   ...
3! 5!

y ( 0)  0

y0.2  0.2 
0.23  2 0.25  0.2  0.008  2 0.00032  0.4213
3 15 3 15

y0.6  0.6 
0.63  2 0.65  0.6720
3 15

dy
Exact solution y 2
1
  dx  tan 1 y  x  c  y  tan x

tan 0  0, tan o.2  0.2927, tan 0.6  0.6942

Let us compare the actual value with the approximate value

Values of x 0 0.2 0.4 0.6

Actual value of y 0 0.2027 0.4228 0.6841

Approximate value of y 0 0.2027 0.4213 0.6720

Error 0 0 0.0015 0.0121

Percentage of error 0 0 0.35 1.77

This table shows that when the distance of x from x 0 increases the error also increases.

In this example, we have expanded y (x ) in the neighborhood of x  0 and used the same

result to find y (x ) when x  0.2 , and x  0.6 .

Now instead of doing this, after getting y( x1 )  y(0.2) , expand y (x ) again in the

neighborhood of x  0.2 and use this result to get y (0.6) . In doing so, we can minimize the

error.

Thus in the neighborhood of x  0.2

36
 y' ' 0.2x  0.2 2 y' ' ' 0.2x  0.2 3
yx   y0.2  y' 0.2x  0.2    ...
2! 3!

y (0.2)  0.2027

y (0.2) = ( y (0.2)) 2  1 = (0.2027)2 + 1 = 1.0411

y' ' 0.2  2 y0.2y' 0.2  20.20271.0411  0.4221

y' ' ' 0.2  2y' 0.2  2 y0.2y' ' 0.2

2

 21.0411  20.20270.4221  2.3389, …, etc.

2

Putting these values and using x  0.4 in (1), we get

y 0.4   0.2027  1.04110.4  0.2   0.4  0.22  2.3389 0.4  0.23  ...
0.4221
2 6

 0.2027  0.2(1.0411)  (0.04)(0.21105)  (0.008)(0.389871)  ...

 0.422480536 0.4225

When we compare this value with the actual one i.e. 0.4228, we see that the error is only 0.0003,

nearly 0.07%

The error has decreased from 0.35% to 0.07%

Therefore, to reduce the error, each time obtain the power series of y (x ) at x  xi 1 and use this

to get y ( xi  2 ) and so on.

This method is called the method of starting the solution.

Point Wise Methods

Consider the previous example y '  y 2  1, y (0)  0

x3
First we got y ( x)  x   ... in terms of x and then we substituted
3
37
x  x1  0.2 . Instead, without getting y (x ) as a function of x we can directly get y ( x1 ) = y (0.2)

as

0.22  y' ' ' 00.2  ...

3
y0.2  y(0)  y' (0)0.2 
y' ' (0)
2! 3!

That is we get ( x1 , y1 ), ( x2 , y2 ) directly. So, a point wise solution is a series of points

( x1 , y1 ), ( x2 , y 2 ),... which satisfy approximately a pre-assigned but not known particular solution.

Solution Using Taylor Series

AIM- To find the numerical solution of the equation

 f  x, y  given the initial condition yx0   y0 ……….. (1)

dy
dx

Now, we expand y (x ) about the point x  x0 using Taylor’s series in powers of x  x0 .

That is,

yx   yx   x  x  y' x  

 x  x0 
2
y' ' x   ...
0 0 0 0
2!

 dry 
Where y r  x0    r 
 dx  x  x
0

h 2 ' ' h3 ' ' '

y1  y x1   y0  hy01  y0  y0  ... (2)
2! 3!

, where h  x1  x0 or x1  x0  h

To find y0' , y0'' ,... we use (1) and its derivatives at x  x0 .

Even though the series in (2) is an infinite series, we can truncate it at any convenient term, if h

is small, and the accuracy can be obtained.

38
 Now, once if we get y , we can calculate y1' , y1'' ,... , etc by using

y   f ( x, y )

Again, expanding y (x ) , in a Taylor’s series about the point

x  x1 , we get

h 2 '' h3 '''
y2  y1  hy1'  y1  y1  ...
2! 3!

Proceeding in the same way, we get

h 2 '' h3 '''
yn 1  yn  hyn'  yn  yn  ...
2! 3!

Since this is an infinite series, to get an approximate value we have to truncate it at some term to

have a calculated numerical value.

Now, let us consider the terms up to and including h n and neglect terms involving h n 1 and

higher powers of h . The Taylor algorithm used this way is said to be of nth order.

Thus, the truncation error is O(h n 1 ) . If h is small enough we can neglect terms after the nth

term and get the error as

h n n 
f   where x0    x1 if x1  x0  h
n!

dy
Example 1. Solve  x  y, given y (1)  0 , and find y (1.1), y (1.2) by using
dx

Taylor’s method.

Solution

Here x 0  1 and y 0  0, h  0.1

y0  y x  1  0
dy
 y'  x  y
dx

39
 y ' '  y '1 y0'  x0  y0  1  0  1

y' ' '  y' ' y' ' 0  1  1  2

y' ' ' 0  2, etc

Thus by Taylor’s method, we have

h ' h 2 '' h3 '''

y1  y0  y0  y0  y0  ...
1! 2! 3!

y1.1  0  0.11 
0.12 2  0.13 2  0.14 2  0.15 2  ...
2! 3! 4! 5!

 0.1  0.01 0.00033 0.00000833 0.000000166 ...

 0.11033847

Now again take x0  1.1 and h  0.1

y1'  x1  y1  1.1  0.11033847  1.21033847

y1''  1  y1'  2.21033847


y1' ' '  y1  2.2103384

y1' ' '  y14   y15  ...

 y2  0.11033847 0.11.21033847 
0.1
2
2.21033847  0.1 3
2.21033847
2 3!


0.1 2.21033847  ...
4

4!

 0.11033847 0.121033847 2.21033847 (0.005)  (0.0016667) (2.21033847 )  ...

 0.2461077

Let us check for its exact solution

40
 dy
 x y
dx

dv dy
Let x  y  v   1
dx dx

dv dv
  v  1 and  dx 
dx v 1

Integrating both sides,

x  c  n v  1

 x  c  n x  y  1

 x  y  1  exc

Since y (1)  0

 2  0  e c 1

 c  1  n 2  c  n 2  1

 y   x  1  2e x 1

Hence y (1.1)  1.1  1  2e1.11  0.110341836

y (1.2)  1.2  1  2e1.21  0.242805516

Example 2. Using Taylor series method, compute y (0.1) correct to

 x 2  y 2 and y0  1 .
dy
four decimal places, given
dx

Solution

dy
 y '  x 2  y 2 here x 0  0 and y0  1, h  0.1
dx

y ' 0  02  12  1

y' '  2 x0  2 y0 y0'  2

41
y0'' '  2  2 y0'  2 y0 y0' '  8

y04  4 y0' y0''  2 y0' y0''  2 y0 y0'''  6 y0' y0''  2 y0' y0''  28

1  0.1 2  0.1 8  0.1 28  ...

2 3 4
y0.1  1 
0.1
1 2 3! 4!

 1  0.1  0.01 0.0013333333  0.000116667  1.11144999

 1.11145

Example 3. Using Taylor method, compute y (0.2) and y (0.4)

correct to four decimal places given

 1  xy and y0  0
dy
dx

Solution

y '  1  2 xy y0'  1  2 x0 y0  1

y' '  2 y  xy' y0''  20  01  0

y' ' '  2 y' y' xy' ' y0'''  22  0  4

y 4   23 y ' ' xy' ' ' y04  0


y 5  2 4 y ' ' ' xy4  y05  32, etc

By Taylor’s series, we have

h 2 ' ' h 3 ' ' ' h 4 4 

y1  y0  hy0'  y0  y0  y0  ...
2! 3! 4!

y0.2  0  0.21 
0.22 0  0.23  4  0.24 0  0.25 33  ...
2! 3! 4! 5!

 0.2  0.00533333 0.000085333

42
  0.194752003

Now again starting with x  0.2 , we have

x  0.2 , y 0  0.194752003, h  0.2

y0'  1  2 x0 y0  1  20.20.194752003  0.9220992

 
y0''  2 x0 y0'  y0  20.20.9220992  0.194752063  0.758343806

 
. y0' ' '  2 2 y0'  x0 y0''  20.2 0.7583436806  20.9220992

 3.38505933

y04  2(0.2)(3.38505933)  (0.758343686  5.90408585

 y 0.4  0.194752003 (0.2)(0.9220992) 

0.22 (0.758343686) 
2
0.2  3.38505933  0.24 5.90408585  ...
3

3! 4!

 0.359883723

Example 4. Using Taylor series method, find y at x  0.1, x  0.2

 x 2  y, y 0   1 (correct to five decimal places)

dy
and x  0.4 given
dx

Solution

Here x0  0, y 0  1, h  0.1, x1  0.1, x 2  0.2,...

y'  x 2  y y0'  1

y' '  2 x  y' y0''  1

y' ' '  2  y' ' y0'''  2  1  1

43
y 4    y ' ' ' y04  1

y 5   y 4  0 y05  1, etc

y (0.1)  1  (0.1)(1) 
0.1
2

 1  0.1 1 
0.12 1  0.1 3 1  0.1 4  1  0.1 5 1  ...
2 3! 4! 5!

= 1+ (-0.1) + 0.005 + 0.0001666 + (-0.0000416) +…

= 0.905125

Now again using x1  0.1 and y1  0.905125, we have

y1'  x12  y1  0.01  0905125  0.895125

y1''  2 x1  y1'  0.2  (0.895125)  1.095125

y1'''  2  y1''  2  0905125  0.904875

y14    y1'''  0.904875, etc.

 y0.2  0.905125 0.1 0.895125 

0.1
2
1.095125 
2

0.13 0.904875  0.14  0.904875  ...

= 0.821235167
3! 4!

Similarly, y (0.3)  0.7492 and y (0.4)  0.6897

Taylor series method for simultaneous first order differential equations

The equation of the type

 f  x, y , z  and  g  x, y , z 
dy dz
dx dx

with initial conditions

y ( x0 )  y 0 , z x0   z 0

44
 can be solved by Taylor series method as given below:

 y  x with y0  1, z(0)  1 by taking

dy dz
Example 5. Solve  z  x,
dx dx

h  0.1 , to get y (0.1) and z (0.1) .

Solution

y'  z  x and z'  x  y

x 0  0 and y0  1 x0  0 , z 0  1 and h  0.1

y1  y(0.1)  ? z1  z0.1  ?

y'  z  x z' '  1  y'

y ' '  z '1 z' ' '  y' ' etc

y' ' '  z' '

Using Taylor series, for y1 and z1 , we have

h 2 '' h 3
y1  y 0.1  y 0  hy  '
0 y0  y 0  ... …… (1)
2! 3!

h 2 '' h3 '''
and z1  z (0.1)  y0  hy0'  z0  z0  ... ..... (2)
2! 3!

y0  1 z0  1

y01  z0  x0  1 z0'  x0  y0  0  1  1

y0''  z0'  1  1  0 z0' '  1  y0'  1  1  2

y0'''  z0''  2 z0'''  y0''  0

z04  y0'''  2

Substituting these in (1) and (2), we get

45
 y1  1  0.1 
0.01
0  0.0012  0.00010  ...
2 6 24

= 1+0.1+0.000333+…

= 1.1003

z1  1  0.1 
0.01
2  0.0010  0.00012  ...
2 6 24

=1+ 0.1 + 0.01+0.0000083+… = 1.1100

Example 6. Find y0.1, y0.2, z0.2 given

 x  y 2 and y0  2, z0  1

dy dz
 x  z,
dx dx

Solution

Here x 0  0, y 0  2 and z 0  1

y'  x  z z'  x  y 2

y' '  1  z' z ' '  1  2 yy'

y' ' '  z' ' z' ' '  2 y'  2 yy' '
2

y 4   z ' ' ' etc z 4   4 y ' y ' '  2 y ' y ' '2 yy' ' '  6 y ' y ' '2 yy' ' etc

y0'  y ' 0  x0  z0  1 z0'  0  22  4

y0'  1   4  3 z0''  1  221  3

y0''  3 z0'''  21  22 3  10

2

y04  10 z04  61 3  22 3  30

y05  30 z05  74

y(0.1)  2  (0.1)(1) 
0.12  3  0.13  3  0.14 10  0.15 30  ...
2 6 24 120

46
 = 2.084544167

z (0.1)  1  (0.1)(4) 
0.1
2
 3  0.1
3
10  0.1
4
30  0.1
5
 74  ...
2 6 24 120

= 2.084544167+0.267132966-0.016226621-0.0003105449027

- (0.000003091501458) +…

y(0.2)  0.5867855 (0.1)(4.245324384) 

0.1
2
 1.863269416  0.1
3
 0.74196035  0.1
4
 ...
2 6 24
=0.152817221

Taylor series method for second order differential equation

Any differential equation of the second order or higher can be solved by reducing it to

a lower order differential equation. A second order differential equation can be reduced to a first

order differential equation by transforming y  z and then the given equation can be solved so

easily.

d2y  dy 
Suppose  f  x, y,  …… (1)
 dx 
2
dx

y' '  f x, y, y' is the given differential equation together with the given initial

conditions

y ( x0 )  y0 and y( x0 )  y0 where y 0 , y0 are known values at x 0 .

Setting y  p , we get y  p and (1) becomes

p '  f ( x, y , p )

with initial condition y ( x0 )  y0 and p( x0 )  p0 = y0

Using Taylor series method, we get

h 2 '' h3 '''
p1  p0  hp0'  p0  y0  ... where p1  px  x1  & x1  x0  h …… (*)
2 3!

47
 h 2 '' h3 '''
y1  y0  hy  '
0 y0  y0  ... becomes
2 3!

h2  h3 
y1  y 0  hp0  p0  p0  ... ...(**)
2! 3!

Since p'  f ( x, y, p ) and taking the derivative of p  again and again with respect to x, we get

p, p , etc.

Hence p0' , p0'' , p0''' ,... can be solved using (**) and (*), so that we can get

y1 & p1 .

Once knowing y1 and p1 we can get p1' , p1'' , p1''' ,... at x1 , y1 
Again using

h 2 '' h3 '''
p2  p1  hp1'  p1  p1  ... we get some value for p2 and using
2! 3!

h 2 '' h3 '''
y2  y1  hy1'  y1  y1  ... we get still some value for y2 .
2! 3!

Example 7. Evaluate the values of y (0.1) and y (0.2) given

y' ' ' x y'  y 2  0, y0  1, y' 0  0 by using Taylor series method.
2

Solution

First put y   z and hence the equation reduces to z ' xz 2  y 2  0

 z '  xz2  y 2

Using the initial condition y0  1, z0  y0  0

'

Now z '  xz  y can be solved given that z0  z 0   0, and x0  0

2 2

h 2 ''
Here, z1  z0  hz  z0  ... ...1
'
0
2

48
z '  xz2  y 2 y'  z

z ''  z 2  2 xzz'2 yy' y ' '  z ' , y ' ' '  z ' ' ,....


z' ' '  2zz'2 xzz' ' zz' xz'  2 yy' ' y'
2
  2

 z0'  x0 z02  y02  1

z0''  0 and z0  2
'''

Substituting, these values in (1), we get

z1  0  0.1 1 
0.01 0  0.001 2  ...
2 6

= -0.0997

By Taylor series for y1 , we get

h 2 ''
y1  y (0.1)  y0  hy0'  y0  ...
2

 1  (0.1)( z ) 
0.1 ' 0.1 ''
2
z  
3
z   ...
0 0 0
2 6

 1  0.10  
0.01 0.001
(1)  (0)  ... = 0.995
2 6

Similarly,

2
y2  y  x2   y1  hy1' 
h
y1''  ...
2

 0.995  0.1 z1   ...2 

0.01 ' 0.001 ''
z1  z1  ...
2 6

z1'  x1z12  y12  0.1 0.0997  0.995  0.9890

2 2

z1''  0.1687
49
Using (2),

y2  0.995  0.1 0.0997   0.9890  0.001 0.1687  ....

0.01
2 6

= 0.9801

Example 8. Solve y' '  y  xy' given y0  0 and calculate y (0.1) .

Solution

Here x0  0, y0  1, y0  0 and y ' '  y  xy'

'

Differentiating with respect to x,

y ' ' '  y ' y ' xy' '  2 y ' xy' ' y0''  y0  x0 y0'  1

y 4   2 y ' ' y ' ' xy' ' '  3 y ' ' xy' ' ' y0'''  2 y0'  x0 y0''  0

y 5  4 y ' ' ' xy4  y04  3

y 6  5 y 4  xy5 y05  0 y06  15,...

x 2 '' x 3 '''
Here y x   y0  xy0'  y0  y0  ...
2 3!

x2 4
 1 0  1  0  x 3  ...
2! 4!

x2 x4 x6
 1    ...
2 8 48

 1
0.12  0.14  0.16  ....  1.00501252
2 8 48

Exercise

Using Taylor series method, find the values required in each problem.

dy
1. Find y (0.1) given  x  y , y ( 0)  1
dx

50
 2. Find y (0.1) given y   x 2 y  1, y (0)  1

dy 1
3. Obtain y ( 4.2) and y ( 4.4) given  2 , y(4)  4 taking h  0.2 .
dx x  y

x 3  xy 2
4. Find y (0.1), y (0.2), y (0.3) given y   , y (0)  1 .
ex

5. Find y (0.1), y (0.2), z (0.1), z (0.2) given

dy dz
 x  z,  x  y 2 and y (0)  2, z (0)  1
dx dx

6. Evaluate x(0.1), y (0.1), x(0.2), y (0.2) given

dx dy
 ty  1,  tx given x  0, y  1 at t  0
dt dt

x dy
7. Solve for x and y  x  y  t,  2 x  t given x  0, y  1 at t  1 .
dt dt

8. Find y at x  1.1, x  1.2, x  1.3 given y   y 2 y   x 3 , y (1)  1, y (1)  1

9. Express y as a power series given, y   (0.1)( x 2  y 2 ), y (0)  1

Picard’s Method of Successive Approximations

 f  x, y  subject to y ( x0 )  y0
dy
AIM: To solve
dx

 f  x, y   dy  f  x, y  dx
dy
Now
dx

Integrating, y   f x, y  dx  c
x

Setting x  x0 , we have

y0   f x, y dx  c
x0

 y  y 0   f x, y dx
x

x0

51
  y  y0   f x, y dx this type of equation is called an integral equation.
x

x 0

As the integration is not possible as it is, we will solve it numerically by successive

approximation. Now substitute the initial values of y namely y 0 in the integrand f ( x, y ) in

place of y and then integrate it to get an approximate value of y .

i.e. y 1  y0  x f x, y0 dx

x

After getting the first approximation y (1) for y , use this value of y (1) in place of y in

f ( x, y ) and then again integrate to get the second approximation of y namely y ( 2) .

Thus y  y0  x f x, y dx

2  1 x

Repeating this procedure again and again, we eventually get


y n   y0   f x, y n 1 dx
x

x0

This equation is called Picard’s iteration formula. This formula gives the general iterative

formula for y.

The sequence y (1) , y ( 2) ,…, y (n ) should converge to y (x ) ; otherwise the process is not valid.

f
The condition for the convergence of the sequence is both f ( x, y ) and are continuous.
y

f
i.e. f x, y   k1 a n d  k2 in a region containing the point ( x0 , y0 ) where k1 , k2 are
y

constants.

Example 1. Solve y  y  x2 , y (0)  1 , by Picard’s method up to the

third approximation. Hence find the value of y (0.1), y (0.2) .

52
 Solution

y'  y  x 2

 y  y0  
x

x0
y  x dx....
2

Hence, x0 = 0, y0 = 1 and f ( x, y )  y  x 2

 y  1    y  x 2 dx .....1
x

y 1 x

 1   1  x dx  1  x 
0
x3
3
2

Using y(1) again in (1), we get

y 2   1    y 1  x 2 d x
x

 x x3 
 1   1  x   x 2 dx

0
 3 

x 2 x3 x 4
 1 x   
2 3 12

Using again this result, we get

y 3   1    y 2   x 2 d x
x

x x 2 x3 x 4 
 1   1  x     x 2 dx
0
 2 3 12 

x 2 x3 x 4 x5
 1 x    
2 6 12 60

Now put x  0.1

y0.1  1  0.1 
0.12  0.13  0.14  0.15
2 6 12 60

53
  1.104824833

y0.2  1  0.2 
0.22  0.23  0.24  0.25
2 6 12 60

= 1.218528

Note- In getting the value y (0.2) we could have started with x0 = 0.1 and

y 0 = 1.104824833 to get a closer value of y (0.2)

x
i.e. y  1.104824833 0.1 y0  x dx
2
 
x
1  x3 
y  1.104824833  y0 x  
 3 0.1

x3 0.13
 1.104824833 1.104824833x   1.104824833(0.1)  
3 3

x3
 0.99467574 1.104824833x 
3

y 2   1.104824833 
x

0.1
y    x dx
1 2

x x3 
 1.104824833   0.99467574 1.104824833x   x 2  dx
0
 3 

 1.104824833 0.99467574( x  0.1)  1.104824833x2  0.1  2



12
 4 1
3
 
x  0.1  x 3  0.1
1 4 3

y 2  0.2  1.2184066

dy
Example 2. Solve  x  y, given y (0)  1 . Obtain the values of
dx

y (0.1) and y (0.2) using Picard’s method and check your

54
 result with the exact solution.

Solution

Here f ( x, y )  x  y, x0  0 and y 0  1

 y  y0   f x, y dx
x

x0

 1   f  x, y dx
x

y 1  1   f  x, y dx
x

 1    x,1dx
x

x2
= 1   x  1 dx  1  x 
x

0 2

x x2 
y 2   1   1  x   x  dx
0
 2 

x3
 1 x  x  2

x x3 
y 3  1    x  1  x  x 2   dx
0
 6

x3 x 4
 1  x  x2  
3 24

x3 x 4
 y x   1  x  x 2    ...
3 24

Now put x  0.1

y0.1  1  0.1  0.1 

2 0.13  0.14  ...
3 24

55
 0.001 0.0001
 1  0.1  0.01    ...
3 24

= 1.1103374

To find y (0.2) we use x0 = 0.1 and y 0 = 1.1103374

y  y0   f  x, y dx
x

0. 1

x  y dx
x
 1.1103374 
0.1

y 1  1.1103374  x  1.1103374dx

x

0.1

 1.1103374 1.1103374( x  0.1)  ( x 2  (0.1) 2 )

= 0.98930366 1.1103374x  x 2


y 2   1.1103374  x  0.98930366 1.1103374x  x 2 dx
x

0.1


 1.1103374 0.98930366( x  0.1)  2.1103374( x 2  (0.1) 2 )  x  0.1

1 3
3
3

1
 0.989970326 0.98930366x  2.1103374x 2  x 3
3

x  1 
y 3  1.1103374   0.989970326 1.98930366x  2.1103374x 2  x3 dx
0.1
 3 

 1.1103374 0.989970326( x  0.1)  1.98930366( x 2  (0.1) 2 )


 2.1103374 x 3  0.1 
3
 12
x  0.1
1 4 4

 y 3 0.2  1.1103374 0.989970326(0.1)  1.98930366(0.03)

 2.1103374(0.007)  0.0015
1
12

= 1.283910904

56
 dy
Solving  x  y analytically, we get y  2e x  x  1
dx

y (0.1) = 1.110341836 and

y (0.2) = 1.242805516

 x 2  y 2 , y 0   1 by Picard’s method
dy
Example 3. Solve
dx

Solution

Here x0  0, y0  1

 y  y0   f x, y  dx = 1   x 2  y 2 dx
x

x0
x

0
 
Now


y 1  1   x 2  1 dx  1  x 
x

0
 x3
3

  x3  
2

y 2   x
 1   x  1  x   dx
2
0  3 
  

2 x3 1 4 1 7 2 5
 1  x  x2   x  x  x
3 6 63 15

2 3 1 4 2 5 1 7
 1  x  x2  x  x  x  x  ...
3 6 15 63

Example 4. Solve y   y  e x , y (0)  0 , by Picard’s method

Solution

By Picard’s method
x

x0
x

0

y  y0   f ( x, y )dx   e x  y dx

 
y 1   e x  0 dx  e x  1
0
x

57
  
y 2    e x  e x  1 dx  x
0
x


0
x

y 3   e x  x dx  e x   x2
2
1

x   x2 
y 4     e x   e x   1  dx
0
  2 

x3
 x
6

5  x  x  x3 
y   e    x   dx
 6 
0
  

x2 x4
e  
x
1
2 24

x2 x4
 y x   e x   1
2 24

Exercise

1. Using Picard’s iterative formula,

dy
a) Solve  x  y 2  1 , given y (0)  0 .
dx

yx
b) Obtain y (0.1) given y   and y (0)  1
yx

c) Solve y   1  2 xy given y (0)  0

2. Find the values of y for x  0, x  0.1 and x  0.5 given y   x 2  y 2 and which passes

through (0,1) .

x2
3. Given y   and y (0)  0 , find y (0.25), y (0.5)
1 y2

58
Euler’s Method

In solving a first order differential equation by numerical methods, we have two types of
solutions:
i) Values of y at specified values of x
ii) A series solution of y in terms of x , which will yield the value of y at a particular
value of x by direct substitution in the series solution.
Taylor and Picard’s method studied so far belong to the first category in finding the numerical
solution of differential equations; the methods due to Euler, Runge-kutta, Adam-Bash Forth and
Milne come under the second category.
The methods of second category are called step-by-step methods because the values of y are
calculated by short steps ahead of equal interval h of the independent variable x .
Euler’s method

 f  x, y  with initial conditions y ( x0 )  y0

dy
Suppose we want to solve
dx
Let us take the points

x  x0 , x1 , x2 ,... where xi  xi 1  h
y
i.e. xi  x0  ih , i  0,1,2,... P2

P1
Q2
P0
(x1,y1) Q1

M2
M0 M1
x0 x1 x2 x
0

Let the actual solution of the differential equation be denoted by the graph Pi  Pi ( xi , yi ) lies on

the curve. We require the value of y the curve at x  xi .

The equation of the tangent line at ( x0 , y0 ) to the curve is

59
 y  y0  y' x0 , y0  x  x0 

= f x0 , y0 x  x0 

 y  y0  f x0 , y0 x  x0 

In the interval ( x0 , x1 ) , the curve is approximated by the tangent.

 The value of y on the curve is approximately equal to the value of y on the tangent

at ( x0 , y0 ) corresponding to x  xi .

 y1  y0  f x0 , y0 x1  x0 

i.e. y1  y0  hy0'

Again, we approximate the curve by the line through ( x1, y1 ) and whose slope is f ( x1, y1 ) , we

get

y2  y1  f x1, y1 h

 y1  hy1'

Thus, yn 1  yn  hf xn , yn , n  0,1,2,...

This formula is called Euler’s algorithm.

In other words, y ( x  h)  y ( x)  hf ( x, y )

In this method, the actual curve is approximated by a sequence of short straight lines.

As the intervals increase the straight line deviates much from the actual curve. Hence the

accuracy cannot be obtained as the number of intervals increase.

Referring to the above graph

Q1 P1  error at x  x1

x1  x0 2 y' ' x , y   h2 y' ' x , y 

= 1 1 1 1 it is of order h2.
2! 2

60
 Example 1. Given y   y and y (0)  1 , determine the values of y at

x  0.01, x  0.02 and x  0.04 by Euler method.

Solution

y    y and y (0)  1 ; f ( x, y )   y

Here x0  0, y0  1, x1  0.01, x2  0.02, x3  0.03, x4  0.04

We have to find y1 , y 2 , y3 , y 4 . Take h  0.01

By Euler algorithm,

yn1  yn  hyn'  yn  hf xn , yn 

y1  y0  hf x0 , y0   1  0.01 1  1  0.01  0.99

y2  y1  hy1'  0.99  (0.01)( y1 )

 0.99  0.01 0.99  0.9801

y3  y2  hf x2 , y2   0.9801 0.01 0.9801  0.9606

y4  y3  hf x3 , y3   0.9703  0.01 0.9703  0.9606

Let us compare the results

x 0 0.01 0.02 0.03 0.04

y 1 0.9900 0.9801 0.9703 0.9606

y  ex 1 0.9900 0.9802 0.9704 0.9608

Example 2. Using Euler’s method, solve numerically the equation,

y   x  y, y (0)  1 for x  0, x  0.2 and x  1.0 .

Check your answer with the exact solution

61
 Solution

Here h  0.2 , f ( x, y )  x  y, x0  0, y 0  1

x1 = 0.2 , x 2 = 0.4 , x3 = 0.6 , x 4 = 0.8 , x5 = 1.0

By Euler algorithm,

y1  y0  hf x0 , y0   y0  hx0  y0 

= 1  0.2(0  1)  1.2

y2  y1  h( x1  y1 )  1.2  (0.2)(0.2  1.2)

= 1.48

y3  y2  hx2  y2   1.48  (0.2)(0.4  1.48)

= 1.856

y4  y3  h( x3  y3 )  1.856  (0.2)(0.6  1.856)

= 2.3472

y5  2.3472  0.2(0.8  2.3472) = 2.94664

Exact solution is y  2e x  x  1

x 0 0.2 0.4 0.6 0.8 1.0

Euler y 1 1.2 1.48 1.856 2.3472 2.94664
Exact y 1 1.2428 1.5836 2.0442 2.6511 3.4366

As you can see from the table, the values of y deviate from the exact values as x increases. To

avoid this discrepancy we need to improve Euler’s method.

Improved Euler Method

Let the tangent at ( x0 , y0 ) to the curve be P0A. In the interval ( x0 , x1 ) by the previous

Euler’s method, we approximate the curve by the tangent P0A.

62
  y11  y0  hf x0 , y0  where y11  M1Q1
y

D
P1 B
(x0,y0) C
A
Q1 (x1y )
P0

x0 x1
x
M0 M1

Let Q1C be the line at Q1 whose slope is f x1 , y11  

Now take the average of the slopes at P0 and Q1

i.e.
1
2
 
f  x0 , y0   f x1 , y11 

Draw a line P0D through P0 ( x0 , y0 ) with this as the slope.

That is,

y  y0 
1
2
 
f  x0 , y0  f x1 , y11  x  x 
0

and this line intersects x  x1 at

y1  y0 
1
2
 
h f  x0 , y0   f x1 , y11 

h f  x0 , y0   f  x1 , y0  hf  x0 , y0 
1
 y0 
2

In general,

h f  xn , yn   f  xn  h, yn  hf  xn , yn 
1
y n 1  y n 
2

This is what is known as improved Euler’s method.

63
 Notice that the difference between Euler’s method and the improved Euler’s method

 
is that in the improved one we take the average of the slopes at ( x0 , y0 ) and x1 , y11 instead of

the slope at ( x0 , y0 ) in the former method.

Modified Euler’s Method

In the improved Euler method, we arranged the slopes, whereas in modified Euler method, we

will average the points.

Let P0 ( x0 , y0 ) be the point on the solution curve

Let P0A be the tangent at ( x0 , y0 ) to the curve. Now let this tangent meet the ordinate at

h at N 1 and y coordinate of N1  y0  hf  x0 , y0  . Calculate the slope at N1. i.e.

1 1
x  x0 
2 2

 
f  x0  h, y0  hf x0 , y0  .
1 1
 2 2 

k1(x1,y1)
B
(x0,y0)
P0 A
N1

M0 x
0 x0 x0+ h x1=x0+h


Let this line meet x  x1 at k1 x1 , y11 . 
This y11 is taken as the approximate value of y at x  x1

  
 y11  y0  h f  x0  h, y0  hf x0 , y0 
1 1
  2 2 

64
In general,

  
yn 1  yn  h f  xn  h, yn  hf xn , yn 
1 1
  2 2 

  
or yx  h  yx   h f  x  h, y  hf x, y 
1 1
  2 2 

This is called modified Euler’s formula.

Note: The Euler predictor is yn 1  yn  hyn' and the corrector is yn 1  yn 

h '
2
 
yn  yn' 1 in the

improved Euler method.

When you read some literature there is some confusion among the authors. Some take the

improved Euler method as the modified Euler method and the modified Euler method is not

mentioned at all.

Example 3 Solve numerically y '  y  e , y 0  0 for x  0.2, x  0.4 by improved

x

Euler’s method.

Solution

y'  y  e x , y0  0 x0  0, y0  0, x1  0.2, x2  0.4 and h = 0.2.

By improved Euler method,

h f  x0 , y0   f  x1 , y0  hf  x0 , y0 
1
y1  y0 
2

 y1  0 
0.2
2
  
y0  e x 0  y0  h y0  e x 0  e x 0  h 

 0.1 0  1  0  0.20  1  e0.2 
y0.2  y1  0.24214

65
 y 2  y1  h f  x1 , y1   f  x1  h, y1  hf  x1 , y1 
1
2

Here f  x1 , y1   y1  e  0.24214  e0.2  1.46354

x1

y1  hf x1, y1   0.24214 0.21.46354  0.53485

f x1  h, y1  hf x1 , y1   0.53485  e0.4  2.02667

y2  0.24214 0.11.46354 2.02667

y(0.4)  y2  0.59116

Example 4. Compute y at x  0.25 by Modified Euler’s method

y  2 xy , y (0)  1 .

Solution

f ( x, y )  2 xy; x0  0, y 0  1

Take h  0.25, x1  0.25

By modified Euler method,

  
yn 1  yn  h f  xn  , yn  hf xn , yn 
h 1
  2 2 

  
 y1  y0  h f  x0  , y0  hf x0 , y0 
h 1
  2 2 

f ( x0 , y0 )  f (0,1)  0

y1  1  (0.25) f (0.125,1)

= 1  (0.25)(2)(0.125)(1)

66
 y (0.25)  1.0625

dy
y '  2 xy   2 xdx
y

Using y (0)  1

y ( 0)  c  1

 y  ex
2

 y0.25  e 0.25  1.064494459

2

The error is only 0.001999445891.

dy
Example 5. Solve the equation  1  y, given y (0)  0 using modified
dx

and Euler’s method and tabulate the solutions at x  0.1, x  0.2 ,

and x  0.3 . Compare these results with the exact solutions.

Solution

Here x0  0, y0  0, x1  0.1, x2  0.2, x3  0.3, h  0.1

y  1  y

 f x, y   1  y,  f x0 , y0   1  y0  1

i) Modified Euler method

 
y1  y0 h f  x0  , y0  h f x0 , y0 
h 1
 2 2 

h  0.1  0.05 and y0  hf  x0 , y0   0.05

1 1 1
x0 
2 2 2

 y1  0  0.1 f 0.05,0.05  0.11  0..05  0.095

67
Now again f ( x1 , y1 )  1  y1 = 1- 0.95 = 0.905

 
y2  y1  h f  x1  h, y1  h f x1 , y1 
1 1
 2 2 

 0.095  0.1 f 0.150.14025

 0.095  0.11  0.14025

= 0.180975

 
y3  y2  h f  x2  h, y2  h f x2 , y2 
1 1
 2 2 

 
 0.180975 0.1 f  0.2  f 0.2,0.180975
0.1 0.1
,0.180975
 2 2 

= 0.180975 (0.1) f(0.25,0.22192625)

= 0.258782375

ii) Improved Euler method

h  f  xn , yn   f  xn  h, yn  h f  xn , yn 
1
yn 1  yn 
2

h  f  x0 , y0   f  x0 h, y0  h f  x0 , y0 
1
 y1  y0 
2

f x0 , y0   1  y  1

f x0  h, y0  h f x0 , y0   f 0.1,0  0.11  f 0.1,0.1

= 1 - 0.1  0.9

 y1  y 0.1  0  1  0.9  0.095

0.1
2

h  f  x1 , y1   f  x1  h, y1  h f x1 , y1 
1
y2  y1 
2

68
 f x1 , y1   1  y1  1  0.095  0.905

f x1  h, y1  hf x1 , y1   f 0.2,0,905  0.11  0.095

= f (0.2,0.1855) = 0.8145

 y2  0.095 
0.1
0.8145  0.905  0.180975
2

y3  y2 
1
h  f x2 , y2   f x2  h, y2  h f x2 , y2 
2

 0.180975
0.1
1  0.180975  f 0.3,0.180975 0.11  0.180975
2

 0.180975 0.05(0.819025 1  0.2628775)

 0.258782375

Exact solution

dy dy
 1 y   dx
dx 1 y

 n1  y   x  c1  1  y  ce  x

At x  0, 1  0  ce 0  c  1

 y  1  e x

 y0.1  0.095162581

y (0.2)  0.181269246

y (0.3)  0.259181779

x Modified improved exact

0.1 0.095 0.095 0.095162581

0.2 0.180975 0.180975 0.181269246

69
0.3 0.258782375 0.258782375 0.259181779

Example 6. Find, correct to four decimal places, the value of y (0.1)

Given y  x 2  y , y (0)  1 by using improved Euler method.

Solution

Here, h  0.1 , x 0  0 , x1  0.1 and y 0  1

By the improved Euler method,

h f  xn , yn   f  xn  h, yn  hf  xn , yn 
1
y n 1  y n 
2

h f  x0 , y0   f  x0  h, y0  hf  x0 , y0 
1
y1  y0 
2

 1
0.1 2
2
 
x0  y0  f 0.1, y0  0.1 x02  y0  

 
 1  0.05  1  0.1  1  0.1 1  0.9055
2

Example 7. Using improved Euler method find y at x  0.1 and at

 y  , y0  1
dy 2x
x  0.2 given
dx y

Solution

By improved Euler theorem,

h  f  xn , yn   f  xn  h, yn  h f  xn , yn 
1
y n 1  y n 
2

Here x0  0, y 0  1, h  0.1

h f  x0 , y0   f  x0  h, y0  hf  x0 , y0 
1
 y1  y0 
2

0.1  20 0.112  0

 1 1 1 
2  1 1.1 

70
 f x0  h, y 0  hf xn , y n   f 0.1,1  0.11  f (0.1,1.1)

20.1
 1.1   0.918181818...
1.1

= 0.918182

 y1  1 
0.1
1  0.918182  1.095909091
2

y2  y1 
h
 f  x1 , y1   hf x1  h, y1  hf  x1 , y1 
2

20.1
f  x1 , y1   1.095909091
1.09590909

 0.913412201

f x1  h, y1  hf x1, y1   f 0.2,1.095909091 0.10.913412201

 f 0.2,1.187250311

20.2 
 1.187250311
1.187250311

 0.850337365

Thus y2  1.095909091
0.1
0.913412201 0.850337365
2

 1.184096569

Example 8. Using modified Euler method, find y (0.1), y (0.2)

 x 2  y 2 , y 0   1
dy
given
dx

Solution

Here, x0  0, y0  1, h  0.1, x1  0.1 , f ( x, y )  x  y

2 2

71
By modified Euler method,

 
y1  y0  hf  x0  h, y0  f x0 , y0 
1 h
 2 2 


 1  0.1 0.052  1  0.051
2

y (0.1) = 1.1105

 
y2  y1  hf  x1  h, y1  f x1 , y1 
1 h
 2 2 


 1.1105  0.1 0.152  1.1726605132 
 1.250263268

Exercise

1. Use Euler’s method to find

a) y (0.4) given y   xy, y (0)  1

b) y (1.5) taking h  0.5 given y   y  1, y (0)  1.1

dy 2x
2. Compute y (0.3) taking h  0.1 given  y  , y (0)  1 using improved Euler
dx y

method.

3. Find y (0.6), y (0.8) and y (1) given y   x  y, y (0)  0 taking h  0.2 by improved

Euler method.

yx
4. Use improved Euler method to find y (0.1) given y   , y (0)  1 .
yx

dy
5. Using improved Euler method find y (0.2), y (0.4) given  x y , y ( 0)  1 .
dx

6. Use improved Euler method to calculate y (0.5) ,taking h  0.1 and

y   y  sin x, y (0)  2

72
 7. Use modified Euler method and obtain y (0.2) given y   log( x  y ), y (0)  1, h  0.2 .

dy y
8. Use improved and modified Euler method, to get y (1.6) if  y 2  , if y (1)  1 .
dx x

9. Solve y   3x 2  y given y (0)  4 if h  0.25 to obtain y (0.25), y (0.5) .

y 5 2 3 1
10. Given y    x y , y (1)  find y ( 2) if h  0.125 .
x 2 2

11. Find y (0.2) by improved Euler method, given y    xy 2 , y (0)  2

if h  0.1 .

Runge-Kutta Method

i) Second order Runge-kutta method

 f  x, y  given y  x0   y0 .......... 1

dy
Suppose
dx

By Taylor series, we have

y  x  h   y  x   hy'  x  
h2
2!
 
y ' ' x   0 h3 .......... ( 2)

Differentiating (1) with respect to x,

f f dy
y' '    f x  y ' f y  f x  ff y ........ 3
x y dx

Using the values of y and y  derived from (1) and (3), in (2) we get

y  x  h   y  x   hf  h  f x  ff y   O h 3
1 2
2
 

 y  hf  h  f x  ff y   O h 3
1 2
2
  ........ 4 

Let 1 y  hf x, y   f x, y x  k1

2 y  hf x  mh, y  mk1   k2

73
and let y  ak1  bk2 where a, b and m are constants to be determined to get the better

accuracy of  y.

Now expand k 2 and y in powers of h, by Taylor series for two variables, we have

k2  h f x  mh, y  mk1 

    
2

  mh  mk1  f 
      x y  
 h  f  x, y    mh  mk1  f   ...
 x y  2!
 
 

    
2

 
 mh  mk  f 
x y 
1
   
 h  f  mhf x  mhff y   ...
2!
 

 


= h f  m h f x  mh ff y  ... higher powers of h .

2 2

Substituting k1 , k2 in y , we get

  
y  ahf  b hf  mh 2  f x  ff y   O h3

 
 (a  b)hf  bmh2  f x  ff y   O h3  ak1  bk2 .....5

Equating (4) and (5), we have

1
a  b  1, mb 
2

1
 a  1  b and m 
2b

 y  1  bk1  b k2

74
Where k1  hf ( x, y)

 h hf 
k2  h f  x  ,y 
 2b 2b 

Now y  yx  h  yx 

 hf 
 yx  h  yx   1  bhf  bhf  x  , y  
h
 2b 2b 


i.e. yn 1  yn  1  b h f xn , yn   b h f  xn 
h
, yn 
h 
 
f xn , yn   O h3
 2b 2b 

1
From this general second order Runge-kutta formula, setting a  0 , b  1 , m  , we get the
2

second order Runge-kutta algorithm.

k1  hf ( x, y )
 1 1 
k2  hf  x  h, y  k1 
 2 2 

and y  k2 where h  x

ii) Third order Runge-kutta Method

For n = 3, a similar derivation to the one as the second-order method can be performed. Since

the derivation is tedious we state simply the formula.

y n 1  y n 
1
k1  4k2  k3 
6

, where k1  h f ( x, y)

 1 1 
k2  hf  x  h, y  k1h 
 2 2 

k3  hf ( x  h, y  k  2k2 )

iii) Fourth Order Runge-kutta method

75
The most popular and commonly used form is the classical fourth-order Runge-kutta method.

y n 1  y n 
1
k1  2k2  2k3  k4 
6

, where k1  hf x, y 

 1 1 
k2  hf  x  h, y  k1h 
 2 2 

 h 1 
k3  hf  x  , y  k2 h 
 2 2 

k4  hf x  h, y  k3h 

Note 1. The second order Runge-kutta method,

 h 1 
y0  k2  hf  x0  , y0  k1 
 2 2 

 
 hf  x0  , y0  hf x0 , y0 
h 1
 2 2 

 
 y1  y0  hf  x0  , y0  hf x0 , y0 
h 1
 2 2 

this is exactly the Modified Euler method.

Thus, the second order Runge-kutta method is simply the modified Euler method.

2. If f ( x, y )  f ( x) , a function of x alone, then the fourth order Runge-kutta method

reduces to

k1  hf x0 

1   h 
y  h f x0   4 f  x0    f x0  h
6   2 

76
 h
  h 
 2  f x0   4 f  x0    f x0  h 
3   2 

= the area of y  f (x) between x  x0 and x  x 0  h with two equal intervals of

h
length by Simpson’s one-third rule. i.e. y reduces to the area by Simpson’s one-third rule.
2

Example 1. Apply the fourth order Runge-kutta method to find

y (0.2) given that y   x  y, y (0)  1

Solution

Since h is not mentioned, we can take h  0.1

 f x, y   x  y, x0  0, y0  1, x1  0.1, x2  0.2

By fourth-order Runge-kutta method, for the first interval

k1  hf x0 , y0   0.10  1  0.1

 1 
k2  hf  x0  h, y0  k1   0.10.05  1  0.05  0.11
1
 2 2 

 
k3  hf  x0  h, y0  hk2   0.10.05  1  0.055  0.1105
1 1
 2 2 

k4  hf x0  h, y0  k3   0.10.1  1  0.1105  0.12105

 y 0.1  y1  y0  k1  2k2  2k3  k4 

1
6

 1
1
0.1  20.11  20.1105  0.12105
6

= 1.110341667

77
Now starting from ( x1 , y1 ) we get ( x2 , y2 ) . Again apply Runge-kutta algorithm replacing

( x0 , y0 ) by ( x1 , y1 ) .

k1  hf x1, y1   0.10.1  1.110341667  0.121034166

 k 
k2  hf  x1  , y1  2   0.10.15  1.176384605
h
 2 2

 0.13263846

 y 0.2   1.110341667
1
0.121034166 20.132085875 0.13263846  0.144298012
6

= 1.242805142

Remember that the exact solution is y  2e  x  1

x

y (0.2)  1.242805516

The difference between the exact solution and the fourth order Runge-kutta method is

0.00000037432 .

As compared to other methods, this method is the best one.

Example 2. Find the values of y at x  0.1, 0.2 using Runge-kutta

method of i) second order ii) third order and

iii) fourth order for the given that y    y and y (0)  1 .

Solution

f ( x, y )   y, x0  0, y 0  1, x1  0.1, x 2  0.2, h  0.1

i) Second order

k1  hf x0 , y0   0.1 1  0.1

78
  1 
k2  hf  x0  h, y0  k1   0.1 1  0.05  0.095
1
 2 2 

 y1  y0  k2

 1  0.095  0.905

Again let ( x1 , y1 )  (0.1, 0.905)

k1 = (0.1) (-0.905) - 0.0905

  
k2  0.1  0.905   0.0905  0.085975
1
  2 

y (0.2)  0.905  (-0.085975)  0.819025

ii) Third order

k1  hf ( x0 , y 0 )  0.1

k 2 = - 0.095

k 3  hf ( x0  h, y 0  2k 2  k1 )

 (0.1) (-(1 2(-0.095) (-0.1))  - 0.091

y 0.1  1  k1  4k2  k3 

1
6

 1
1
 0.1  4(0.095)  (0.091
6

 0.904833333

Now take ( x1 , y1 )  (0.1,0.904833333) and repeat the process

k1  hf x1 , y1   0.1,0.904833333  0.0904833333

79
  h 1 
k2  (0.1) x1  , y1  k1 
 2 2 

  
 0.1   0.904833333  0.0904833333 
1
  2 

 - 0.085959166

k 3  0.1 0.904833333  2 0.085959166   0.0904833333

 - 0.082339833

 y 2  y1 
1
k1  4k 2  k 3 
6

 0.904833333
1
 0.0904833333 4 0.082339833   0.082339833
6

 0.821136249

iii) Fourth order

k1  0.1

k2  0.095

 1 
k3  hf  x0  h, y0  k2   0.1 1   0.0475  0.09525
1
 2 2 

k4  hf x0  h, y0  k3   0.1 1   0.09525  0.090475

 y1  1 
1
k1  2k2  2k3  k4 
6

 0. 9048375

Again taking ( x1 , y1 )  (0.1,0.9048375) and repeating the process, we get

k1  (0.1) (-0.9048375)  - 0.09048375

80
  1 
k2  0.1 f  x0  h, y0  k1 
1
 2 2 

 (0.1) (-(0.9048375 - 0.045241875))  - 0.085959562

 1 1 
k3  (0.1) f  x0  h, y0  k2 
 2 2 

  
 0.1   0.9048375  0.085959562   0.086185771
1
  2 

k4  hf x1  h, y1  k3 

 (0.1)(-0.9048375 0.086185771  - 0.0818655172

 y2  y1 
1
k1  2k2  2k3  k4 
6

 0.9048375
1
 0.09048375 2 0.085959562 0.086185771
6
  0.081865172  0.818730902

x 2nd order 3rd order 4th order exact y  e  x

0.1 0.905 0.904833333 0.9048375 0.904837418

0.2 0.819025 0.821136249 0.818730902 0.818730753

As we can see from the table fourth order values are closer to the exact values.

 y   xy 2 , y 0   1 by taking h  0.1
dy
Example 3. Compute y (0.3) given
dx

using Runge-kutta method of fourth order.

Solution

y    y  xy 2  f ( x, y )  ( y  xy 2 )

81
 x0  0, y 0  1, h  0.1, x1  0.1, x 2  0.2, x3  0.3

Now k1 = (0.1) (-1) = -0.1

  (0.1)     0.1 2  0.1   

k2  0.1  1     1      
  2   2   2  
  1


 0.1 0.95  0.050.95
2

 - 0.0995125

 h 1 
k3  hf  x0  , y0  k2 
 2 2 

 0.1 f 0.05,0.95024375


 0.1 0.95024375 0.050.95024375
2

 - 0.09953919

k4  hf x0  h, y0  k3 

 0.1 f 0.1,0.900460809

 0.098154377

 y1  1 
1
 0.1  2 0.0995125  2 0.09953919   0.098154377
6

y0.1  0.900623707

Now again taking ( x1 , y1 ) in place of ( x0 , y 0 ) , we get

k1  hf x1, y1   0.1 f 0.1,0.900623707

 0.098173601

82
  h k 
k2  hf  x1  , y1  1 
 2 2

 0.1 f 0.15,0.851536906

 0.096030417

 h k 
k3  hf  x1  , y1  2 
 2 2


 0.1 0.852608498 0.150.852608498
2

 0.096164968

k4  hf x1  h, y1  k3 

 0.1 f 0.2,0.804458738

 - 0.093388951

 y2  0.900623707
1
 0.098173601 2 0.096030417  2
6
 0.096164968   0.093388951 

 0.804631486

Again, using ( x2 , y 2 ) = (0.2, 0.804631486) we find y 3  y (0.3)

k1  hf x2 , y2 


 0.1 0.804631486 (0.2)(0.804631486) 2 
= - 0.093411785

 
  (0.1) f 0.25,0.757925593
h k
k2  hf  x2  , y2  2
 2 2 

= - 0.090153839

83
  k 
k3  hf  x2  h, y2  2   0.1 f 0.25,0.759554566
 2

 -0.090378535

k4  hf x2  h, y2  k3   0.1 f 0.25,0.714252951

= - 0.084179227

 y3  0.804631486
1
 0.093411785 2 0.090153839  2
6
 0.090378535   0.084179227

 0.714855526

Example 4. Using Runge-kutta method of fourth order, find y (0.8)

if y   y  x , y (0.6)  1.7379
2

Solution

Here x0  0.6, y0  1.7379, h  0.1, x1  0.7, x2  0.8 and f ( x, y )  y  x

2

 
k1  hf x0 , y0   0.1 1.7379  0.62  0.13779

 h k 
k2  hf  x0  , y0  1 
 2 2

 0.10.65,1.806795  0.1384295

 h k 
k3  hf  x0  , y0  2 
 2 2

 
 0.1 1.80711475 0.652  0.38461475

k4  hf x0  h, y0  k3 

84
 y 0.7   1.7379 
1
0.13779 20.1384295  20.138461475  0.138636147
6

= 1.876268016

Now using ( x1 , y1 ) = (0.7,1.876268016) we continue to find y (0.8)

 
k1  hf x1 , y1   0.1 1.876268016 0.7 2  0.138626801

 h k 
k2  hf  x1  , y1  1 
 2 2

 h k 
k3  hf  x1  , y1  2 
 2 2

 0.1 f 0.75,1.945422087  0.138292208


k4  hf x1  h, y1  k3 
 0.1 f 0.8,2.014560224   0.137456022

y 0.8  1.876268016
1
0.138626801 20.138308141  2
6
0.138292208  0.13746022
 2.014481936

Example 5. Using Runge-kutta method of fourth order, solve

dy y 2  x 2
 given y (0)  1 at x  0.2,0.4
dx y 2  x 2

Solution

x0  0, y0  1, h  0.2, x1  0.2, x2  0.4 and f x, y   y 2  x 2

2 2

y x

85
 1 0 
k1  hf x0 , y2   0.2   0.2
1 0 

 h k 
k2  hf  x0  , y0  1 
 2 2

 1.12  0.12 
 0.2 2   0.196721311
2 
 1. 1  0.1 

 h k 
k3  hf  x0  , y0  2 
 2 2

 1.0983606562  0.12 
 0.2   0.196711597
2 
 1.098360656  0.1 
2

k4  hf x0  h, y0  k3 

 1.1967115982  0.22 
 0.2    0.18913131
2 
 1.196711598  0.2 
2

y 0.2   1  0.2  20.196721311  20.196711597  0.18913131

1
6

= 1.195999521

To find y (0.4) we use ( x1 , y1 ) = (0.2,1.195999521)

 1.1959995212  0.22 
k1  0.2 f x1 , y1   0.2 
2 
 1.195999521  0.2 
2

= 0.189118717

 k 
k2  0.2 f  x1  , y1  1   0.2 f 0.3,1.29055888
h
 2 2

 0.179493515

86
  k 
k3  hf  x1  , y1  2   0.2 f 0.3,1.285746279
h
 2 2

 0.179347655

k4  hf x1  h, y1  k3   0.2 f 0.4,1.375347176

 0.168804528

 y 0.4  1.195999521
1
0.189118717 20.179493515  2
6
0.179347655  0.168804528
= 1.375267119

Runge-kutta method for simultaneous first order differential equations

Aim- To solve numerically the simultaneous equations

 f1  x, y, z  and  f 2  x, y, z  given the initial conditions

dy dz
dx dx

y ( x0 )  y 0 , z ( x0 )  z 0

Now starting from ( x0 , y 0 , z 0 ) the increments y and z in y and z respectively are given

by formulae

k1  hf1 x0 , y0 , z0  m1  hf2 x0 , y0 , z0 

 h k m   h k m 
k2  hf1  x0  , y0  1 , z0  1  m2  hf2  x0  , y0  1 , z0  1 
 2 2 2  2 2 2

 h k m   h k m 
k3  hf1  x0  , y0  2 , z0  2  m3  hf2  x0  , y0  2 , z0  2 
 2 2 2   2 2 2 

k4  hf1 x0  h, y0  k3 , z0  m3  m4  hf2 x0  h, y0  k3 , z0  m3 

and y 
1
k1  2k2  2k3  k4  z 
1
m1  2m2  2m3  m4 
6 6

87
 y1  y0  y and z1  z0  z

By repeating this algorithm once again we can find ( x2 , y2 , z 2 ) starting from ( x1 , y1 , z1 ) .

Example 6. Find y (0.1) , z (0.1) from the system of equations,

dy dz
 x  z,  x  y 2 given y (0)  2, z (0)  1 using Runge-kutta
dx dx

method of fourth order.

Solution

Here x0  0, y0  2, z0  1, h  0.1, f1 ( x, y, z)  x  z, f 2 ( x, y, z)  x  y
2

So k1  hf1 x0 , y0 , z0  m1  hf2 x0 , y0 , z0 

 (0.1) (0  1)  0.1  (0.1) (0 - 2 2 )  - 0.4

 h k m  h k m 
k2  hf1  x0  , y0  1 , z0  1  m2  hf2  x0  , y0  1 , z0  1 
 2 2 2  2 2 2

 (0.1)(0.05 1  (-0.2))  (0.1) (0.05 - 4.2025)

 0.085  - 0.41525

 h k m   h k m 
k3  hf1  x0  , y0  2 , z0  2  m3  hf2  x0  , y0  2 , z0  2 
 2 2 2   2 2 2 

 (0.1) (0.05  0.792375)  (0.1) (0.05 - (2.0425)2 )

 0.0842375  - 0.412180625

k4  hf1 x0  h, y0  k3 , z0  m3  m4  hf2 x0  h, y0  k3 , z0  m3 

 (0.1) (0.1  0.587819375)  (0.1) (0.1 - (2.0842375) 2 )

 0.068781937  - 0.424404595

 y1  2 
1
0.1  2(0.085  2(0.0842375)  0.068781937)
6

88
 = 2.084542823

z1  1 
1
 0.4  2 0.41525  2 0.412180625   0.424404595
6

= 0.586789025

Runge-kutta method for second order differential equation

Aim – To solve y   f x, y, y' given y ( x0 )  y 0 , y ( x0 )  y 0

Now, let y '  z and so y   z 

 y   z'  f x, y, y' and y   z are the two simultaneous equations with f1 ( x, y, z )  z

and f 2 ( x, y, z)  f ( x, y, z) itself.

Once again by applying Runge-kutta method for solving simultaneous first order differential

equations we can get the solution of the given problem.

Example 7. Given y   xy   y  0, y (0)  1, y (0)  0, find y (0.1) by using

Runge-kutta method of fourth order.

Solution

y ' ' xy' y  0  y ' '   y  xy'

Let y '  z  z '   y  xz

 z  f1  x, y, z 
dy

dx

  y  xz  f 2  x, y, z 
dz
dx

given y0  1, z0  0, y0 , h  0.1, x0  0
'

k1  hf1 x0 , y0 , z0   0.10  0 m1  hf x0 , y0 , z0   0.1

89
  h k m   h k m 
k2  hf1  x0  , y0  1 , z0  1  m2  hf2  x0  , y0  1 , z0  1 
 2 2 2  2 2 2

= (0.1) (0.05) = - 0.005 = (0.1) (0.05,1,-0.05)

 0.005 
k3  hf1  0.05,1   ,0.049875  (0.1) [-1 - (0.05)(-0.05)]
 2 

 (0.1) (-0.049875)  - 0.09975

 - 0.0049875 m3  hf2 0.05,0.9975,0.049875

k4  hf1 x1 , y0  k3 , z0  m3   (0.1) ((0.05)(0.049875)0.9975)

 (0.1) (-0.099500625)  - 0.099500625

= - 0.0099500625 m4  hf2 xi , y0  k3 , z0  m3 

 (0.1)(0.1,0.09950125,0.099500625)

 (0.1)((0.1)(0.09950125)  0.9950125)

 - 0.098506243

 y1  1 
1
0  2(0.005)  2(0.0049875)  (0.0099500625)
6

 0.995012489

Exercise

In the exercise below, unless specified use fourth order Runge-Kutta method.

1. Find y (0.2) given y   y  x, y (0)  2 taking h  0.1 .

dy
2. Obtain the value of y at x  0.2 if y satisfies  x 2 y  x, y (0)  1
dx

taking h  0.1 .

dy
3. Solve  xy for x  1.4 , taking y (1)  2, h  0.2 .
dx

90
 dy y  x
4. Solve  given y (0)  1 , to obtain y (0.2) .
dx y  x

du
5. Solve the initial-value problem  2tu 2 , u (0)  1 with h  0.2 on the interval (0,0.6) .
dt

1
6. Evaluate y (0.1), y (0.2), y (0.3) given y   ( x  1) y 2 , y (0)  1 .
2

dy y 1
7. Solve   2 , y (1)  1 for y (1.1) taking h  0.05 .
dx x x

8. Find y (0.1), y (0.2) given y   x  2 y, y (0)  1 taking h  0.1 by

i) second order ii) third order

iii) fourth order Runge-Kutta methods.

9. Solve y   xy  1 as x  0.2, x  0.4, x  0.6 given y (0)  2 taking h  0.2 .

dy dz
10. Solve the system:  xz  1,   xy for x  0.3, x  0.6, x  0.9 taking
dx dx

x  0, y  0, z  1 .

dy dz
11. Solve  x  z,  x  y , given y (0)  0, z (0)  1 for x  0.0 to 0.2 taking h  0.1.
dx dx

dy dz xy
12. Evaluate y (1.1), z (1.1) given  xyz,  , y (1)  0.5, z (1)  1 .
dx dx z

13. Using Runge-Kutta method determine x(0.1), y (0.1)

dx dy
given  xy  t , x(0)  1,  ty  x, y (0)  1 .
dt dt

14. Solve y   x( y ) 2  y 2  0 using Runge-Kutta method for x  0.2 given

y (0)  1, y (0)  0 taking h  0.2 .

15. Find y (0.1) given y   y 2 , y (0)  10, y (0)  5 by Runge-Kutta method.

16. Find y (0.1, y (0.2) given y   x 2 y   2 xy  1, y (0)  1, y (0)  0 .

91
 d 2 x tdx
17. Obtain the value of x (0.1) given   4 x, x(0)  3, x (0)  0 .
dt 2 dt

18. Compute the value of y (0.2) given y    y, y (0)  1, y (0)  0 .

Predictor-Corrector Method

The methods so far discussed are called single step methods because they use only the

information from the last step computed. But now we try to discuss multi-step methods.

 f  x, y , y  x0   y0
dy
Consider
dx

We have used Euler’s formula to solve differential equation of the form

yi 1  yi  hf ' xi , yi , i  0,1,2, ...1

, and we have improved the Euler method by

h f  xi , y i   f  xi 1 , y i 1  ...2 
1
y i 1  y i 
2

In (2), to get the value of y i 1 we require y i 1 on the right hand side.

To overcome this difficulty, we calculate y i 1 using Euler’s formula (1) and then we use it on the

right hand side of (2), to get the left hand side of (2). This y i 1 can be used further to get refined

y i 1 on the left hand side. Here, we predict the value of y i 1 from the rough formula (1) and use

in (2) to correct the value. Every time, we improve using (2). Hence equation (1) Euler’s

formula is a predictor and (2) is a corrector.

A predictor formula is used to predict the value of y at y i 1 and a corrector formula is used

to correct the error and to improve that value of y i 1 .

Milne’s Predictor- Corrector Formulae

 f  x, y , y  x0   y0 numerically.
dy
Suppose our aim is to solve
dx

92
 y1  y0 x0  h   y x1   y x1 , y2  y x0  2h , ..., yn  y x0  nh , where h is a suitable

accepted spacing, which is very small.

By Newton’s forward interpolation formula, we have

u u  1 2 u u  1u  2  3
y  y 0  u y 0   y0   y0  ...
2! 3!

x  x0
, where u   x  x0  uh
h

Changing y to y 

u u  1 2 ' u u  1u  2  3 '

y '  y0'  uy0'   y0   y0  ... (2)
2! 3!

Integrating both sides from x 0 to x 4 ,

x4 x0  4 h  ' uu  1 2 ' 

 y' dx    y0  uy0   y0  ... dx  y |xx00  4 h
'
x0 x0
 2! 

4 uu  1 2 ' 
 h  y0'  uy0'   y0  ... du
0
 2 

Since x  x0  uh  dx  hdu

 u 2 1 2 '  u3 u 2  1 3 '  u 4  
Now y4  y0  h uy0'  y0'   y0      y0   u 3  u 2   ...
 2 2 3 2 6  4  

 1  64  
 h4y0'  8y0'    8 2 y0'  3 y0' 64  64  16  ...
1
 2 3  6 

 20 8 
 h4y0'  8y0'  2 y0'  3 y0'  ...
 3 3 

 
 h4y0'  8( E  1)y0'  E  1 y0'  E  1 y0'  4 y0'  ...
20 2 8 3 14
 3 3 45 


 h 4y0'  8( y1'  y0' )  
20 '
 
8
 14 
y0  2 y1'  y0'  y3'  3 y2'  3 y1'  y0'  4 y0'  ...
 3 3 45 

93
  20 8  '  40  '  20  8  14
 h 4  8    y0   8   8  y1    8  y2  y3'   h4 y0'  ..
 3 3  3   3  3  45

8 4 8  14
 h y1'  y2'  y3'   h4 y0'  ...
3 3 3  45


4h
3

2 y1'  y2'  2 y3'    y0  ... ....3
14h 4 '
45

Taking into account only up to the third order equation, (3) gives

y4  y0 
4h
3

2 y1'  y 2'  2 y 3' 
…………… (4)
 y0 
4h
2 f1  f 2  2 f 3 
3

y   , where
14h 4 ' 14h v 
The error happened in (4) is  y0  ... and this can be proved to be
45 45

  x0 , x4  .

Since   E  1  e hD  1  hD for small values of h .

14h 5 v 
 The error is y   and hence (3) becomes
45

y 4  y0 
4h
3

2 y1'  y2'  2 y3'  45
y  
14h v 
……..(5)

In general,

yn 1  yn  3 
4h
3

2 yn'  2  yn' 1  2 yn'  
14h5 v 
45
y 1  , where 1  xn  3 , xn 1  ……(6)

This equation is called Milne’s predictor formula.

To get Milne’s corrector formula, integrate equation (2) between the limits x 0 to x 0  2h .

x0  2 h x0  2 h  ' uu  1 2 ' 

 y' dx    y0  uy0 
'
 y0  ... dx
x0 x0
 2 

94
 2 uu  1 2 ' 
y2  y0  h  y0'  uy0'   y0  ... du
0
 2 

2 u2 1  u3 u 2  
 h   y0' u  y0'    2 y0'  ... 
2
0
0
 2 2 3 2 

 18  
 h 2 y0'  2E  1y0'    2 E  1 y0'   4 y0'  ... 
2 4 1
 23  15 24 


 1
  1
 
 h 2 y0'  2 y1'  y0'  y2'  2 y1'  y0'  4 y0'  ... 
 3 90 


3

h '

h
y0  4 y1'  y2'  4 y0'  ...
90
…….…(7)

Again here taking into account only up to third order, we get

y 2  y0 
h '
3

y0  4 y1'  y2'  ...(8)

h 5 v 
y  
h 4 '
Here the error is   y0  ... and this can be proved to be 
90 90

where x0    x2 ; and thus (7) becomes

y2  y0 
h '
3

y0  4 y1'  y2'  
h 5 v 
90
y  

In general,

y n 1  y n 1 
h '
3
 
y n 1  4 y n'  y n 1 ' 
h 5 v 
90
y  2  where xn 1   2  xn 1 .

This equation is called Milne’s corrector formula.

  x  y  given
dy 1
Example 1. Find y ( 2) if y (x ) is the solution of
dx 2

y (0)  0, y (0.5)  2.636, y (1)  3.595 and y (1.5)  4.968

95
 Solution

Here x0  0, x1  0.5, x 2  1, x3  1.5, x 4  2, h  0.5

y 0  2, y1  2.636, y 2  3.595, y3  4.968

f  x, y   x  y   y'
1
2

By Milne’s predictor formula,

y n 1  y n  3 
4h
3

2 yn'  2  yn' 1  2 yn' 

 y 4  y0 
4h
3

2 y1'  y2'  2 y3 

Now, y1' 
1
x1  y1   1 0.5  2.636  1.568
2 2

y2' 
1
x2  y2   1 1  3.595  2.2975
2 2

y3' 
1
x3  y3   1 1.5  4.968  3.234
2 2

40.5
 y4  2  21.568  2.2975  23.234
3

= 6.871

Using Milne’s corrector formula, we get,

y4  y2 
h '
3

y2  4 y3'  y4' 

y4' 
1
x4  y4   1 2  6.871  4.4355
2 2

 y4  3.595 
0.5
2.2975  43.234  4.4355
3

= 6.873166667

96
Example 2. Using Milne’s method find y ( 4.4) given 5 xy   y  2  0
2

given y (4)  1, y (4.1)1.0049, y (4.2)  1.0097 and y (4.3)  1.0143 .

Solution

2  y2
y'  , x0  4, x1  4.1, x2  4.2, x3  4.3, x4  4.4, h  0.1
5x

y0  1, y1  1.0049, y2  1.0097, y3  1.0143

2  y12 2  1.00492
y  '
  0.048301267
54.1
1
5x1

2  y12 2  1.00972
y2'    0.046690757
5x2 54.2

2  y32 2  1.01432
y 
'
  0.045171884
54.3
3
5x2

By Milne’s predictor formula,

y 4  y0 
4h
3

2 y1'  y2'  2 y3' 
40.1
 1 20.048301267 0.046690757 20045171884
3

= 1.018700739

2  y42 2  1.0187007392
y4'    0.043738581
5x4 54.4

Using the corrector formula, we get

y4  y2 
3

h '
y2  4 y3'  y4' 

97
  1.0097 
0.1
0.046690757 40.045171884  0.043738581
3

= 1.018737229

Example 3. Determine the value of y (0.4) using Milne’s method given y '  xy  y , y 0  1;
2

use Taylor series to get the values of y (0.1), y (0.2), and y (0.3)

Solution

Here x0  0, x1  0.1, x2  0.2, x3  0.3, x4  0.4, y0  1, h  0.5

y'  xy  y 2  y0'  x0 y0  y02  1

y' '  xy' y  2 yy'  y0''  x0 y0'  y0  2 y0 y0'  3

y' ' '  xy' '2 y'2 yy' '2( y' )2  y0'''  10

h 2 '' h3 '''
y1  y0  hy0'  y0 y0  ...
2 3!

 1  0.11  3  0.00110  ...  1.116666667

0.01
2 6

Now y1  x1 y1  y1  1.358611111
' 2

y1''  x1 y1'  y1  2 y1 y11  4.28675926

 
y1'''  x1 y1''  2 y1  2 y1 y11'  2 y1'
2
 16.4113088

 y 2  1.11666667 0.11.358611111  
0.01 4.28675926   0.00116.4113088   ...
2 6

= 1.276696793

Once again y2'  x2 y2  y22  1.885294059

y2''  x2 y2'  y2  2 y2 y2'  6.467653362

98
  
y2'''  x2 y2''  2 y2'  2 y2 y2''  2 y2'
2
 28.68725078

 y3  1.276696793 0.11.885294059  6.467653362  0.00128.68725078

0.01
2 6

= 1.502345674

By Milne’s predictor formula

y 4  y0 
4h
3

2 y1'  y2'  2 y3' 

Since y1'  1.358611111

y2'  1.885294059

y3'  x3 y3  y32  2.707746226

40.1
 y4  1  21.358611111  1.885294059 22.70774626
3

= 1.832989424

y14  x4 y4  y42  4.093046

Now using Milne’s corrector formula,

y4  y2 
3

h '
y2  4 y3'  y4' 

 1.276696793
0.1
1.885294059 42.70774626  4.093046
3

 1.83700763

Example 4. Given
dy 1

 1  x 2 y 2 and
dx 2
 y0  1, y0.1  1.06, y(0.2)  1.12, y(0.3)  1.21 ,

evaluate y (0.4) by Milne’s predictor-corrector method.

99
 Solution

x0  0, x1  0.1, x2  0.2, x3  0.3, x4  0.4

y0  1, y1  1.06, y 2  1.12, y3  1.21, h  0.1

y' 
1
2
1  x2 y2 

y0' 
1
1  x0 2 y02  1 1  01  1
2 2 2

y1' 
1
2
 2 1
2
  2

1  x12 y12  1  0.12 1.06  0.567418

 1 1
y 3  (1  x32 ) y 32  (1  0.3 2 )(1.21) 2  0.7979345
2 2

By Milne’s predictor formula

y 4  y0 
4h
3

2 y1'  y2'  2 y3' 
40.1
 1 20.367418  0.652288 20.7979345
3

= 1.277122267

By Milne’s corrector formula

y4  y2 
h '
3

y2  4 y3'  y4'  y '
4 
 0.946003944

 1.12 
0.1
0.652288 4(0.7979345)  0.946003944
3

= 1.279667665

Once again if we use this value of y 4 , we get

y4' 
1
2
 
1  0.4 2 1.279667665  0.949778612
2

100
 So y4  1.12 
0.1
0.652288 40.7979345  0.949778612
3

= 1.279793487

If we repeat this procedure once again using this y 4 ,

y 4' 
1
2
 
1  0.4 2 (1.279793487) 2  0.949965394

By Milne’s corrector formula

y4  1.12 
0.1
0.652288 40.7979345  0.949965394
3

= 1.279799713

Continuing this process again and again, after some steps we get

y4  1.279800037

Example 5. Given y   1  y , and y (0)  0 , and

i) y (0.1) by Euler method; using that value obtained

ii) y (0.2) by Modified Euler method

iii) Obtain y (0.3) by Improved Euler method and find

iv) y (0.4) by Milne’s method.

Solution

By Euler method y1  y0  hf x0 , y0   0  0.11  0  0.1

By Modified Euler

 
y2  y1  hf  x1  , y1  hf x1 , y1 
h 1
 2 2 

  
 0.1  0.11   0.1  0.11  0.1   0.1855
1
  2 

101
By improved Euler method

h f  x2 , y2   f  x3 , y2  hf  x2 , y2 
1
y3  y2 
2

 0.1855 
0.1
1  0.1855  1  0.1855  0.11  0.1855
2

= 0.2628775

Now using Milne’s predictor formula

y 4  y0 
4h
3

2 y1'  y2'  2 y3' 

 y0 
4h
21  y1   1  y2   21  y3 
3

40.1
 0 21  0.1  1  0.1855  2(1  0.2628775
3

= 0.327966

y4'  1  y4  1  0.327966  0.672034

By Milne’s corrector formula

y4  y2 
3

h '
y2  4 y3'  y4' 

 y2 
h
3

1  y2   41  y3   y4' 
0.1
 0.1855  ((1  0.1855)  4(0.2628775)  0.672034
3

= 0.333334133

Exercise

1. Using Milne’s method, find y (0.2) given

dy
 0.2 x  0.1 y, y (0)  2, y (0.05)  2.0103, y (0.1)  2.0211, y (0.15)  2.0323
dx

102
 2. Find y (0.8) given y   y  x 2 , y (0)  1, y (0.2)  1.12186, y (0.4)1.46820, y (0.6)  1.7379 .

3. Using Runge-Kutta method of fourth order, find y at x  0.1, x  0.2, x  0.3 given

y   xy  y 2 , y (0)  1. Continue your work to get y (0.4) by Milne’s method.

1
4. Solve y   (1  x) y 2 , y (0)  1 by Taylor series method at x  0.2, x  0.4, x  0.6 and
2

hence find y (0.8) and y (1) by Milne’s method.

dy
5. If  2e x  y, y (0)  2, y (0.1)  2.010, y (0.2)  2.040, y (0.3)  2.090 , find y (0.4) and
dx

y (0.5) by Milne’s method.

6. Estimate y (0.8) and y (1) using Milne’s method correct to three decimal places, given

y   1  y 2 , y (0)  0, y (0.2)  0.2027, y (0.4)  0.4228, y (0.6)  0.6841 .

7. Solve y   x  y 2 , y (0)  1 to obtain y (0.4) by Milne’s method. Obtain the data you

require by any method you like.

dy y 1
8. Using both predictor-corrector methods, estimate y (1.4) if y satisfies   2 , and
dx x x

y (1)  1, y (1.1)  0.996, y (1.2)  0.986, y (1.3)  0.972 .

1
9. Given y   , y (0)  2, y (0.2)  2.0933, y (0.4)  2.1755, y (0.6)  2.2493 , find
x y

y (0.8) by Milne’s predictor-corrector method.

10. Compute y (0.6) by Milne’s method given y   x  y, y (0)  1 with h  0.2 . Obtain the

required data by Taylor series method.

11. Given y   3e x  2 y, y (0)  0 , find y (0.1) by Euler method; y (0.2) by Taylor series

method; y (0.3) by Runge-Kutta method and y (0.4) by Milne’s method.

103
12. Find y (0.2) by Taylor series method; y (0.4) by modified Euler method; y (0.6) by

Runge-Kutta method and y (0.8) by Milne’s method, given y   1  y 2 , y (0)  0 .

13. Given y   1  xy, y (0)  1 obtain y (0.1) by Picard’s method; y (0.2) by modified Euler

method; y (0.3) by Runge-Kutta method; y (0.4) by Milne’s predictor-corrector method.

14. Given y   2  xy 2 , y (0)  10, obtain power series by Picard’s method; using Milne’s

method, estimate and show that y (1)  1.6505, h  0.2 .

Estimate y (0.5), y (0.4) given y   x  y 2 , y (0)  1 using h  0.1 .

Boundary-Value Problem (BVP)

BVPs can be solved numerically by using either the Shooting Method or the Finite-

Difference Method (FDM). Here we consider the numerical solution of PVP using FDM. The

former method is left for the students as a reading assignment.

Some simple examples of two-point linear BVPs are:

y ''  x   f  x  y '  x   g  x  y  x   r  x  (1)

with boundary conditions

y  x0   a and y  xn   b (2)

yi 1  2 yi  yi 1
yi''  2
 O  h2  y iv  x   p  x  y  x   q  x 
h

(3)

with boundary conditions

y  x0   y '  x0   A and y  xn   y '  xn   B . (4)

Problems of the type in Eq. (3) and Eq. (4), which involve the fourth-order differential

equation, are much involved and will not be discussed here. There exist many methods for

104
solving second-order BVPs of the type in Eq. (1) and Eq. (2). Of these, the Finite-Difference

Method is popular one and will be discussed.

Finite-Difference Method (FDM)

The FDM for the solution for a two-point BVP consists in replacing the derivatives occurring

in the differential equation (and the boundary conditions as well) by means of their finite-

difference approximations and then solving the resulting linear system of equations by a

standard procedure.

To obtain the approximate finite-difference approximations to the derivatives, we proceed as

follows:

Expanding y  x  h  in Taylor’s series, we have

y  x  h   y  x   hy '  x   h2 y ''  x   h6 y '''  x   ...

2 3
(5)

from which we obtain

y  x  h   y  x  h ''
y'  x    y  x   ...
h 2

Thus we have

y  x  h  y  x
y'  x    O h (6)
h

Which is forward difference approximation for y '  x  .

Similarly, expansion of y  x  h  in Taylor’s series gives

y  x  h   y  x   hy '  x   h2 y ''  x   h6 y '''  x   ...

2 3
(7)

from which we obtain

y  x  y  x  h
y'  x    O h (8)
h

105
Which is the backward difference approximation for y '  x  .

A central difference approximation for y '  x  can be obtained by subtracting Eq. (7) from

Eq. (5). We thus have

y  x  h  y  x  h
y'  x    O  h2  (9)
2h

It is clear that Eq. (9) is better approximation to y '  x  than either Eq. (6) or Eq. (8). Again,

adding Eq. (5) and Eq. (7) we get an approximation for y ''  x  as

y  x  h  2 y  x  y  x  h
y ''  x   2
 O  h2  (10)
h

In a similar manner, it is possible to derive finite-difference approximations to higher

derivatives.

To solve the BVP defined by Eq. (1) and Eq. (2), we divide the range  x0 , xn  into n equal

subintervals of width h so that

xi  x0  ih, i  0,1, 2,..., n.

The corresponding values of y at these points are denoted by

y  xi   yi  y  x0  ih  , i  0,1, 2,..., n.

From Eq. (9) and Eq. (10), values of y '  x  & y ''  x  at the point x  xi can now be written as

yi 1  yi
yi'   O  h2 
2h

yi 1  2 yi  yi 1
and yi''  2
 O  h2 
h

106
In many applied problems; however, derivative boundary conditions may be prescribed, and

this requires a modification of the procedures described above. The following examples

illustrate the application of the FDM.

Example 1: A boundary-value problem is defined by

y ii  y  1  0, 0  x  1,

where y  0   y 1  0, and h  0.5

Use the FDM to determine the value of y(0.5). Its exact solution is given by

1  cos1
y  x   cos x  sin x  1 ,
sin1

from which, we obtain

y  0.5  0.139493927 .

Here nh  1 . The difference equation is approximated as

yi 1  2 yi  yi 1
 yi  1  0
h2

, and this gives after simplification

yi 1   2  h2  yi  yi 1  h2 , i  1, 2,..., n  1 ,

which together with the boundary conditions y0  yn  0 , comprises a system of  n  1

equations for  n  1 unknowns y0 , y1 ,..., yn .

1
Choosing h  (i.e. n  2 ), the above system becomes
2

107
  1 1
y0   2   y1  y2   .
 4 4

With y0  y2  0 , this gives

y1  y  0.5  0.142857142 .

Comparing with exact solution given above shows that the error in the computed solution is

0.00336 .

1
On the other hand, if we choose h  (i.e. n=4), we obtain the three equations:
4

31 1
y0  y1  y2  
16 16
31 1
y1  y2  y3  
16 16
31 1
y2  y3  y4  
16 16

Where y0  y4  0 . Solving the system we obtain y2  y  0.5  0.140311804 , the error in

which 0.00082 . Since the ration of the error is about 4 , it follows that the order of

convergence is h 2 .

These results show that the accuracy obtained by the finite difference method depends upon the

width of the subinterval chosen and also on the order of approximations. As h is reduced, the

accuracy increases but the number equations to be solved also increases.

Example 2: Solve the boundary-value problem

d2y
 y  0, 0  x  2,
dx 2

with y  0   0 & y  2   3, 62686 .

108
The exact solution of this problem is y  sinh x . The finite-difference approximation is given by

1
 yi 1  2 yi  yi 1   yi i 
h2

We subdivide the subinterval  0, 2 the four equal parts so that h  0.5 . Let the values of y at the

five points be y0 , y1 , y2 , y3 , & y4 . We are given that y0  0 , and y4  3.62686 .

Writing the difference equations at the three interval points (which are the unknowns), we obtain

4  y0  2 y1  y2   y1 

4  y1  2 y2  y3   y2   ii 

4  y2  2 y3  y3   y3 

, respectively. Substituting for y0 & y4 and rearranging, we get the system

 9 y1  4 y2  0 

4 y1  9 y2  4 y3  0   iii 
4 y2  9 y3  14.50744


The solution of (iii) is given in the table below.

x Computed solution of Exact value Error

y y  sinh x

0.5 0.52635 0.52110 0.00525

1.0 1.18428 1.17520 0.00908

1.5 2.13829 2.12928 0.00901

Exercise

1. Solve the boundary value defined by

y ''  y  0, y  0   0, y 1  1 ,

109
 by finite-difference method. Compare the computed solution at y  0.5 with the exact

value. Take h  0.25 .

2. Project work. Shooting Method. This is a popular method for the solution of two-point

boundary-value problems. If the problem is defined by

y ''  x   f  x  , y  x0   0 and y  x1   A ,

Then it is first transformed into the initial value problem

y '  x   z, z '  x   f  x  with y  x0   0 and z  x0   m0 , where m0 is a guess

for the value y '  x0  . Let the solution corresponding to x  x1 be Y0 . If Y1 is the value

obtained by another guess m1 for y '  x0  , then Y0 and Y1 and are related linearly. Thus,

linear interpolation can be carried out between the values  x0 , y0  and  m1 , y1  .

Obviously, the process can be repeated till we obtain the value for y  x1  is close to A .

Apply the Shooting method to solve the BVP

y ''  x   y  x  , y  0   0, y 1  1

110