lOMoARcPSD|52342570

Numerical Analysis For Engineer

Numerical Computing (Arba Minch University)

Scan to open on Studocu

Studocu is not sponsored or endorsed by any college or university

Downloaded by Yemane Fsaha ([email protected])
lOMoARcPSD|52342570

Numerical Methods For Engineers

CH APT E R O N E

BASI C CO N CE PT S I N E RRO R E ST I M AT I O N

I ntroduction

N umerical technique is widely used by scientists and engineers to solve their problems. A major
advantage for numerical technique is that a numerical answer can be obtained even when a problem
has no analytical solution. H owever, result from numerical analysis is an approximation, in general,
which can be made as accurate as desired. The reliability of the numerical result will depend on an
error estimate or bound, therefore the analysis of error and the sources of error in numerical
methods is also a critically important part of the study of numerical technique.
1.1. E rrors and Approximations in Computation
A computer has a finite word length and so only a fixed number of digits are stored and used
during computation. This would mean that even in storing an exact decimal number in its
converted form in the computer memory, an error is introduced. This error is machine
dependent and is called machine epsilon. I n general, we can say that

= − .
Approximation of E rror numbers:
 Exact number: number with which no uncertainly is associated to no approximation is taken.
E xample: 5, , . are exact numbers.

 Approximate number: There are numbers which are not exact.

E xample: = 2.7182 … , √2 = 1.41421 …, . They contain infinitely many non-
recurring digits. Therefore, the numbers obtained by retaining a few digits are called
approximation numbers.
E xample: ≈ 2.718 and ≈ 3.142
 Significant digits (Figures): are the numbers of digits used to express the number. The digits
1,2,3,…,9 are significant digits and ′0′ is also a significant figure except when it is used to
fix the decimal point or used to fill the place of discarded digits.
E xample: 5879, 0.4762 contains four significant digits, 0.00486, 0.000382 contains

three significant digits and 2.0682 contains five significant digits

Page 1

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

1.2. T ypes of E rrors

1.2.1. I nherent E rror: The inherent error is that quantity which is already present in the
statement of the problem before its solution. The inherent error arises either due to the
simplified assumptions in the mathematical formulation of the problem or due to the errors
in the physical measurements of the parameters of the problem. I nherent error can be
minimized by obtaining better data, by using high precision computing aids and by
correcting obvious errors in the data.
1.2.2. Round-off error: is the quantity, which arises from the process of rounding off numbers. I t
sometimes also called numerical error.
The round-off the number is the process of dropping unwanted digits is called round-off.

E xample: The number can be written as 0.29, 0.286, 0.2857

 N umber round-off according the following rules

a. Less than 5 in the ( + 1) place, leave the ( ) unaltered.
E xample: 8.893 to 8.89

b. Greater than 5 in the ( + 1) place, increase the ( ) digit by unity.

E xample: 5.3456 to 5.346

c. Exactly 5 in the ( + 1) place, increase the ( ) digit by unity if it is odd

otherwise, leaves it unchanged.
E xample: 11.675 to 11.68

11.685 to 11.68
1.2.3. T runcation E rror: These types of errors caused by using approximate formulae in
computation or on replace an infinite process by a finite one.
E xample: The Taylor’s series formula for finite terms

( )= ( )+ ( − ) ( )+ ( − ) ( )+⋯

1.2.4. Absolute error: is the numerical difference between the true value of a quantity and its
approximate value. Thus if is the approximate value of quantity then | − | is called
the absolute error and denoted by . Therefore = | − | . The unit of exact or unit
of approximate values expresses the absolute error.

1.2.5. Relative E rror: The relative error defined by = | |

= . Where is the

approximate value of quantity . The relative error is independent of units.

Page 2

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

1.2.6. Percentage error: The percentage error in which is the approximate value of is given

by = 100 ∗ = 100 ∗ | |
. The percentage error is also independent of units.

E xample: Find the absolute, relative and percentage errors if is rounded-off to three

decimal digits. Given = 0.005998

Solution: I f is rounded-off to three decimal places we get = 0.006. Therefore,
= −
= 0.005998 − 0.006 = −0.000002.
= = | | = 0.000002
.
= = = = 0.0033344 and
.

= = ∗ 100 = 0.33344.

Remark: I f the number is rounded to decimal places, then ∆ = (10 )

E xample: I f = 0.51 and correct to two decimal places then ∆ = 0.005 and the relative
∆ .
accuracy is given by = × 100 = × 100 ≅ 0.98%
.

E xercise:

a. An approximate value of is given by = = 3.1428571 and its true

value is = 3.1415926 then find the absolute and relative errors.

b. Find the relative error of the number 8.6 if both of its digits are correct.

c. Evaluate the sum = √3 + √5 + √7 to four significant digits and find its

absolute and relative errors.
d. Perform the following computations
I. Exactly
II. Using 3-digit chopping

III. Using 3-digit round-off

I V. Compute the relative error in parts of (I I ) and (I I I )

i. − +

ii. −

Page 3

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

1.3. General Formula for errors
I n this section, we derive a general formula for the error committed in using a certain formula for a
functional relation. Let = ( , ,…, ) ------------------------------------------------------------------ (1.1)
be the function of several variables ( , ,…, ) , and let the error in each be ∆ . Then the error ∆ in

is given by +∆ = ( + ∆ , + ∆ ,…, +∆ ) -------------------------------------------- (1.2)

Using Taylor’s series for more than two variables, to expand the R.H .S of above, we get

+∆ = ( , ,…, )+ ∆ +∆ + ⋯+ ∆ +

(∆ ) + (∆ ) + ⋯ + (∆ ) + ⋯ ---------------------------------------- (1.3)
∆
Assuming the errors ∆ , ∆ ,…., ∆ in all are small and that ≪ 1 so that the terms containing

( ∆ ) , ( ∆ ) ,…, ( ∆ ) and higher powers of ∆ , ∆ ,…., ∆ are being neglected.

Therefore, +∆ = ( , ,…, )+ ∆ +∆ + ⋯+ ∆ ------------------------- (1.4)

I mplies that, ∆ = ∆ + ∆ + ⋯+ ∆ --------------------------------------------------- (1.5)

Equation (5) represents the general formula for E rrors. I f equation (5) divided by we get relative
∆ ∆ 1 ∆ 2 ∆
error = = + + ⋯+ ------------------------------------------------------------ (1.6)
1 2

On taking modulus both of the sides, we get maximum relative error,

∆ ∆ ∆ ∆
≤ + + ⋯+

Also from equation (5), by taking modulus we get maximum absolute error,

|∆ | ≤ ∆ 1 + ∆ 2 + ⋯+ ∆ --------------------------------------------------------- (1.7)
1 2

1.3.1. E rror in Addition of N umbers

Let = ( + + ⋯+ )
∴ +∆ = ( +∆ )+( + ∆ ) + ⋯+ ( +∆ )
= ( + + ⋯+ ) + (∆ +∆ + ⋯+ ∆ )
Therefore, ∆ = ∆ 1 + ∆ 2 + ⋯+ ∆ this is an absolute error.
∆ ∆ 1 ∆ 2 ∆
D ividing by we get, = + + ⋯+ ; which is a relative error.
∆ ∆ 1 ∆ 2 ∆
Again, ≤ + + ⋯+ is the maximum relative error in the numbers.

1.3.2. E rror in Subtraction of N umbers

Page 4

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

∆
Let = − then we have ; which is a maximum relative error. Therefore, it shows

that when the given numbers are added then the magnitude of absolute error is the sum of
the magnitudes of the + ∆ = ( + ∆ ) −( +∆ )
or +∆ = ( − ) + (∆ −∆ )
∆ ∆ 1 ∆ 2
∴ Absolute error is given by ∆ = ∆ −∆ and Relative error is = − .

Therefore, on taking modulus of relative errors and absolute errors to get its maximum
value, we have | ∆ | ≤ | ∆ 1 | + | ∆ 2 | which is the maximum absolute error and
∆ ∆ 1 ∆ 2
≤ + which gives the maximum relative error in subtraction of numbers.

1.3.3. E rror in Product of N umbers

Let = × × × …× then using general formula for error
∆ ∆ 1 ∆ 2 ∆
∆ = ∆ 1 + ∆ 2 + ⋯+ ∆ We have = + + ⋯+
1 2 1 2

1 … 1 1 … 1 1 … 1
N ow, = = , = = ,… , = =
1 … 1 2 … 2 …
∆ ∆ 1 ∆ 2 ∆
Therefore, = + + ⋯+
1 2

∆ ∆ ∆ ∆
⇒ = ≤ + + ⋯+ and

∆ ∆
= = ×( … )
1.3.4. E rror in D ivision of N umbers

Let = then again using general formula for errors

∆ = ∆ 1 + ∆ 2 + ⋯+ ∆
1 2

∆ ∆ 1 ∆ 2 ∆ 1 1 ∆ 2 ∆ 1 ∆ 2
We have = + = 1 × + 1 × − 1
= −
1 2 2 2 1 2
2 2

∆ ∆ ∆
Therefore, the relative maximum error is ( ) = ≤ + and absolute error is
∆
= |∆ | ≤

E xample 1: Let = , find the maximum absolute and relative errors.

Page 5

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

5 2 10 −15 2
Solution: = 3 , = 3 , = 4 and

∆ = ∆ +∆ +∆

I n general, the errors ∆ , ∆ , ∆ may be negative or positive and hence we take the absolute
values of the terms of the right side.

This gives, ( ∆ ) = ∆ + ∆ + ∆

N ow let ∆ = ∆ = ∆ = 0.001 and = = = 1 then the relative maximum error is given by

(∆ ) .
( ) . = = = 0.006

E xample 2: Evaluate the sum = √3 + √5 + √7 to four significant digits and find its absolute
and relative errors.

Solution: We have √3 = 1.732, √5 = 2.236, √7 = 2.646 hence = 6.614

= 0.0005 + 0.0005 + 0.0005 = 0.0015
The total absolute error shows that the sum is correct to three significant figures only.
.
H ence we take = 6.61 and thus = = 0.0002.
.

1.4. I nverse Problems of the theory of errors

To find the error in the function = ( , ,…, ) is to have a desired accuracy and to evaluate

errors ∆ , ∆ ,…, ∆ in , ,…, we have ∆ = ∆ +∆ + ⋯+ ∆ , Using the

principles of equal effects, which states ∆ = ∆ = ⋯= ∆ this implies that

∆ ∆ ∆ ∆
∆ = ∆ or ∆ = . Similarly we get ∆ = , ∆ = ,… , ∆ = and so on.

E xample: The percentage error in which is given by = + , is not allowed to exceed 0.2%, find

allowable error in and ℎ where = 4.5 and ℎ = 5.5 .

∆
Solution: The percentage error in = × 100 = 0.2
. . . . . × .
Therefore ∆ = × = × + =
× .

Page 6

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

. × .
∆ ∆ ∆ (∆ ) ( . )
Percentage in = × 100 = × = × = = = 0.12
( . )

Percentage error in ℎ (exercise)

E xercise:
a. I f = 3 − 6 find the in at = 1, if the error in is 0.05.

b. D efine the term absolute error in = 10.00 ± 0.05, = 0.0356 ± 0.0002

= 15300 ± 100 = 62000 ± 500 Find the M aximum value of in

i. + + +
ii. + 5 −
iii.
c. Find the product of the numbers 56.54 and 12.4 which both correct to the significant digits
given.

d. Find the quotient = , where = 4.536 and = 1.32 both and being correct to the

digit given. Find also in the result.

Page 7

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

CH APT E R T WO
N U M E RI CAL SO L U T I O N S O F N O N -L I N E AR E Q U AT I O N S
M ethod for finding the root of an equation can be classif ied as direct and I terative methods.
i. D irect method: in some cases, roots can be found by using direct/ analytical methods.

For a quadratic equation + + = 0 the roots of the equation obtained by

√ √
= and = . These are called closed form solution. Similar

formula also available for cubic and bi-quadratic polynomial equations but we rarely
remember them. For higher order polynomial equations and non-polynomials it is
difficult and in many cases impossible to get closed form solutions.
ii. I terative method: those methods also known as trial and error methods are based on
the idea of successive approximation starting with one or more initial approximation to
the value of the root to obtain the accurate solution. The above quadratic equation can
be found by iterative method in the following manner.

a. = , = 0,1,2, …

b. = , = 0,1,2,…

c. = , = 0,1,2,…

This method is the best way to solve non-linear equations.

2.1. Bisection (Bolzano) method

This is one of the simplest iterative methods and is strongly based on the property of intervals
(bracketing). Finding a root using this method, let ( ) be continuous between and .

For definiteness, let ( ) be –ve and ( ) be + ve.

Then there is a root of ( ) = 0, lying between and .
T he procedure for the Bisection method
Step1: Choose two initial guess values (approx .) and (where( < ) ) such that, ( ) ∗ ( ) < 0.

Step2: E valuate the mid-point of and by = ( + ) and also evaluate ( ).

Step3: I f ( )∗ ( ) < 0 then set = else = then apply the formula of step2.

Page 8

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Step4: stop evaluation when the difference of two successive values of obtained from step2, is numerically less

than the prescribed accuracy . = ( )

( )

( )

I f an error tolerance has been prescribed in advance, it is possible to determine the number of steps

required in the bisection method. Let be a root in ( , ) and | − |< implies < by
( ) ( )
taking logarithms (with any convenient base), we obtain > ----------------------- (2.1)

Flow chart of Bisection M ethod

Start: Given , &

= ( ); = ( )

+
= ; = ( )
2
No

Is

. < 0
Yes No
Is

− Yes Stop
= ; = = ; = <
2

Page 9

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample:

a. Using Bisection method find the real root of the equation ( ) = 3 − √1 + with the

initial interval (0,1) .

b. Find the minimum number of iteration that the bisection method requires to get the

approximate root of the given function with a maximum absolute error of 10 .

Solution:

a. ( ) = 3 − √1 + = 0 is transdental equation then,

( 0) = −1 and ( 1) = 1.643. Therefore, a root lies between 0 and 1.

1st approximation = = 0.5, ( 0.5) = 0.2837 > 0 so the root lies between 0 and 0.5.
.
2nd approximation = = 0.25, ( 0.25) = −0.3669 < 0 so the root lies between 0.25

and 0.5.
. .
3rd approximation = = 0.375, ( 0.375) = −0.0439 < 0 so the root lies between

0.375 and 0.5.

( ) ( )
b. > with = 10 then to find ,

( ) ( )
> = 16.6096, this implies that the minimum number of iteration is 17.

2.2. False Position (Regula-falsi) method

The false position method retains the main features of the Bisection method, that the root is
trapped in a sequence of intervals of decreasing size.
This method uses the point where the secant lines intersect the -axis.

( = , ( ))

( = , ( ))

Page 10

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

The secant line over the interval [ , ] is the chord between ( , ( )) and ( , ( )) . The two right
angles in the figure are similar, which mean that = .
( ) ( )

( ) ( ) ( )( )
This implies that = ( )
= − ( )
------------------------------------------------------ (2.2)
( ) ( )
then we can compute ( ) and repeat the process with the interval [ , ] , if ( ) . ( ) < 0 or to
the interval [ , ] , if ( ) . ( ) < 0.

Procedures for the false position

Step1: Choose two initial guess values (approx imations) and (where < ) such that

( )∗ ( ) < 0.
( ) ( )
Step2: F ind the nex t approx imation = ; and also evaluate ( ).
( ) ( )

Step3: I f ( )∗ ( ) < 0 set = then go to the nex t step if not, rename = and then go to the
nex t step.

( ) ( )
Step4: E valuate the successive approx imation using the formula = ( )
, = 2,3,…
( )

But before applying the formula for , ensure that whether ( ). ( ) < 0, if not rename as
and proceed.

Step5: stop the evaluation when | − |< , where is the prescribed accuracy.

E xample: Find a real root of the equation ( ) = − 3 using Regula-falsi method correct to
three decimal places.

Solution: We have ( ) = − 3 = 0 assume = 1 and = 1.5 as initial guess values then

( 1) = −0.2817 and ( 1.5) = 3.7225 therefore the root lies between = 1 and = 1.5.
( ) ( ) ( . ) ( . )
1st approximation: = = = 1.0352 and
( ) ( ) ( . ) ( . )

( )= ( 1.0352) = −0.0852 thus, the root lies between = 1.0352 and = 1.5.

The function assignment is ( )= ( 1.0352) = −0.0852 and ( )= ( 1.5) = 3.7225

( ) ( ) . ( . ) ( . )
2nd approximation: = = = 1.0456 and
( ) ( ) ( . ) ( . )

( )= ( 1.0456) = −0.0252. I n similar manner, = 1.0487, = 1.0496, and

= 1.0498. H ence the approximate root is 1.0498 correct to three decimal places.

Page 11

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

2.3. Fixed-point iteration (I teration) method

This method is also known as substitution method or method of fixed iterations.
To find the root of the equation ( ) = 0 by successive approximation, we rewrite the given
equation in the form of = ( ) now first we assume the approximate value of the root
(let ) then substitute in ( ) to have the first approximation given by = ( ).
Similarly, the second approximation given by = ( )

I n general = ( ), = 1,2,3,…------------------------------------------------------------- (2.3)

= ( )

= ( )

= ( )

= ( )

Procedure for iteration method

Step1: T ak e an initial approx imation

Step2: F ind the nex t (first) approx imation by using = ( )

Step3: F ollow the above procedure to find the successive approx imation

= ( ), = 1,2,3,…

Step4: Stop evaluation where relative error less than the prescribed accuracy .

N ote: The iteration method = ( ) is convergent if | ( ) | < 1,

When | ( )| > 1 ⟹ ( ) > 1 or ( ) < −1, then it is divergent, and then the iteration

not recommended to continue.

Page 12

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample1: Find the root of the equation = 3 − 1 correct to three decimal places

using iteration method.

Solution: H ere we have ( ) = − 3 + 1 assume = 0 and = and then

( 0) = 2 = + and = −3 + 1= − so root lies b/ n 0 and .

N ow the given equation can be rewrite as = ( + 1) = ( ) (say) then we can check

that ( )= and | ( )| < 1 in [0, ] , hence iteration method can be applied.

Take the first initial guess = 0 then

= ( )= ( 0 + 1) = 0.667

= ( )= [cos(0.667) + 1] = 0.5953

= 0.6072 → = 0.6071

N ow, & being almost the same. H ence the required root is given by 0.6071.

E xample2: Find the root of the equation ( ) = + − 1 = 0 by using iteration

method.

Solution: ( 0) = −1 & ( 1) = 1 so a root lies b/ n 0 & 1.

N ow = ( )= & ( )= − we have, | ( )| < 1 for 0 ≤ ≤1

√ ( )

H ence iterative method can be applied, take = 0.5,

We get = ( )= = 0.81649
√ .

= ( )= = 0.7419, … , = ( ) = 0.75487
√ .

Page 13

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

2.4. N ewton’s method

I t is also called N ewton-Raphson iteration and it is the general root finding method.
I f the curve = ( ) cuts the -axis at , mean is the exact root of ( ) = 0. Let be an
approximation root of ( ) .
D raw the tangent to the curve at where ( ) = .
I f this tangent meets the -axis at , then is the better approximation to the root.

Tangent line

( )=

= ( )

( ) ( ) ( )
Form the figure = ( )= , ∴ − = , ⇒ = −
( ) ( )

( ) ( )
Still better approximation value = − ( )
in general, = − ( )
, = 0,1,2,….

The other way of interpreting N ewton’s method is the Taylor series at . Thus, we could write

( )+ ℎ ( )+ ( ) + ⋯ = 0 the N ewton’s method can be found by ignoring the

!
second derivative and the above ⇒ ( )+ℎ ( )= 0

( )
⇒ℎ= −
( )

( )
⇒ − = −
( )

( )
⇒ = −
( )

Page 14

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( )
Then, generally = − ( )
, = 0,1,2,…. ----------------------------------------------------- (2.4)

And lim ⟶ = where is the desired root.

Procedures in N ewton’s method

step1: T ak e initial guess as , find ( ) and ( ).

( )
step2: F ind nex t (first) approx imation by using the formula = − ( )
.

step3: F ollow the above procedure to find the successive approx imation using the formula

( )
= − , = 1,2,3,….
( )

step4: Stop the process when | − |< where is the prescribed accuracy.

E xample: Apply N ewton Raphson method to determine a root of the equation

∗ ∗
( )= − = 0 such that | ( ) | < 10 , where is the initial approximation to the
root.

Solution: Take the initial approximation = 1, we get ( )= ( 1) = −2.17797952

( )
= − ( )
= 0.65307940 , ( ) = −0.4606

( )
= − ( )
= 0.53134337 , ( ) = −0.4180 × 10

( )
= − ( )
= 0.51790991 , ( ) = −0.4641 × 10

( )
= − ( )
= 0.51775738 , ( ) = −0.5926 × 10

( )
= − ( )
= 0.51775736 , ( ) = −0.2910 × 10

⇒ | −0.2910 × 10 | < 10

Page 15

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

2.5. Secant M ethod

The N ewton’s M ethod is defined by the equation
( )
= − ( )
, = 0,1,2,…. --------------------------------------------------------------- (2.5)
One of the drawbacks of the N ewton’s M ethod is that it involves the derivative of the
function whose zero is sought.

Secant line

= ( )

Geometrical representation of secant method

Replace ( ) in equation (1) by the difference quotient, such as

( ) ( )
( )≈ ------------------------------------------------------------------------------ (2.6)
When this replacement is made, the resulting algorithm is called Secant method.
( )
I ts formula is = − ( )

= − ( )
( ) ( )
---------------------------------- (2.7)
, , ,…
This method is also called the modification of Regula-Falsi method.
Given function ( ) and two initial points , the formula is the same with Regula-Falsi
( ) ( )
method can be written as = − ( ) ( )
( ) ( )
This implies that the Secant formula is, = − ( ) ( )
, = 1,2,3,…

I n general the Secant formula can be written as = − ( )

( ) ( ) , , ,…

Page 16

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Procedures for Secant M ethod to find the root of ( )=

Step1: Choose the interval [ , ] in which ( ) = 0 has a root, where > .
Step2: F ind the nex t approx imation of the required root using the formula
( )
= − ( ) ( )
( )
Step3: F ind the successive approx imations of the required root using the formula
= − ( )
( ) ( ) , , ,…
Step4: Stop the process when the prescribed accuracy is obtained.

Flow Chart of Secant M ethod

, , =

( − ) ( )
= −
( )− ( )

= +

| − |< Yes Stop

No

Page 17

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample: Find a real root of the equation ( ) = − 3 using Secant method correct to
three decimal places.
Solution: We have ( ) = − 3 = 0 assume = 1 and = 1.5 as initial guess
values then ( 1) = −0.2817 and ( 1.5) = 3.7225 therefore the root lies b/ n 1 and 1.5.
( ) ( ) ( . ) .
= − ( ) ( )
= 1.5 − = 1.0352 ,
. ( . )
( )= ( 1.0352) = −0.0852
( ) ( ) ( . . )( . )
= − ( ) ( )
= 1.0352 −
(
= 1.0456,
. ) .
( )= ( 1.0456) = −0.0252
( ) ( ) ( . . )( . )
= − ( ) ( )
= 1.0456 − ( ) (
= 1.0500,
. . )
( )= ( 1.0500) = 0.0005
( ) ( ) ( . . )( . )
= − ( ) ( )
= 1.0500 − ( . ) ( . )
= 1.0499,
( )= ( 1.0499) = 0.0000
Therefore, the approximate real root of the equation ( ) = − 3 is 1.0499.

E xercise
1. Find the root of log = cos correct to two decimal places using Bisection method.
2. Use the Bisection method to find a root of the equation − 4 − 9 in the interval (2,3) ,
accurate to four decimal places.
3. H ow many number of iteration is needed is to find the solution of the equation
cos( )− + 2 = 0 on the interval [ 0,1] having with an accuracy less than = 0.0001
4. Given these three equations
i. − − 10 = 0
ii. − = 0
iii. +4 = 0

D etermine the initial approximation and use those to find the roots correct to three
decimal places with the following methods

a. Regula-Falsi method
b. Secant method
c. N ewton Raphson method

5. Find the approximate value of ( 17) Using N ewton Raphson method, correct to six decimal
places starting with = 2.
6. Find the root of the equations
i. − = 1
ii. −3 − + 2= 0

by method of iteration and correct to three decimal places.

Page 18

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

CH APT E R T H RE E
SO L VI N G SYST E M O F L I N E AR AN D N O N -L I N E AR E Q U AT I O N S
3.1. D irect M ethods
Consider a system of linear equations in unknowns.
+ + ⋯+ =
+ + ⋯+ =
------------------------------------------------------------------- (3.1)
⋮
+ + ⋯+ =
Where ( , , = 1,2,3,…, ) , are the known coefficients.
( , = 1,2,3,…, ) , are the known values
( , = 1,2,3,…, ) , are unknown to be determined
I n the matrix notation, the system (3.1) can be written as
= ------------------------------------------------------------------------------------------------------ (3.2)
The matrix [ | ] is called the augmented matrix. Where is the column matrix and is the
× matrix. Therefore, the system of equation = can be directly solved in the following
cases.
i. = , the equation (3.2) become
=
=
----------------------------------------------------------------- (3.3)
⋱ ⋮
=
The solution is given by = , ≠0
, , ,…,
ii. = , the equation (3.2) become
=
+ =
+ + = --------------------------------------------------- (3.4)
⋮ ⋱ ⋮
+ + … + =
Solving the first equation and then successively solving the 2nd, 3rd and so on.
We obtain
( − ) ( − − )
= , = , = , …,

∑
= where ≠ 0, = 1,2,3,…,
This method of solving equation is called forward substitution method.
iii. = , the equation (3.2) become

Page 19

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

+ + … + =
+ … + =
+ … + =
----------------------------------------- (3.5)
⋱ ⋮
, + , =
=
Solving for the unknowns in the in the order , , …, , we get
− ,
−∑
= , = ,
, …, =

This method of solving equation is called back ward substitution method.

Therefore, matrix is solvable if it can be transformed in to any one of the
forms , .

3.1.1. Cramer’s Rule

= ⇒ = | | ⁄| | , = 1,2,3,…, ------------------------------------------------ (3.6)
| | - D eterminant of the matrix and | | ≠ 0 (non-singular)
| | - D eterminant of the matrix obtained by replacing the column of by the right
hand vector .
This process is called Cramer’s rule.
+ 2 − = 2
E xample: Solve the system 3 + 6 + = 1 using Cramer’s rule.
3 + 3 + 2 = 3
Solution: | | = 35 , | | = −13 , | | = −15 , | | = 12
= , = , =
3.1.2. I nverse M ethod
Solving system of equations using the inverse of the coefficient of matrix to find the
solutions of = multiplying both sides by
⇒ = =
First, solving =
( )
Therefore, = --------------------------------------------------------------------------------- (3.7)
E xample: Solve the system of equation
3 −3 + 4 = 0
2 −3 + 4 = 3
− + = 1
Solution:
3 −3 4 1 1 1 −1 0
= 2 −3 4 ⇒ = ⇒ = −2 3 −4
det ( ) 1
0 −1 1 −2 3 −3

Page 20

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

1 −1 0
= −2 3 −4
−2 3 −3
1 −1 0 0 −3
Therefore, = = −2 3 −4 3 = 5
−2 3 −3 1 6
−3
H ence, the solution is = = 5
6
3.1.3. Gaussian E limination M ethod
This is a method reduction of equation (3.1) in to an equivalent upper triangular system
which is then solved by back substitution method.
Consider the 3 × 3 system
+ + =
+ + = --------------------------------------------------------------------- (3.8)
+ + =
I n this first stage of elimination: multiply the first row by & respectively
and subtract from the second and third rows. We get
( ) ( ) ( )
+ =
( ) ( ) ( )
----------------------------------------------------------------------------- (3.9)
+ =
( ) ( )
⎧ = − , = − ⎫
⎪ ( ) ( ) ⎪
Where, = − , = − -------------------------------------- (3.10)
⎨ ⎬
⎪ ( )
= − ,
( )
= − ⎪
⎩ ⎭
( )
Second stage of elimination: multiply the first row in (3.4) by ( ) and subtract from

the second row in (3.9). We get

( ) ( )
= ---------------------------------------------------------------------------------------- (3.11)
Collecting the equations from (3.8), (3.9) and (3.11) we obtain
( ) ( ) ( ) ( )
+ + =
( ) ( ) ( )
+ = --------------------------------------------------------------- (3.12)
( ) ( )
=
( ) ( ) ( )
Where, = , = , , = 1,2,3,…
The system (3.11) is an upper triangular system can be solved by back substitution method
[ ⁄ ] ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯ [ ∕ ]
( ) ( ) ( )
The elements , & which have been assumed to be non-zero are called pivot
elements.
Therefore, the elimination procedure describe above is called Gauss-elimination method.

Page 21

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Partial pivoting
I n the 1st stage elimination: Search an element from the first column which is largest element
in magnitude and brought as the first pivot by interchanging the first equation with the
equation having largest element in magnitude, and the same for the 2nd stage elimination
among the − 1 elements leaving the first equation and the same for the 3rd stage
elimination among the − 2 elements leaving the first and second equations … and so on.
Additionally, to reduce round-off error; it is often necessarily to perform row interchanges
even when the pivot elements are not zero.
E xample: Using the four decimal places computer, solve
0.729 + 0.81 + 0.9 = 0.6867
+ + = 0.8338
1.331 + 1.21 + 1.1 = 1.0000
I ts exact solution rounded to four decimal places, is
= 0.2245 , = 0.2814 & = 0.3279
Solution without pivoting
0.729 0.81 0.9 0.6867
1 1 1 0.8338
1.331 1.21 1.1 1.0000

= = 1.372 and ← −
= = 1.826 & ← −

0.729 0.81 0.9 0.6867

0.0 −0.1110 −0.2350 −0.1084
0.0 −0.2690 −0.5430 −0.2540

( )
= ( ) = 2.423 and ← −

0.729 0.81 0.9 0.6867

0.0 −0.1110 −0.2350 −0.1084
0.0 0.0 0.0264 0.0087
Using back substitution the solution becomes
= 0.2251 , = 0.2790 & = 0.3295
Solution with pivoting
To interchange the rows & we will use the notation ↔
0.729 0.81 0.9 0.6867
1 1 1 0.8338
1.331 1.21 1.1 1.0000

Page 22

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

↔ , = = 0.7513 and ← −
= = 0.5477 ← −

1.331 1.21 1.1 1.0000

0.0 0.0909 0.1736 0.0825
0.0 0.1473 0.2975 0.1390

( )
↔ , = ( ) = 0.6171 and ← −

1.331 1.21 1.1 1.0000

0.0 0.1473 0.2975 0.1390
0.0 0.0 −0.010 −0.0033
The solution by back substitution is = 0.2246 , = 0.2812 & = 0.3280
Comparing this solution with the exact solution rounded to four decimal places, we observe
that this is much accurate solution than the solution without pivoting.
PI T FAL LS OF E L I M I N AT I O N M E T H O D S
All nonsingular systems of linear algebraic equations have a solution. I n theory, the solution
can always be obtained by Gauss elimination. H owever, there are two major pitfalls in the
application of Gauss elimination (or its variations)
a. Round-O ff E rrors
Round-off errors occur when exact infinite precision numbers are approximated by finite
precision numbers.
E xample: Consider the following system of linear algebraic equations
0.0003 + 3 = 1.0002
+ = 1
0.0003 3 | 1.0002
Solve by Gauss-elimination. Thus, − ⁄0.0003
1 1 |1
0.0003 3 | 1.0002
⇒ = 0.3333 …
0 −9999 | −3333
We can see the difference using table

2 0.33 34
3 0.333 4
4 0.3333 1
5 0.33333 0.7
6 0.333333 0.67
7 0.3333333 0.667
8 0.33333333 0.6667
The exact solution of this problem is = and = .
Round-off errors can never be completely eliminated. H owever, they can be minimized
by using high precision arithmetic and pivoting.

Page 23

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

b. System Condition
A system of equations = is said to be ill-conditioned or unstable if it is highly sensitive to
small changes in and i.e., small change in or causes a large change in the solution of
the system. On the other hand if small changes in and give small changes in the solution,
the system is said to be stable, or, well conditioned. Thus in an ill-conditioned system, even the
small round off errors effect the solutions very badly. Unfortunately it is quite difficult to
recognize an ill-conditioned system.
E xample: Consider the following two almost identical systems
− = 1 − = 1
and
− 1.00001 = 0 − 0.99999 = 0
The respective solutions are: (100001,100000) and (– 99999,– 100000)
Obviously the two solutions differ very widely as the result of small change in the value of .
Therefore, the system is ill conditioned.

3.1.4. M atrix D ecomposition M ethod

i. decomposition (Factorization) method
= is equivalent to solving = if = & =
& Symbolically can be represent as
0 0 … … 0 … …
⎡ ⎤ ⎡ 0 … … ⎤
0 … … 0
⎢ ⎥ ⎢ 0 0 ⎥
0 … 0
= ⎢ ⎥ and = ⎢ ⎥-- (3.13)
⎢ ⋮ ⋱ ⋮ ⎥ ⋮ ⎢ ⋮ ⋱ ⎥
⎢ ⋮ ⋱ ⋮ ⎥ ⋱ ⋮ ⎢ ⋮ ⎥
⎣ … … ⎦ 0 … … 0 ⎣ 0 ⎦
To produce unique solution it is convenient to choose either = 1 or = 1
a. = 1 the method is called D oolittle’s method
b. = 1 the method is called Crout’s method

E xample: Consider the equations

+ + = 1
4 +3 − = 6
3 + 5 + 3 = 4

Solution: We can choose = 1 (Crout’s method) and write as

1 1 1 0 0 1
4 3 −1 = 0 0 1
3 5 3 0 0 1

= + +
+ + +

Page 24

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

On computing the corresponding elements, we obtain

= 1, = 4, = 3, = −1 , = 2, = −1 and

= 1, = 1, = 5

Using forward substitution we get = 1, −2,

Using backward substitution we get = 1, ,

Choose = 1 (D oolittle’s method) (E xercise)

ii. Cholesky method

For a symmetric, positive definite matrix (thus = , > 0 ≠ 0)
we can rewrite = in to = , ⇒ = and also from the LU
decomposition = , thus = (but cannot impose conditions on the main
diagonal entries). The popular method of solving = based on this factorization
is called Cholesky’s method. I n terms of the entries of = the formulas are
= √
⎧
⎪
⎪ = = 2,3,…,
---------------------------- (3.14)
⎨ = −∑ = 2,3,…,
⎪
⎪
= −∑ = + 1,…, ; ≥ 2
⎩
If is symmetric but not positive definite, this method could still be applied, but
then leads to a complex matrix , so that the method becomes impractical.
E xample: For this system of equations solve using Cholesky’s method
4 +2 + 14 = 14
2 + 17 − 5 = −101
14 − 5 + 83 = 155
Solution: From the form of the factorization
4 2 14 0 0
⇒ = ⇒ 2 17 −5 = 0 0
14 −5 83 0 0
× = 4 ⇒ = √4 = 2
× = 2 ⇒ = 1
× = 14 ⇒ = 7
× = 2 ⇒ = 1
× + × = 17 ⇒ ( ) = 16 ⇒ = 4
× + × = −5 ⇒ 4 = −12 ⇒ = −3
× = 14 ⇒ = 7

Page 25

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

× + × = −5 ⇒ 4 = −12 ⇒ = −3
× + × + × = 83 ⇒ ( ) = 25 ⇒ = 5
N ow we have to solve = , that is,
2 0 0 14
1 4 0 = −101
7 −3 5 155
7
Solving the system, then the solution is = −27
5
At the second step, we have to solve = = , that is,
2 1 7 7
0 4 −3 = −27
0 0 5 5
3
Solving the system using backward substitution, then the solution is = −6
1

3.2. I ndirect (I teration) M ethods

All the previous methods seen in solving the system of simultaneous algebraic linear equations
are direct methods. N ow we will see some indirect methods or iterative methods. This iterative
method is not always successful to all systems of equations. I f this method is to succeed, each
equation of the system must possess one large coefficient and the large coefficient must be
attached to a different unknown in that equation. This condition will be satisfied if the large
coefficients are along the leading diagonal of the coefficient matrix. When this condition is
satisfied, the system will be solvable by the iterative method.
This system in equation (3.1) will be solvable by this method if
| |> | | + | | + ⋯+ | |
⎧ | ⎫
|> | | + | | + ⋯+ | |
(D iagonally dominant)
⎨ ⋮ ⎬
⎩| |> | |+ | | + ⋯+ , ⎭
I n other words, the solution will exist (iterating will converge) if the absolute values of the
leading diagonal elements of the coefficient matrix of the system = are greater than the
sum of absolute values of the other coefficients of that row. The condition is sufficient but not
necessary.
Under the category of iterative method, we shall describe the following two methods:
i. Jacobi’s method (Gauss-Jacobi method)
ii. Gauss-Seidel method

Page 26

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

3.2.1. Jacobi I teration M ethod

Solving of the system of equations, we assume that the quantities in (3.1) are pivot
elements.
The equations (3.1) may be written as
= −( + + ⋯+ )+
⎫
= −( + + ⋯+ )+
------------------------------------ (3.15)
⋮ ⎬
= − + + ⋯+ , + ⎭
The Gauss-Jacobi iteration method may now be defined as
( ) ( ) ( ) ( )
= − + + ⋯+ − ⎫
( ) ( ) ( ) ( ) ⎪
= − + + ⋯+ − -------------------------- (3.16)
⋮ ⎬
( ) ( ) ( ) ( ) ⎪
= − ( + + ⋯+ , − )⎭
= 0,1,2,…process is continued till convergence in secure.
Where ( ) = ( )
− ( ) is the error in approximation.
N ote: I n the absence of any better estimates, the initial approximations are taken as
( ) ( ) ( )
= 0, = 0 ,…, = 0
E xample: Solve the system of Equations

4 + + = 2
+ 5 +2 = −6 Using Jacobi iteration method
+ 2 +3 = −4

Take initial approximation as ( )

= [0.5, −0.5, −0.5] and perform three iterations in each case.
( ) ( ) ( )
⎧ = − + − 2 = 0.75
⎪ ( ) ( ) ( )
Solution: = − +2 + 6 = −1.1
⎨
⎪ ( )
= −
( )
+2
( )
+ 4 = −1.1667
⎩
( ) 1 ( ) ( )
⎧ = − + − 2 = 1.0667
⎪ 4
( ) 1 ( ) ( )
= − +2 + 6 = −0.8833
⎨ 5
⎪ ( ) 1 ( ) ( )
⎩ = − +2 + 4 = −0.8500
3
( ) 1 ( ) ( )
⎧ = − + − 2 = 0.9333
⎪ 4
( ) 1 ( ) ( )
= − + 2 + 6 = −1.0733
⎨ 5
⎪ ( ) 1 ( ) ( )
⎩ = − + 2 + 4 = −1.1000
3

Page 27

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

3.2.2. Gauss-Seidel M ethod
We write the Gauss-Seidel method for the equation (3.1)
( ) ( ) ( ) ( )
= − + + ⋯+ −
( ) ( ) ( ) ( )
= − + + ⋯+ −
------------------------------ (3.17)
⋮
( ) ( ) ( ) ( )
= − ( + + ⋯+ , − )
Which may be rearrange in the form of
( ) ( )
= −∑ +
( ) ( ) ( )
+ = −∑ + ------------------------------------------------------- (3.18)
⋮
( ) ( ) ( )
+ + ⋯+ , =
This process is continued till the values of , ,…, are to the desire degree of accuracy.
The error in Gauss-Seidel method can be evaluated by ( ) = ( )
− ( ).
N ote: Gauss–Seidel method is a modification of Jacobi method. The convergence is Gauss–Seidel
method is more rapid than in Jacobi Method.
N ote: I n the absence of any better estimates, the initial approximations are taken as
( ) ( ) ( )
= 0, = 0 ,…, = 0

E xample: Solve by Gauss-Seidel method of iteration the equations

10 + + = 12
2 + 10 + = 13
2 + 2 + 10 = 14
Solution: from the given equation we have
( ) ( ) ( )
⎧ = 12 − −
⎪ ( ) ( ) ( )
= 13 − 2 − Take initial approximation as ( )
= [0,0,0]
⎨
⎪ ( )
= 14 − 2
( )
−2
( )
⎩
The first iteration becomes
( ) 1 ( ) ( ) 12
⎧ = 12 − − = = 1.2
⎪ 10 10
( ) 1 ( ) ( ) 13 − 2(1.2)
= 13 − 2 − = = 1.06
⎨ 10 10
⎪ ( ) 1 ( ) ( ) 14 − 2( 1.2) − 2(1.06)
⎩ = 14 − 2 −2 = = 0.948
10 10
The second iteration we have
( ) 1 ( ) ( ) 12 − 1.06 − 0.948
⎧ = 12 − − = = 0.9992
⎪ 10 10
( ) 1 ( ) ( ) 13 − 2( 0.9992) − 0.948
= 13 − 2 − = = 1.00536
⎨ 10 10
⎪ ( ) 1 ( ) ( ) 14 − 2( 0.9992) − 2(1.00536)
⎩ = 14 − 2 −2 = = 0.9990988
10 10

Page 28

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

This iteration process is continued. The results are tabulated as follows correcting to four decimal
places. (where is the number of iteration)
( ) ( ) ( )

0 0.0000 0.0000 0.0000

1 1.2000 1.0600 0.9480
2 0.9992 1.0054 0.9991
3 0.9996 1.0010 1.0010
4 1.0000 1.0000 1.0000
5 1.0000 1.0000 1.0000

3.3. E igenvalue Problem

A square matrix takes the set of -dimensional vectors into itself, which gives a linear function
from ℝ to ℝ . I n this case, certain non-zero vectors might be parallel to , which means
that a constant exists with = . For these vectors, we have ( − ) = 0. There is a
close connection between these numbers and the likelihood that an iterative method will
converge. We will consider this connection in this section.
Given a square × matrix
…
…
⋮ ⋮ ⋮ -------------------------------------------------------------------------------- (3.19)
…
Consider the vector equation = ⇒ ( − ) = 0, ≠ 0 where is a scalar quantity. I n
other words, we seek all those values of and non-zero vectors of such that the above
equation is satisfied. All those values of that satisfy the above equation are called the eigen-
values of and their corresponding vectors are called eigenvectors.
Recall that the Eigenvalues and eigenvectors are defined as the solutions to ( − ) = 0. This
is a homogeneous system with a coefficient matrix − .
Where ( ) = det ( − ) = 0 --------------------------------------------------------------------- (3.20)
is degree polynomial of , is called the characteristic polynomial of the matrix . Thus, we
can compute the eigenvalues and eigenvectors of a given matrix as follows:
Step 1: Compute the characteristic polynomial of .
Step 2: Find all roots for the polynomial obtained in Step 1 and label them , , …, .
Step 3: Find corresponding eigenvectors , ,…, by solving the homogeneous system
( − ) = 0, = 1,2,…, --------------------------------------------------------------------- (3.21)
E xample: Consider a 3 × 3 matrix
1 2 3
= 0 5 6
0 6 5

Page 29

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Solution:
Step 1: Compute the characteristic polynomial
1− 2 3
det ( − )= 0 5− 6 = (1 − ) [ (5 − ) − 36]
0 6 5−
= (1 − )( −1 − )(11 − )
Step 2: The Eigen-values are = 11, = 1, = −1.
Step 3: Compute corresponding eigenvectors
For = 11, we have
1− 2 3 −10 2 3
( − ) = 0 5− 6 = 0 −6 6 = 0
0 6 5− 0 6 −6
1
⇒ = = 2
2
For = 1, we have
1− 2 3 0 2 3
( − ) = 0 5− 6 = 0 4 6 = 0
0 6 5− 0 6 4
1
⇒ = = 0
0
For = −1, we have
1− 2 3 2 2 3
( − ) = 0 5− 6 = 0 6 6 = 0
0 6 5− 0 6 6
1
⇒ = = 2
−2

3.3.1. T he Power M ethod

Consider the linear eigenvalue problem:
= ----------------------------------------------------------------------------------------------- (3.22)
( )
The power method is based on repetitive multiplication of a trial eigenvector by matrix
with a scaling of the resulting vector , so that the scaling factor approaches the largest
eigenvalue and the scaled vector approaches the corresponding eigenvector . The
power method and several of its variations are presented in this section.
a. T he D irect power method
When the largest (in absolute value) eigenvalue of is distinct, its value can be found
using an iterative technique called the direct power method. The procedure is as follows:
1. Assume a trial value ( ) for the eigenvector . Choose one component of
to be unity. D esignate that component as the unity component.

Page 30

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( )
2. Perform the matrix multiplication: = ( )
3. Scale ( ) so that the unity component remains unity:
( ) ( ) ( )
=
( )
4. Repeat steps 2 and 3 with = . I terate to convergence. At convergence,
the value is the largest (in absolute value) eigenvalue of , and the vector
is the corresponding eigenvector (scaled to unity on the unity component).
The general algorithm for the power method is as follows:
( ) ( ) ( ) ( )
= = ------------------------------------------------------------- (3.23)
When the iterations indicate that the unity component could be zero, a different unity
component must be chosen. The method is slow to converge when the magnitudes (in
absolute value) of the largest eigenvalues are nearly the same. When the largest
eigenvalues are of equal magnitude, the power method, as described, fails.
E xample: Find the largest (in absolute value) eigenvalue and the corresponding
eigenvector of the given matrix
8 −2 −2
= −2 4 −2
−2 −2 13
Solution: Assume ( ) = [ 1.0 1.0 1.0] . Scale the third component to unity.
8 −2 −2 1.0 4.00 0.444444
( )
= −2 4 −2 1.0 = 0.00 , ( )
= 9.00, ( )
= 0.000000
−2 −2 13 1.0 9.00 1.000000
8 −2 −2 0.444444 1.555555
( )
= −2 4 −2 0.000000 = −2.888888 ,
−2 −2 13 1.000000 12.111111
0.128440
( )
= 12.111111, ( )
= −0.238532
1.000000
To generalize in table

0 1.000000 1.000000 1.000000

1 9.000000 0.444444 0.000000 1.000000
2 12.111111 0.128440 −0.238532 1.000000
3 13.220183 −0.037474 −0.242887 1.000000
4 13.560722 −0.133770 −0.213602 1.000000
5 13.694744 −0.192991 −0.188895 1.000000
… … … … … .... ……………… ………………. ……………… ………………
29 13.870583 −0.291793 −0.143499 1.000000

The results of the first two iterations presented above and subsequent iterations are
presented in Table. These results were obtained on a 13-digit precision computer. The
iterations were continued until changed by less than 0.000001 between iterations. The
final solution for the largest eigenvalue, denoted as , and the corresponding
eigenvector is

Page 31

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

= 13.870584 and = [ −0.291794 −0.143499 1.000000]

This problem converged very slowly (30 iterations), which is a large number of iterations
for a 3 × 3 matrix.

b. T he I nverse power method

When the smallest (in absolute value) eigenvalue of matrix is distinct, its value can
found using a variation of the power method called the inverse power method. Essentially,
this involves finding the largest (in magnitude) eigenvalue of the inverse matrix ,
which is the smallest (in magnitude) eigenvalue of matrix . Recall the original eigen
problem:
= M ultiply by gives = = after rearranging this gives an
eigenproblem for .
Thus = = ----------------------------------------------------------- (3.24)
The eigenvalues of matrix , that is, , are the reciprocals of the eigenvalues of
matrix . The eigenvectors of matrix are the same as the eigenvectors of matrix .
The power method can be used to find the largest (in absolute value) eigenvalue of
matrix , . The reciprocal of that eigenvalue is the smallest (in absolute value)
eigenvalue of matrix . I n practice the LU method is used to solve the inverse
eigenproblem instead of calculating the inverse matrix . The power method applied
( ) ( )
to matrix is given by = --------------------------------------------- (3.25)
( ) ( ) ( )
M ultiply eq. (3.25) by gives = =
( )
This can be written as = ( ) ----------------------------------------------------- (3.26)
Equation (3.26) is in the standard form = , where = ( ) and = ( ) .
Thus, for a given ( ) , ( )
can be found by the D oolittle LU method. The
procedure is as follows:
1. Solve for and such that = by the D oolittle method.
( )
2. Assume . D esignate a component of to be unity.
3. Solve for by forward substitution using the equation = ( )
4. Solve for ( ) by back substitution using the equation ( )
=
( ) ( ) ( ) ( )
5. Scale so that the unity components unity. Thus, =
( )
6. Repeat steps 3 to 5 with . I terate to convergence. At convergence, =
and ( )
is the corresponding eigenvector.
The inverse power method algorithm is as follows:
( )
=
( )
=
( ) ( ) ( )
=
E xample: Find the smallest (in absolute value) eigenvalue and the corresponding
eigenvector of the matrix given by
8 −2 −2
= −2 4 −2
−2 −2 13

Page 32

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( )
Solution: Assume = [ 1.0 1.0 1.0] . Scale the first component of to
unity. The first step is to solve for and by the D oolittle method. The results are
1 0 0 8 −2 −2
= 1 0 and = 0
1 0 0
( )
Solve for by forward substitution using =
1 0 0
⎡−1 ⎤
1.0
⎢ 1 0⎥
⎢4 ⎥ = 1.0
⎢−1 −5 ⎥ 1.0
⎣4 1⎦
7
= 1.0
−1 5
= 1.0 − ( 1.0) =
4 4
−1 −5 5 15
= 1.0 − ( 1.0) − =
4 7 4 7
( ) ( )
Solve for =
by back substitution using
8 −2 −2 1.0
⎡ 7 − 5⎤ ⎡5⎤
⎢0 ⎥ ⎢ ⎥
⎢ 2 2 ⎥ = ⎢4⎥
⎢ 75 ⎥ ⎢ 15 ⎥
⎣0 0 7 ⎦ ⎣7 ⎦
= [ 1.0 − ( −2.0)( 0.5) − ( −2) (0.2) ] ⁄8 = 0.30
5 5
= − − (0.2) (7⁄2) = 0.50
4 2
15 75
= ( ) ( ) = 0.20
7 7
( )
Scale so that the unity component is unity.
0.30 1.00000
( ) ( ) ( )
= 0.50 = 0.300000 = 1.66667
0.20 0.66667
The results of the first iteration presented above and subsequent iterations are presented
in the Table below. These results were obtained on a 13-digit precision computer. The
iterations were continued until changed by less than 0.000001 between
iterations. The final solution for the smallest eigenvalue and the corresponding
eigenvector is = = = 2.508981 and
.
= [ 1.000000 2.145797 0.599712] Therefore, to generalize in table

0 1.000000 1.000000 1.000000

1 0.300000 1.000000 1.666667 0.666667
2 0.353333 1.000000 1.981132 0.603774
3 0.382264 1.000000 2.094439 0.597565
4 0.393346 1.000000 2.130396 0.598460

Page 33

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

⋮ ⋮ ⋮ ⋮ ⋮
12 0.398568 1.000000 2.145796 0.599712
13 0.398568 1.000000 2.145796 0.599712

3.4. System of N on-L inear E quations

Solving a system of nonlinear equations is a problem that is avoided when possible, customarily
by approximating the nonlinear system by a system of linear equations. When this is unsatisfactory, the
problem must be tackled directly. The most straightforward approach is to adapt the methods from
Chapter 2, which approximate the solutions of a single nonlinear equation in one variable, to apply when
the single-variable problem is replaced by a vector problem that incorporates all the variables.
The principal tool in Chapter 2 was N ewton’s method, a technique that is generally, quadratically
convergent. This is one of the technique we modify to solve systems of nonlinear equations. N ewton’s
method, as modified for systems of equations, is quite costly to apply.

3.4.1. N ewton’s M ethod

Consider the system of two non-linear equations in two unknowns
( , )= 0
----------------------------------------------------------------------------------------------------- (3.27)
( , )= 0
Let ( , ) be a suitable approximation to the root ( , ) of the system (3.27).
Let ∆ & ∆ an increment & respectively, such that ( + ∆ , + ∆ ) is the exact solution,
that is
( +∆ , +∆ )= 0
Expand in Taylor series about the point ( , ) we get
( +∆ , +∆ )= 0
⎧ ( , )+ ∆ 1
⎪ + ∆ ( , ) + ∆ +∆ ( , )+ ⋯= 0
2!
⎨ 1
⎪ ( , )+ ∆ +∆ ( , )+ ∆ +∆ ( , )+ ⋯= 0
⎩ 2!
N eglecting 2nd and higher powers of ∆ & ∆ , we obtain
( , )+ ∆ ( , )+ ∆ ( , )= 0
------------------------------------------------------ (3.28)
( , )+ ∆ ( , )+ ∆ ( , )= 0
Since = + ∆ and = + ∆ writing the equation (3.28) in matrix form, we get
( , ) ( , ) ∆ ( , )
= −
( , ) ( , ) ∆ ( , )
Or ∆ = − ( , ) --------------------------------------------------------------------------------------- (3.29)
The solution of the system (3.29) is ∆ = − ( , ) -------------------------------------------- (3.30)
−
where = =
( , ) − ( , )
and = − ( , )
∆ ( , )
Therefore equation (3.30) can be written as = −
∆ ( , )
∆ ( , )
⇒ = + = − , = 0,1,2,…------------------------------ (3.31)
∆ ( , )

Page 34

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Or ( )
= ( )
− ( )
----------------------------------------------------------------------- (3.32)
( ) ( ) ( ) ( )
where ( )
= ( )
, ( )
, = , , , ( )
The method given by (3.32) is called N ewton’s method to the system of 2 × 2 equations.
This method can be easily generalized for solving a system of equations in unknowns.
(
, ,…, ) = 0
(
, ,…, ) = 0
--------------------------------------------------------------------------------------- (3.33)
……………
( , ,…, ) = 0
Or ( ) = 0 where = [ , ,…, ] and = [ , ,…, ]
( ) ( ) ( ) ( )
If = , ,…, is an initial approximation to the solution vector then we can write
the method as ( )
= ( )
− ( )
, = 0,1,2,…------------------------------------------ (3.34)
⁄ ⁄ ,…, ⁄
⁄ ⁄ ,…, ⁄
where = is the Jacobian matrix of the functions
⋮ ,⋯ , ⋮
⁄ ⁄ ,…, ⁄ ( ( ))
( )
, ,…, evaluated at .I f the method converges, then its rate of convergence is two.
The iterations are stopped when ( )
− ( )
< , where is the error tolerance.
E xample1: Perform three iterations of the N ewton-Raphson method to solve the system of
+ + = 7
equations Take initial approximation as = 1.5 & = 0.5
+ = 9
Solution:
( )= + + −7= 0
( )= + −9= 0
( , ) ( , ) 2 + + 2
= =
( , ) ( , ) 3 3
3 −( +2 )
= where = | |= 3 (2 + ) −3 ( + 2 )
−3 2 +
We can write now the method as
3 −( + 2 ) + + −7
1
⇒ = − , = 0,1,2,…
−3 2 + + −9
Using ( , ) = (1.5,0.5) , we get
1.5 1 0.75 −2.5 −3.75 2.2675
= − =
0.5 −14.25 −6.75 3.5 −5.50 0.9254
2.2675 1 2.5691 −4.1183 1.0963 2.0373
= − =
0.9254 −49.4951 −15.4247 5.4604 3.4510 0.9645
2.0013
=
0.9987
E xample2: Solve the system by N ewton’s method
10 + sin( + ) = 1
8 − ( − )= 1
12 + = 1
Take one step form a suitable starting point. To obtain a suitable starting point, we use the
approximation sin( + ) ≈ 0 ,cos( − ) ≈ 1 & ≈ 0 for the system of equations

Page 35

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( , , ) = 10 + sin( + ) − 1 = 0
( , , )= 8 − ( − ) −1 = 0
( , , ) = 12 + −1 = 0
and obtain = , = & =
( ) ( ) ( ) 10 + cos( + ) cos( + ) 0
= ( ) ( ) ( ) = 0 8 − sin 2( − ) sin 2( − )
( ) ( ) ( ) 0 0 12 +
1 1 1 10.939373 0.939373 0
= , , = 0 8.327195 −0.327195
10 4 12
0 0 12.996530
⎡ , , ⎤
0.091413 −0.010312 −0.000260 0.342898
⎢ ⎥
= 0 0.120089 0.003023 and = ⎢ , , ⎥ = 0.027522
0 0 0.076944 ⎢ ⎥ 0.083237
⎣ , , ⎦
Using the N ewton’s method ( )
= ( )−
We obtain for = 0 ( ) = ( ) −
H ence, = 0.0689 , = 0.246443 & = 0.076929 are the first iteration solutions of the
equation using N ewton’s method.

3.4.2. General I teration M ethod

Consider the solution of the following system of equations
( , )= 0
---------------------------------------------------------------------------------------------------- (3.35)
( , )= 0
We may write this system in an equivalent form as
= ( , )
---------------------------------------------------------------------------------------------------- (3.36)
= ( , )
Let ( , ) be exact solution of this system. Therefore, ( , ) satisfies the equations
= ( , )
---------------------------------------------------------------------------------------------------- (3.37)
= ( , )
Let ( , ) be a suitable approximation to ( , ) . Then we write the general iteration method for
= ( , )
the solution of (3.34) as -------------------------------------------- (3.38)
= ( , ) , = 0,1,2,…
I f the method is converges, then lim → = and lim → =
The functions and are called the iteration functions but not all lead to converge.
− = ( , )− ( , )
Subtracting (3.38) from (3.37) we get, ------------------------ (3.39)
− = ( , )− ( , )
Let = − and − be errors in the
= iteration. I f ( , ) is a closed
approximation to the root ( , ) then we use usually test the condition, maximum absolute row sum
| ( , )| + ( , ) < 1
norm, ----------------------------------------------------------------- (3.40)
| ( , )| + ( , ) < 1
at the initial approximation ( , ).
The method can be easily generalized to a system of equations in unknouns.

Page 36

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample 1: The system of the equation
( , )= + 3 + −5 = 0
has solution (1,1) .
( , )= −3 −4 = 0
D etermine the iteration functions ( , ) and ( , ) so that the sequence of iteration obtained
= ( , )
from
= ( , ) , = 0,1,2,…
Assume ( , ) = (0.5,0.5) converge to the root, perform 5 iterations.
Solution: We write the equivalent form as
= + ( + 3 + − 5) = ( , )
are arbitrary parameters.
= + ( + 3 − 4) = ( , )
I f the maximum absolute row sum norm, we required that
| ( , )| + ( ) < 1
,
| ( , )| + ( , ) < 1
( , ) = 1 + ( 2 + 3) ( 0.5,0.5) = 1 + 4
( , )= ( 0.5,0.5) =
( , )= 2 ( 0.5,0.5) =
( , ) = 1+ 6 ( 0.5,0.5) = 1 + 3
Therefore, the condition of convergence become
|1 + 4 | + | | < 1
& any value of which satisfy this.
| | + |1 + 3 | < 1
Take = , = we obtain the iteration method
−1 −1
= ++3 + −5 = − + −5 = ( , )
4 4
−1 −1
= + + 3 −4 = + 3 −6 −4 = ( , )
6 6
Starting with ( , ) = (0.5,0.5) we get
( , ) = ( 1.1875,1.0) ( , ) = ( 0.944336,0.931641)
( , ) = ( 1.030231,1.015702) ( , ) = ( 0.988288,0.989647)
( , ) = ( 1.005482,1.003828)
E xample 2: The iterative (explicit) method for the above second example on N ewton’s method is
can be obtain a solution correct to six decimal points. Then we can explicit the method in the form
1
= [1 − sin( + )] = ( , , )
10
1
= [1 + ( − )] = ( , , )
8
1
= [1 − sin( )] = ( , , )
12
( )
Starting with the initial approximation = , , we obtain
( )
= [0.065710 ,0.246560 ,0.076397]
( )
= [0.069278 ,0.246415 ,0.076973]
( )
= [0.068952 ,0.246445 ,0.076925]
( )
= [0.068978 ,0.246442 ,0.076929]
( )
= [0.068978 ,0.246442 ,0.076929]

Page 37

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

H ence, the solution correct to six decimal places is obtained after five iterations. This implies that
using N ewton’s method can be obtaining more accurate solution with in small number of iteration.

E xercise
1. Solve the following system of equations using Gauss-Elimination method
− + = 1 + 3 +6 = 2
(a) −3 + 2 − 3 = −6 (b) −4 + 2 = 7
2 −5 + 4 = 5 3 − +4 = 9
5 + + + = 1
+ 7 + + = 12
c)
+ + 6 + = −5
+ + + = −6
2. Solve the following set of simultaneous linear equations by the matrix inverse method
2 + 3 − = 1
a. − +2 + = 8
− 3 − 2 = −13
− +3 − = 1
−3 + 5 = 2
b.
− + = 0
+ 2 − = −5
3. Solve the following set of simultaneous linear equations using the LU decomposition method
+ + = 9
i. 2 −3 + 4 = 13
3 + + 5 = 40
2 + − = 6
ii. −3 + 5 = 11
− +5 + 4 = 13
4. Solve the following set of simultaneous linear equations using Jacobi and Gauss-Seidal iteration
method
2 − + 5 = 15
a. 2 + + = 7
+ 3 + = 10
5 + 2 + = 12
b. + 4 + 2 = 15
+ 2 + 5 = 20
4 + 2 = 4
2 +8 +2 = 0
c.
2 + 8 + 2 = 0
3 + 4 = 14
5. Solve the following systems of non-linear equations by any suitable method
+ = 11
a.
+ = 7
= 3 −7
b.
= 2( + 1)
+ + =
c. + = 1
+ =

Page 38

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

CH APT E R FO U R

I N T E RPO L AT I O N
 The word interpolation denotes the method of computing the intermediate
(approximate) values of the function = ( ) for an between different - values
, ,…, at which the values of ( ) are given.
 Or, given( , ) , ( , ) , ( , ) , …,( , ) , find the values of ‘ ’ at a values of ‘ ’
in ( , ) is called interpolation.
 The process of computing the value of the function outside the given range is called
ex trapolation.
 Polynomials are most common choice of interpolation because they are easy to
evaluate, differentiate and integrate.
 I nterpolation is required in many engineering and science applications that use tabular
data as input.

The polynomial through a specific set of points may take many different forms, but all forms are
equivalent. Any form can be manipulated into any other form by simple algebraic rearrangement
The Taylor series is a polynomial of infinite order. Thus,
( )= ( )+ ( )( − )+ ( )( − ) + ( )( − ) + ⋯----------- (4.1)
! !
I t is of course, impossible to evaluate an infinite number of terms. The Taylor polynomial of degree
is defined by
( )= ( )+ ( ) ------------------------------------------------------------------------------------- (4.2)
where the Taylor polynomial ( ) , and the remainder term ( ) are given by
( )= ( )+ ( )( − )+ ( )( − ) + ⋯+ ( )( )( − ) --------- (4.3)
! !

( )= ( )( )( − ) ≤ ≤ ----------------------------------------------- (4.4)
( )!

When a polynomial of degree , ( ) is fit exactly to a set of ( + 1) discrete data points,

( , ) , ( , ) ,…,( , ) the polynomial has no error at the data points themselves. H owever, at
the locations between the data points, there is an error which is defined by
( ) = ( ) − ( ) ----------------------------------------------------------------------------------- (4.5)
I t can be shown that the error term, ( ) , has the form
( )= ( − )( − ) …( − ) ( )( ) where ≤ ≤ ----------------- (4.6)
( )!
Equation (4.13) shows that the error in any polynomial approximation of discrete data.

4.1. Finite D ifference O perators

) , = 1,2,3,…, -------------------------------------- (4.7)
Assume that we have a table of values (
of any function = ( ) where the value of being equally spaced i.e. = + ℎ.
To determine the values of ( ) or ( ) for some intermediate values of x, the following
different types of difference are found useful.

Page 39

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

4.1.1. Forward D ifference

Let = ( ) be any function given by the values , , ,…, which is takes for the
equidistant values , , ,…, of the independent variable , then
− , − , …, − are called the first differences of the function .
They are denoted by ∆ , ∆ , …, .
∆ = −
∆ = −
Therefore, we have ------------------------------------------------------------- (4.8)
⋮
∆ = −
The symbol ∆ is called the forward difference operator. The differences of the first
differences denoted by ∆ , ∆ ,…, ∆ are called the second differences, where
∆ = ∆[ ∆ ] = −2 +
∆ = ∆[ ∆ ] = −2 + ----------------------------------------------------------- (4.9)
⋮
∆ = ∆[ ∆ ] = −2 +
∆ is called second forward difference operator.
∆ = ∆ −∆ = −3 + 3 −
th
Similarly for the n order, ⋮ ------------ (4.10)
∆ = ∆ −∆

The forward difference table (T able 4.1)

∆ ∆ ∆ ∆ ∆

∆
∆
∆ ∆
∆ ∆
∆ ∆ ∆
∆ ∆
∆ ∆
∆
∆

The above table is called a diagonal difference table. The first difference term in the table
,∆ ,∆ ,∆ …, ∆ are called the leading differences.
Alternative N otation of forward difference operator is ∆ ( ) = ( + ℎ) − ( )
This follows that ∆ ( + ℎ) = ( + ℎ + ℎ) − ( + ℎ) = ( + 2ℎ) − ( + ℎ)
Similarly, ∆ ( ) = ∆[ ∆ ( ) ]
= ∆[ ( + ℎ) − ( )]
= ∆ ( + ℎ) − ∆ ( )

Page 40

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

= ( + 2ℎ) − ( + ℎ) − [ ( + ℎ) − ( )]
= ( + 2ℎ) − 2 ( + ℎ) + ( )
∆ is called the second difference of ( ) and ℎ is the step size.

4.1.2. Backward D ifference

Let =
( ) be any function given by the values , , ,…, which is takes for the
equally spaced values , , , …, of the independent variable , then
− , − , …, − are called the first backward differences of the function .
They are denoted by ∇ , ∇ , …, ∇ respectively.
∇ = −
∇ = −
Thus, we have ----------------------------------------------------------------- (4.11)
…
∇ = −
The symbol ∇ is called the backward difference operator.

The backward difference table (T able 4.2)

∇ ∇ ∇ … ∇ ∇

∇
∇
∇ ∇
∇ ⋮
∇ ⋮ ∇ ∇
⋮ … ∇
⋮ ∇ …
⋮ ⋮ ∇ ∇
∇ ∇
∇

The last difference terms in the table , ∇ , ∇ ,∇ ,…. ∇ y can be considered as the
leading differences of the backward difference.
Alternative N otation of backward difference operator is ∇ ( ) = ( ) − ( − ℎ)
We can observe that ∇ ( + ℎ) = ( + ℎ) − ( ) = ∆ ( )
∇ ( + 2ℎ) = ( + 2ℎ) − ( + ℎ) = ∆ ( + ℎ)
⋮
∇ ( + ℎ = ( +
) ℎ ) − ( + ( − 1) ℎ ) = ∆ ( + ( − 1) ℎ )
Similarly we get, ∇ ( + 2ℎ) = ∇[ ∇ ( + 2ℎ) ]
= ∇[ ∆ ( + ℎ) ]
= ∆[ ∇ ( + ℎ)]
= ∆[ ∆ ( )]

Page 41

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

= ∆ ( )
⋮
∇ ( + ℎ) = ∆ ( )

4.1.3. Central D ifference operator

We now introduce another operator known as the central difference operator to represent
the successive differences of a function in a more convenient way. The central difference
operator, denoted by the symbol δ is defined by
⁄ = −
⁄ = −
… --------------------------------------------------------------------------------- (4.12)
= −
For higher order differences
= ⁄ − ⁄

⁄ = − --------------------------------------------------------- (4.13)
…
= ⁄ − ⁄

The central difference table (T able 4.3)

⁄ ⁄

⁄ ⁄

The alternative notation of central difference operation can be written as

( )= + ℎ − ( − ℎ) where ℎ is the interval of differencing.

The values

⁄ ⁄
…
⁄ ⁄

can be consider as leading coefficients of the central difference table.

Page 42

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

4.1.4. Additional O perators

a. T he Shift operator
Let = ( ) be function of and , + ℎ, + 2ℎ, + 3ℎ, …, etc., be the consecutive
values of , then the operator E is defined as ( ) = ( + ℎ), is called shift operator . I t is
also called displacement operator . ( ) mean the operator is applied twice on ( ) i.e.
( )= [
( )]
= ( + ℎ)
= ( + 2ℎ)
…
Similarly, ( )= ( + ℎ)
And ( )= ( − ℎ)

b. T he D ifferential operator
Let = ( ) be function of and denotes the differential coefficient of with respect
to where = . we have = = ( ) . The derivative of with respect to is
( )
denoted by = = ( ).
c. T he M ean operator
I n addition to the operator ∆, ∇, , & we define the mean operator (averaging operator)
as ( )= + ℎ + ( − ℎ) .
4.1.5. Relationship of the operators
Relationship the operator and ∆: ∆ ( ) = ( + ℎ) − ( ) where ℎ is the interval of
differencing. Using the operator we can write
∆ ( )= ( )− ( )
⇒ ∆ ( ) = ( − 1) ( )
The above relation can be expressed as an identity.
∆= − 1
i.e. = ∆+ 1
Relationship the operator , and ∆: From the definition we have
( ) = ( + ℎ)
= ( )+ ( )+ ( )+ ⋯
! !

= ( )+ ( )+ ( )+ ⋯
! !
(Expanding by Taylor’s series method)
= 1+ + + ⋯ ( )
! !
( )=∴ ( )
H ence, the identity is = .
We have already proved that = ∆ + 1 and = . Apply the algorithm,

Page 43

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

We get
∆ ∆ ∆
⇒ℎ = = [ 1 + ∆] = ∆ − + − +⋯
∆ ∆ ∆
⇒ = ∆− + − + ⋯ .
Relation between & ∇: ∇ ( ) = ( ) − ( − ℎ ) = ( )− ( )
⇒ ∇= 1 −
∇=
Relation between the operators ∆ , , & : From the definition we know that
( )= + ℎ − − ℎ
i. ( )= + ℎ − − ℎ
⁄ ( )− ⁄ ( )
=
⁄ ⁄
= − ( )
⁄ ⁄
∴ = −
⁄ ( − 1) ( ) = ⁄
Further ( )= ∆ ( )
⁄
∴ = ∆
And from the mean operator result we get
⁄
= ∆

ii. ( )= + ℎ + − ℎ

⁄ ⁄ ( )
= +

1 ⁄ ⁄
∴ = +
2

iii. ∇ ( )= [ ( ) − ( − ℎ)]

= ( )− ( − ℎ)

= ( + ℎ) − ( ) = ∆ ( )

∴ ∇= ∆

and ∇ ( ) = ∇ ( + ℎ)
= ( + ℎ) − ( ) = ∆ ( )
⇒ ∇ = ∆
∴ ∇= ∇

Page 44

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample1: Construct a forward difference table for the following data

0 10 20 30

0 0.174 0.347 0.518

Solution: The difference table is

∆ ∆ ∆
0 0
0.174
10 0.174 −0.001
0.173 −0.001
20 0.347 −0.002
0.171
30 0.518

E xample2: Construct a difference table for = ( )= + 2 +1

= 1,2,3,4,5

Solution: The difference table is

∆ ∆ ∆
1 4
9
2 13 12
21 6
3 34 18
39 6
4 73 24
63
5 136

E xample3: Construct the backward difference table for = given that

10 20 30 40 50
1 1.3010 1.4771 1.6021 1.6990
and find the values of ∇ 40 ∇ 50 .

Page 45

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Solution: For the given data, backward difference table as:

∇ ∇ ∇ ∇
10 1
0.3010
20 1.3010 −0.1249
0.1761 0.0738
30 1.4771 −0.0511 −0.0508
0.1250 0.0230
40 1.6021 −0.0281
0.0969
50 1.6990
H ence, ∇ 40 = 0.0738 ∇ 50 = −0.0508

E xample4: Show that

⁄ ⁄ ( 1 + ∆) ⁄
a. + = 2+ ∆
∆ ∇
b. − = ∆+ ∇
∇ ∆

Solution:

⁄ ⁄ ⁄
a. since 1 + ∆= therefore + = + 1= 1+ ∆+ 1= ∆+ 2
∆ ∇ ( )
b. − = − = − = − = ( − )
∇ ∆ ( )

= {( 1 + ∆) − ( 1 − ∇)} = ( ∆ + ∇)
∆ ∇
H ence, − = ∆ + ∇.
∇ ∆

4.2. I nterpolation with equal interval

4.2.1. N ewton I nterpolation Formulas

a. Forward N ewton I nterpolation Formulas
Let = ( ) take the values , , , …, corresponding to the values
, , , …, of the independent variable , when the tabular points ’s are equally
spaced
( . . = + ℎ = 1,2,3,…, ) .
Let ( ) = ( ) where the polynomial ( ) is given in the following form
( )= + ( − ) + ( − )( − ) + ⋯
+ ( − )( − ) …( − ) ------------------------------ (4.14)
The constants , , ,…, can be determined as follows.
Putting = , , ,…, , successively in (4.14) to determine the values of the
coefficients , , , …, respectively.

Page 46

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Assume =
(
)= =
When =
( )= = + ( − )= + ℎ
− ∆
= =
ℎ ℎ
When =
( )= = + ( − )+ ( − )( − )
−
= + ( 2ℎ) + ( 2ℎ)( ℎ) = + 2( − )+ ( 2ℎ )
ℎ

− 2( − )− −2 + ∆
⇒ = = =
2ℎ 2ℎ 2! ℎ
⋮
∆
⇒ = , = 3,4,5,…
!
Substitute the values in (4.14) we get
∆ ∆
( )= + ( − )+ ( − )( − ) + ⋯.
!
∆
+ ( − )( − ) …( − ) ------------------------------- (4.15)
!
Setting = then − = ℎ
− − + = − = − −( − )= ℎ − ℎ = ( − 1) ℎ
Similarly, − = ( − 2) ℎ
− = ( − 3) ℎ
⋮
− = ( − ( − 1) ) ℎ
Equation (4.15) reduces to
( )
( )= + ∆ + ∆ + ⋯.
! !
( )…( ( ))
+ ∆ ------------------------------------ (4.16)
!
The above formula is called N ewton’s forward interpolation formula.
I t is useful for interpolating near the beginning of the tabular values.
The error term for the N ewton forward-difference polynomial can be obtained from the
general error term equation (4.13)
( )= ( − )( − ) …( − ) ( )( )
( )!
We know that − = ℎ, − = ( − 1) ℎ,… , − = ( − ( − 1) ) ℎ
Therefore the error term can be written as
( )= ( − 1)( − 2) …( − ) ℎ ( )( ) ---------------------- (4.17)
( )!

And also can be written as, ( )= ℎ ( )( )

+1

Page 47

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample1: Evaluate = ( )= = 0.05 using the following table

0.00 0.10 0.20 0.30 0.40
1.0000 1.2214 1.4918 1.8221 2.2255

Solution: The difference table is

∆ ∆ ∆ ∆
0.00 1.0000
0.2214
0.10 1.2214 0.0490
0.2704 0.0109
0.20 1.4918 0.0599 0.0023
0.3303 0.0132
0.30 1.8221 0.0731
0.4034
0.40 2.2255
We have = 0.00 , = 0.05 , ℎ = 0.10
− 0.05 − 0.00
∴ = = = 0.5
ℎ 0.10
Using N ewton’s forward formula
( − 1) ( − 1)( − 2)
( )= + ∆ + ∆ ∆ +
1! 2! 3!
( − 1) ( − 2)( − 3)
∆ + ⋯
4!
0.5 0.5( 0.5 − 1)
( 0.05) = 1.0000 + × 0.2214 + × 0.049 +
1! 2!
0.5( 0.5 − 1)( 0.5 − 2)
× 0.0109 +
3!
0.5( 0.5 − 1) (0.5 − 2)(0.5 − 3)
× 0.0023 + ⋯
4!
= 1.0000 + 0.1107 − 0.006125 + 0.00068125 − 0.00008984375
= 1.105166
∴ ( 0.05) ≈ 1.1051678125
To find the maximum error expected, ( ) ( ) = 32
1
( )= ( − 1)( − 2) …( − ) ℎ ( )( )
( + 1)!
1
( ) = 0.5( 0.5 − 1)( 0.5 − 2)( 0.5 − 3)( 0.5 − 4) (0.1) 32
5! .
1
( ) .= ∗ 3.28125 ∗ ( 0.1) ∗ (32 ∗ ∗ . ) = 1.947348 ∗ 10
120

Page 48

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

b. Backward N ewton I nterpolation Formulas

Let = ( ) take the values , , , …, corresponding to the values
, , , …, of the independent variable , when the tabular points ’s are equally
spaced.
( . . = + ℎ = 1,2,3,…, ) .
Let ( ) = ( ) where the polynomial ( ) is given in the following form
( )= + ( − )+ ( − )( − ) + ⋯.
+ ( − )( − ) …( − ) --------------------------------------------- (4.18)
The constants , , ,…, can be determined as follows.
Putting = , , , …, , successively in (4.18) to determine the values of the
coefficients , , , …, respectively.
Assume =
( )= =
When =
( )= = + ( − )= + ( −ℎ)
− −∇ ∇
= = =
−ℎ −ℎ ℎ
When =
( )= = + ( − )+ ( − )( − )
−
= + ( −2ℎ) + ( −2ℎ)( −ℎ) = + 2( − ) + ( 2ℎ )
−ℎ
− 2( − )− −2 + ∇
⇒ = = =
2ℎ 2ℎ 2! ℎ
⋮
∇
⇒ =
!ℎ
Substitute the values in (4.18) we get
∇ ∇
( )= + ( − )+ ( − )( − ) + ⋯.
!
∇
+ ( − )( − ) …( − ) -------------------------------------- (4.19)
!
Setting = then − = ℎ
= − − + − = ℎ + ℎ = ( + 1) ℎ
Similarly, − = ( + 2) ℎ
− = ( + 3) ℎ
⋮
− = ( + ( − 1) ) ℎ
Equation (4.19) reduces to
( )
( )= + ∇ + ∇ + ⋯.
! !
( )…( ( ))
+ ∇ --------------------------------------------------------- (4.20)
!

Page 49

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

The above formula (4.20) is called N ewton’s back ward interpolation formula.
I t is useful for interpolating near the end of the tabular values.
The error term of backward interpolation is also can be written as
( )= ( + 1)( + 2) …( + )ℎ ( )( ) ---------------------- (4.21)
( )!
E xample: The area of a circle of diameter is given for the following values.
80 85 90 95 100
5026 5674 6362 7088 7854

Find for = 105. (Extrapolation)

Solution: First the backward difference table is

∇ ∇ ∇ ∇
80 5026
648
85 5674 40
688 −2
90 6362 38 4
726 2
95 7088 40
766
100 7854
H ence, ℎ = 5 = = 105 = = 105
− 105 − 100
∴ = = = 1
ℎ 5
N ow on applying N ewton’s backward interpolation formula, we have
( + 1) ( + 1) ( + 2)
= + ∇ +∇ + ∇
1! 2! 3!
( + 1) ( + 2)( + 3)
+ ∇
4!
1 1( 1 + 1) 1( 1 + 1)( 1 + 2)
= 7854 + (766) + (40) + (2)
1! 2! 3!
1( 1 + 1)( 1 + 2)( 1 + 3)
+ (4)
4!
= 7854 + 766 + 40 + 2 + 4
= 8666

Page 50

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

c. Central D ifference Formulas

The N ewton forward- and backward-difference formulas are not appropriate for
approximating = ( ) when lies near the center of the table because neither will
permit the highest-order difference to have close to . The main advantage of central
difference formulae is that they give more accurate result than other method of
interpolation. Their disadvantages lie in complicated calculations and tedious expression,
which are rather difficult to remember. These formulae are used for interpolation near
the middle of an argument values. For the centered-difference formulas, we choose
near the point being approximated and label the nodes directly below as , ,…and
those directly above as , , …with this convention, Stirling’s formula is for Odd
centered differences are based on the averages of the centered differences at the half
points ⁄ and ⁄ . This central difference formula is useful when | | ≤ .
The Stirling’s central-difference interpolating polynomial is
( )= +1
+ ⁄ + ⁄ + + +.
1 2 2
+1 +2 +1
⁄ + ⁄ + + + ⋯------ (4.22)
3 4 4
( )( )…( ( ))
Where = and =
!
E xample: Employ Stirling’s formula to evaluate . from the given data below
( = 1+ ).

10 11 12 13 14
10 23967 28060 31788 35209 38368

Solution: From the difference table

= 10 23967
4093
= 11 28060 −365
3728 58
= 12 31788 −307 −13
3421 45
= 13 35209 −262
3159
= 14 38368
= 12.2 − 12 = 0.2 ≤ 0.5
( )
0.2 1 ( 1 1.2 0.2 ( −307)
= 31788 + 3421 + 3728) + +
1 2 2 2 2
1
1.2 ( 1 2.2 1.2 (
+ 45 + 58) + + −13)
3 2 2 4 4

Page 51

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( )
= 31788 + ( 0.1)( 7149) + ( 0.02)( −307) + ( −0.016)( 103)
+ ( −0.0016)( −13)

( ) = 31788 + 714.9 − 6.14 − 1.648 + 0.0208 = 32495

( ) 32495
( ) = 10 ⇒ = ⇒ . = = 0.32495
10 10

4.3. I nterpolation with unequal interval

The interpolation formulae derived before for forward interpolation, backward interpolation and
central interpolation have the disadvantages of being applicable only to equally spaced argument
values. So it is required to develop interpolation formulae for unequally spaced argument values
of . Therefore, when the values of the argument are not at equally spaced then we use two such
formulae for interpolation such as L agrange’s I nterpolation formula and N ewton’s divided
difference formula.
The main advantage of these formulas is, they can also be used in case of equal intervals but the
formulae for equal intervals cannot be used in case of unequal intervals.
A. L agrange’s I nterpolation formula
A straight forward approach is the use of Lagrange polynomial.
The formula used to interpolate between pairs
, ( ) , , ( ) , , ( ) , …, , ( ) is given by ( ) = ∑ ( ) where
the polynomial is given by ( )= ( )∏ ,

1, =
( )= ∏ and ( )=
0, ≠
I n general,
( )( )…( )
⎧ ( ) +
( )( )…( )
⎪ ( )( )…( )
( )= ( ) + ⋯ + ------------------------------------------------ (4.23)
⎨ ( )( )…( )
⎪ ( )
( )( )…( )
⎩ ( )( )…( )
I n short this one can be written as
( )= ( ) ( )+ ( ) ( ) + ⋯+ ( ) ( ) ------------------------------------ (4.24)

Consider the tabular interpolating points we wish to fit

( )
0 ( )
1 ( )
2 ( )
3 ( )

Page 52

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

The linear Lagrange’s interpolation polynomial is

( − ) ( − )
( )= ( )= ( ) + ( )
( − ) ( − )
The Quadratic Lagrange’s interpolation polynomial is
( )= ( )
( − )( − ) ( − )( − )
= ( ) + ( )
( − )( − ) ( − )( − )
( − )( − )
+ ( )
( − )( − )
Similarly, the Cubic Lagrange’s interpolation polynomial is
( )= ( )
( − )( − )( − )
= ( )
( − )( − )( − )
( − )( − )( − )
+ ( )
( − )( − )( − )
( − )( − )( − )
+ ( )
( − )( − )( − )
( − )( − )( − )
+ ( )
( − )( − )( − )
The Lagrange polynomial pass through each of the points used in its construction and data
are not required to be specified with in ascending or descending order. The computations
of ( ) is simple but the method is still not particularly efficient for large value of .

E xample1: Construct the linear Lagrange interpolation polynomial that pass through the
points (2,4) and (5,1) .
Solution: ( )= = ( − 5) and ( )= = ( − 2)
So that ( )= ( ) ( )+ ( ) ( )
= ( − 5)( 4) + ( − 2)( 1)
= − +6
(2,4)

(5,1)

Page 53

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample2:
a. Use the number (called nodes) = 2.75 & = 2, = 4 to find the second
Lagrange interpolating polynomial for ( ) = 1⁄ .
b. Use this polynomial to approximate ( 3) = 1⁄3.

Solution:
a. We first determine the coefficient polynomials ( ), ( )& ( ) . I n nested form
( . )( )
they are ( )= = ( − 2.75)( − 4)
( . )( )
( )( )
( )= = ( − 2)( − 4)
( . )( . )
( )( . )
and ( )= = ( − 2)( − 2.75)
( )( . )
Also, ( ) = ( 2) = 1⁄2 , ( )= ( 2.75) = 4⁄11 ( )= ( 4) = 1⁄4,
So, ( )= ∑ ( ) ( )
= ( − 2.75)( − 4) − ( − 2)( − 4) + ( − 2)( − 2.75)
= − +
b. An approximate to ( 3) = 1⁄3 is
9 105 49
( 3) ≈ ( 3) = − + ≈ 0.32955
22 88 44

B. N ewton’s divided difference interpolation formula

The following two disadvantages of Lagrangian polynomial method lead us to develop a new
method for the interpolation. They are:
a. it involves more arithmetic operations; and
b. We essentially need to start over the computation if we desire to add or subtract a point
from the set of data.
T he basic idea of D ivided D ifference:
The degree polynomial can be written in a special form:
+ ( − ) + ( − )( − ) + ⋯
( )= ------------------------------------ (4.25)
+ ( − )( − ) …( − )
The key idea is to find , , , …, so that interpolates the given
data , ( ) , , ( ) , , ( ) , …, , ( ) .
A divided difference is defined as the difference in the function values al two points, divided
by the difference in the values of the corresponding independent variable.
Thus, = [ ]= ( )
The first divided difference at a point is defined as
( )− ( )
= [ , ]=
−
The second divided difference is given as

Page 54

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

[ , ]− [ , ]
= [ , , ]=
−
The third divided difference is given as
[ , , ]− [ , , ]
= [ , , , ]=
−
[ , ,…, ] [ , ,…, ]
I n general, = [ , , ,…, ]=
The N ewton’s D ivided D ifference table can be construct as (Table 4)
( ) [ , ] [ , , ] [ , , , ] [ , , , , ]
( )
[ , ]
( ) [ , , ]
[ , ] [ , , , ]
( ) [ , , ] [ , , , , ]
[ , ] [ , , , ]
( ) [ , , ]
[ , ]
( )

The advantage of the above method is that there is no need to start all over again if their
additional pairs of data are added. We simply need to compute additional divided
differences.
Since order polynomial interpolation of a given ( + 1) pairs of datais unique, thus the
above polynomial and Lagrangian polynomial are exactly the same.

E xample: Use N ewton’s D ivided D ifference formula and evaluate (3.0) , given
3.2 2.7 1.0 4.8 5.6
22.0 17.8 14.2 38.3 51.7

Solution:
( ) [ , ] [ , , ] [ , , , ] [ , , , , ]
3.2 .
.
2.7 17.8 .
2.118 − .
1.0 14.2 2.012 .
6.342 0.0865
4.8 38.3 2..63
16.750
5.6 51.7

The third degree polynomial fitting all points from = 3.2 to = 4.8 is given by
( )= + ( − )+ ( − )( − ) + ( − )( − )( − )

Page 55

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( ) = 22 + 8.4( − 3.2) + 2.856( − 3.2)( − 2.7)

− 0.528( − 3.2)( − 2.7)( − 1)
The fourth degree polynomial fitting all points from = 3.2 to = 5.6 is also given by
( ) = ( ) + ( − )( − ) ( − ) ( − )
( ) = ( ) + 0.256( − 3.2)( − 2.7) ( − 1) ( − 4.8)
( 3) ≈ ( 3) = 22 + 8.4( 3 − 3.2) + 2.856( 3 − 3.2)( 3 − 2.7)
−0.528( 3 − 3.2)( 3 − 2.7)( 3 − 1)
= 20.2120
( 3) ≈ ( 3) = 20.2120 + 0.256( 3 − 3.2)( 3 − 2.7) (3 − 1) ( 3 − 4.8)
= 20.2120 + 0.055296
= 20.267296

4.4. Spline I nterpolation

A. L inear Spline I nterpolation
Given ( , ) , ( , ) , ( , ) ,…, ( , ), + 1 data points
Then interpolate the data to linear spline.
We have an ascending order < < < ⋯<
By drawing a figure to show linear spline

( , )
( , )
( , )
( , )
( , )

The general formula becomes

−
( )= + ( − ) , ≤ ≤
−
Thus, the first spline is,
−
( )= + ( − ) , ≤ ≤
−
The second,
−
( )= + ( − ) , ≤ ≤
−
The term linear spline is,
( )= + ( − ) , ≤ ≤ ------------------------------ (4.26)

Page 56

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample: The upward velocity of a rocket is given as a function of time in the table
below.
Find the velocity at = 16 seconds.
, , / , , /
0 0 0 0
10 227.4 rewrite in to 10 227.4
20 517.35 15 362.78
15 362.78 ascending 20 517.35
22.5 602.97 22.5 602.97
30 901.67 30 901.67

Solution:
−
( )= + ( − ) , ≤ ≤
−
. .
= 362.78 + ( − 15) , 15 ≤ ≤ 20
= 362.78 + 30.913( − 15) , 15 ≤ ≤ 20
( 16) ≈ 362.78 + 30.913( 16 − 15)
= 393.7 ⁄

B. Cubic Spline I nterpolation

The most common piecewise polynomial approximation uses cubic polynomials between
each successive pair of nodes is called cubic spline interpolation.
( )

…
N ow we have
 + 1 grid points , ( = 0,1,…, )
 intervals ≤ ≤ , ( = 0,1,…, − 1)
 cubic spline ( ) , ( = 0,1,…, − 1)
 − 1 interior grid points , ( = 0,1,…, − 1)
 4 coefficient to be determined

Page 57

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Given a function defined on [ , ] and a set of nodes = < < ⋯< = a

cubic spline interpolation for is a function that satisfies the following conditions.

a. ( ) is a cubic polynomial denoted ( ) on the subinterval [ , ] for each

= 0,1,…, − 1
b. ( ) = ( ) for each = 0,1,…,
c. ( )= ( ) for each = 0,1,…, − 2
d. ( )= ( ) for each = 0,1,…, − 2
e. ( )= ( ) for each = 0,1,…, − 2
f. One of the following boundary condition is satisfied
i. ( )= ( ) = 0 (free or natural boundary)
ii. ( )= ( ) and ( ) = ( ) (clamped boundary)

Let a cubic polynomial for the interval is

( )= ( − ) + ( − ) + ( − )+ , = 0,1,…, − 1 ------------ (4.27)
Since this polynomial is valid for both the points and
Therefore, ( )= ( − ) + ( − ) + ( − )+
This implies ( ) = --------------------------------------------------------------------------- (4.28)
( )= ( − ) + ( − ) + ( − )+
( )= ℎ + ℎ + ℎ + ------------------------------------------------------------ (4.29)
Where ℎ = −
N ow the first and second differentiable of equation (4.27), we get
( ) = 3 ( − ) + 2 ( − ) + ----------------------------------------------------- (4.30)
( ) = 6 ( − ) + 2 ---------------------------------------------------------------------- (4.31)
N ow let = ( ) then from equation (4.31) we get at =
= 2 ⇒ = ------------------------------------------------------------------------------- (4.32)
at = , = ( )
= 6 ( − )+ 2
= 6 ℎ +
⇒ = ----------------------------------------------------------------------------------- (4.33)

N ow substituting the values of , and in (4.29)

1
( )= ( − )ℎ + ℎ + ℎ + ( )
6ℎ 2
ℎ
⇒ ( )− ( )= ( +2 )+ ℎ
6
( ) ( )
⇒ = − ( + 2 ) ------------------------------------------------------ (4.34)
Since the curve has equal slope at the point, hence from (4.30)
( )= 3 ( − ) + 2 ( − )+

Page 58

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

⇒ ( )= --------------------------------------------------------------------------------------- (4.35)
( ) ( )
= ( )= − ( +2 ) ----------------------------------------------- (4.36)
But ( ) the slope for the last subinterval is
( )= 3 ℎ +2 ℎ + ---------------------------------------------------- (4.37)
After using , and
( )= 3 ( − )ℎ +2 ℎ .
( ) ( )
+ − ( +2 )
( ) ( )
⇒ ( )= + ℎ -------------------------------------------- (4.38)

From equation (4.35) and (4.36) = 3 ℎ +2 ℎ +

By substituting the values of , , and
( ) ( ) ( ) ( )
− ( + 2 )= + (2 + )
for = 1,2, …,
( ) ( ) ( ) ( )
⇒ + + = − ; -------------------- (4.39)
Where, ( )=
), ( ) = ( ) = ( ( ) , and ( )= ( )
N ow for equally spaced arguments ℎ = ℎ equation (4.39) becomes
+4 + = [ ( ) −2 ( ) + ( )] ; = 1,2,…, − 1---------------- (4.40)

While, ( ) for the equally spaced becomes,

( )= [( − ) + ( − ) ]

+ ( − ) ( )−

+ ( − ) ( )− ; = 0,1,…, − 1----------------------------------- (4.41)

Equation (4.41) gives the cubic spline interpolation while equation (4.40) gives the condition
for .
Remark: ( )= ( ) = 0 (Free or N atural spline) and this one is the most common
boundary condition.
E xample: Find the natural spline interpolating for the given data
0 1 2 3

( ) 1 2 33 244

with ( 0) = ( 3) hence estimate ( 2.5)

Page 59

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

+ 4 + = [ 33 − 2(2) + 1]
Solution: We have ℎ = 1,
+ 4 + = [ 244 − 2(33) + 2]
4 + = 180
⇒
+4 = 1056
Solving this system and then = −24 and = 276
1
( )= [( − ) + ( − ) ]
6
+( − ) ( )−

+( − ) ( )−
1
( )= [( 1 − ) (0) + ( − 0) ( −24) ]
6
+ ( 1 − ) 1 − (0)
+ ( − 0) 2 − ( −24)
= 4 + 5 + 1, 0≤ ≤1
1
( ) = [( − ) + ( − ) ]
6
+( − ) ( ) −
+( − ) ( )−
1
( )= [( 2 − ) ( −24) + ( − 1) (276) ]
6
+ ( 2 − ) 2 − ( −24)
+ ( − 1) 33 − ( 276)
= 50 − 162 + 167 − 53, 1≤ ≤2
1
( ) = [( − ) + ( − ) ]
6
+( − ) ( )−

+( − ) ( )−
1
( )= [( 3 − ) (276) + ( − 2) (0) ]
6
+ ( 3 − ) 33 − (276)
+ ( − 2) 244 − ( 0)
= −46 + 414 − 985 + 715, 2≤ ≤3
( 2.5) = −46(2.5) + 414(2.5) − 985( 2.5) + 715
( 2.5) = . ≈ (2.5)

Page 60

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xercise
∆
1. Evaluate
2. Show that = +∆ + ∆ + ∆
3. I f ( ) = 2 − + 3 + 1 then show that ∆ ( ) = 12 + 10
4. Find the missing figures in the following table
0 5 10 15 20 25
7 11 18 32

5. The following are data from the steam table

Temperature in ° 140 150 160 170 180
Pressure in ⁄ 3.685 4.84 6.302 8.076 10.225
Using the N ewton’s formula, find the pressure of the steam for a temperature of 142° .
6. Given the following score distribution of statistics
30 − 40 40 − 50 50 − 60 60 − 70
. 52 36 21 14
a. The number of students who scored below 35
b. The number of students who scored above 65
c. The number of students who scored between 35 − 45
7. Find a cubic polynomial which takes the following values
1 0 2 3
( ) 2 1 1 10
And also evaluate (4) .
8. A rod is rotating in a plane. The following table gives the angle (in radians) through
which the rod has turned for various values of time seconds.
0 0.2 0.4 0.6 0.8 1.0 1.2
0 0.12 0.49 1.12 2.02 3.20 4.67

9. Use the following table to evaluate 16° and compare with the exact value.
° 0° 5° 10° 15° 20° 25° 30°
° 0 0.0875 0.1763 0.2679 0.3640 0.4663 0.5774

10. Construct the natural cubic spline for the following data.
( )
0.1 −0.62049958
0.2 −0.28398668
0.3 0.00660095
0.4 0.24842440

11. Use the divided-difference method to obtain a polynomial of least degree that fits the
values shown.
a. 0 1 2 −1 3 b. 1 3 −2 4 5
( ) −1 −1 −1 −7 5 ( ) 2 6 −1 −4 2

Page 61

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

12. The following table gives the normal weights of babies during the first 12 months of life:
Age in months 0 2 5 8 10 12
Weight in lbs. 7.5 10.25 16 18 18 21
Find the weight of babies during 6.1 month of life.
13. A natural cubic spline is defined by
( ) = 1+ ( − 1) − ( − 1) 1≤ ≤2
( )= 3
( ) = 1 + ( − 2) + ( − 2) + ( − 2) 2≤ ≤3
4
If is interpolates the data(1,1),(2,1) , and (3,0) , find , , , and .

Page 62

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

CH APT E R FI VE

L E AST SQ U ARE M E T H O D

5.1. L inear L east Square M ethod

I n experimental, social, and behavioral sciences, an experiment or survey often produces
a mass of data. To interpret the data, the investigator may resort to graphical methods.
For instance, an experiment in physics might produce a numerical table of the form
⋯
-------------------------- (5.1)
⋯
And from it, ( + 1) points on a graph could be plotted.
A reasonable tentative conclusion is that the underlying function is linear and that the
failure of the points to fall precisely on a straight line is due to experimental error.
Assuming that = + --------------------------------------------------------------------- (5.2)
What are the coefficients and ? Thinking geometrically, we ask: W hat line most nearly
passes through the nine points plotted?

= +

To answer this question, suppose that a guess is made about the correct values of
and . This is equivalent to deciding on a specific line to represent the data. I n general,
the data points will not fall on the line = + . I f by chance the datum falls on
the line, then
+ − = 0
I f it does not, then there is a discrepancy or error of magnitude
| + − |≠0
The total absolute error for all ( + 1) points is therefore
( , )= ∑ | + − |
I n practice, it is common to minimize a different error function of and :
( , )= ∑ ( + − ) ------------------------------------------------------------- (5.3)

Page 63

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Let us try to make ( , ) a minimum. By calculus, the conditions become

= 0 = 0
(Partial derivatives of with respect to and , respectively) are necessary at the
minimum.
Taking derivatives in Equation (5.3), we obtain
∑ 2( + − ) = 0
∑ 2( + − )= 0
This is a pair of simultaneous linear equations in the unknowns and . They are called
the normal equations and can be written as
(∑ ) + (∑ ) = ∑
-------------------------------------------- (5.4)
(∑ ) + (∑ 1) = ∑
H ere, of course, ∑ 1= + 1, which is the number of data points.
We solve this pair of equations using Cramer’s rule (possible to use other direct methods).
Therefore, we can find the values of and to linear curve fit of equation (5.2).

E xample1: Find the linear least-squares solution for the following table for the values
4 7 11 13 17
2 0 2 6 7

Plot the original data points and the line using a finer set of grid points.
Solution: To find the values of the normal equations using table
N umber of data
0 4 2 8 16
1 7 0 0 49
2 11 2 22 121
3 13 6 78 169
4 17 7 119 289
+ 1= 5
= 52 = 17 = 227 = 644

644 + 52 = 227
52 + 5 = 17
Solving the system of equations, we get = 0.4864 and = 1.6589
Therefore, the linear line that fit the given data points is = 0.4864 − 1.6589
The sum of the total error squares will be obtained as
( , ) = ∑ ( 0.4864 − 1.6589 − )
= 2.9354 + 3.0482 + 2.8612 + 1.7841 + 0.1522 = 10.7813

Page 64

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample2: Find the best fit values of and so that = + fits the data given in
the table
0 1 2 3 4
1 1.8 3.3 4.5 6.3

Solution: (Exercise) Answer: = 0.72 + 1.33

5.2. N on-linear L east Square M ethod

a. Polynomial least square method
The general of approximating a set of data, ( , ) , ,…, with an algebraic
polynomial ( )= + + ⋯+ + of degree
< − 1, using the least square procedure in handled similarly. We choose the
constants , , …, to minimize the least squares error ( , , …, ) , where
= ∑ ( − ( )) = ∑ −2∑ ( ) + ∑ ( ( ))
= ∑ − 2∑ ∑ +∑ ∑
= ∑ − 2∑ ∑ + ∑ ∑ ∑
As in the linear case, for to minimized is necessarily that ⁄ = 0 , ,…,
Thus, for each , we must have
0= = −2 ∑ + 2∑ ∑
This gives ( + 1) normal equations in the ( + 1) unknowns . These are
∑ ∑ = ∑ --------------------------------------------- (5.5)
( , ,…, )
I t is helpful to write equation (5.5) as follows:
∑ + ∑ + ∑ + ⋯+ ∑ = ∑
⎧
∑ + ∑ + ∑ + ⋯+ ∑ = ∑
⎨ ⋮
⎩ ∑ + ∑ + ∑ + ⋯+ ∑ = ∑
These normal equations have a unique solution provided that the ’s are distinct.

E xample: Fit the data below with the discrete least squares polynomial of degree at
most 2.

1 0 1.0000
2 0.25 1.2840
3 0.50 1.6487
4 0.75 2.1170
5 1.00 2.7183

Page 65

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Solution: For this problem, = 2, = 5, and the three normal equations are

1 0 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000

2 0.25 1.2840 0.0625 0.0156 0.0039 0.3210 0.0802
3 0.50 1.6487 0.2500 0.1250 0.0625 0.8244 0.4122
4 0.75 2.1170 0.5625 0.4219 0.3164 1.5878 1.1908
5 1.00 2.7183 1.0000 1.0000 1.0000 2.7183 2.7183

= 5 = 2.5 = 8.7680 = 1.8750 = 1.5625 = 1.3828 = 5.4515 = 4.4015

∑ + ∑ + ∑ = ∑
∑ + ∑ + ∑ = ∑
∑ + ∑ + ∑ = ∑
5 + 2.5 + 1.875 = 8.7680
2.5 + 1.875 + 1.5625 = 5.4515
1.875 + 1.5625 + 1.3828 = 4.4015
Solving the system of equations gives
{ = 0.9268, = 0.6401}
= 1.1383,
Thus the least squares polynomial of degree 2 fitting the data in the above table is
( ) = 0.9268 + 1.1383 + 0.6401
The total error use this table
1 2 3 4 5
0 0.25 0.50 0.75 1.00
1.0000 1.2840 1.6487 2.1170 2.7183
( ) 0.9268 1.2514 1.6560 2.1406 2.7052
− ( ) 0.0732 0.0326 −0.0073 −0.0236 0.0131

= ∑ − ( )) = 0.00720286 is the least that can be obtained by using a

(
polynomial degree at most 2.

b. E xponential
Suppose an exponential curve of the form =
Taking logarithm on both the sides,
We get log = log + log --------------------------------------------------------- (5.6)
That is, = + where = log , = log , = log
The normal equation for (5.6) is,
∑ = + ∑
---------------------------------------------------- (5.7)
∑ = ∑ + ∑
On solving the above two equations, we get and .

Page 66

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Then = log , =
E xample: Fit a curve = to the following data:
x 2 3 4 5 6
y 144 172.8 207.4 248.8 298.5

Solution: Given equation = ⇒ log = log = log + log =

log + log
reduced to = + , where = log , = log , = log
The normal equations are:
∑ log = log + log ∑
∑ log = log ∑ + log ∑
log log
2 144 4 2.1584 4.3168
3 172.8 9 2.2375 6.7125
4 207.4 16 2.3168 9.2672
5 248.8 25 2.3959 11.9795
6 298.5 36 2.3749 14.8494
= 20 = 90 log = 11.5835 log = 47.0254
Putting the values in the normal equations, we have
11.5835 = 5 log + 20 log
47.0254 = 20 log + 90 log
Solving the system of equations log = 2.0401,log = 0.0691 and taking antilog
we have = 109.6731, = 1.1725.
Therefore, equation of the curve is = 109.6731 ∗ (1.1725) .
c. Fitting of the Curve = + √
The error estimate for the point ( , ) is = −( + )
We have, = ∑
= ∑ −( + )
By the principle of least square, the value of is minimum
Therefore, = 0 and = 0
= 0 ⟹∑ 2 − − (− ) = 0

⇒ ∑ = ∑ + ∑ ---------------------------------------------------- (5.8)

and = 0 ⟹∑ 2 − − (− )= 0

⇒ ∑ = ∑ + ∑ ----------------------------------------------- (5.9)

D ropping the suffix in (5.8) and (5.9) then the normal equations are,

Page 67

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

∑ = ∑ + ∑
√
--------------------------------------------------------------- (5.10)
∑ √ = ∑ + ∑
√
E xample: Use the method of least squares to fit the curve = + √ to the
following table of values:
0.1 0.2 0.4 0.5 1 2
21 11 7 6 5 6

Solution:
√ 1 1 √
√
0.1 21 0.3162 3.1623 100 210 6.6402
0.2 11 0.4472 2.2361 25 55 4.9192
0.4 7 0.6324 1.5811 6.25 17.5 4.4268
0.5 6 0.7071 1.4142 4 12 4.2426
1 5 1 1 1 5 5
2 6 1.4142 0.7071 0.25 3 8.4852
1 1
√ √
√
= 4.2 = 4.5171 = 10.1008 = 136.5 = 302.5 = 33.714

∑ = ∑ + ∑
√ 302.5 = 136.5 + 10.1008
⇒
∑ √ = ∑ + ∑ 33.714 = 10.1008 + 4.2
√
Solving the system of equations we get = 1.9733 and = 3.2815
.
H ence the required curve fit equation is = + 3.2815√ .

Page 68

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

5.3. Continuous L east square M ethod

This approximation is concerns for the approximation of functions. Suppose ∈[ , ]
and that a polynomial ( ) of degree at most is required that will minimize the error.
∫ [ ( )− ( )] -------------------------------------------------------------------------- (5.11)
To determine the least square approximation of polynomial that is polynomial to
minimize this expression.
Let ( )= + + ⋯+ + = ∑ and can be defined as
= ( , , .. , ) = ∫ [ ( ) −∑ ]
The problem is to find real coefficients , , . ., that will minimize . A necessarily
condition for the numbers , ,. ., to minimize is that = 0( , ,…, ) .

Since = ∫ [ ( )] − 2∑ ∫ ( ) + ∫ (∑ )
We have, = −2 ∫ ( ) + 2∑ ∫
H ence to find ( ) , the ( + 1) normal equations
∑ ∫ = ∫ ( ) for each = 0,1,…, --------------------------- (5.12)
We will be solved for ( + 1) unknown . The difficulties in obtaining least squares
polynomial approximation ( + 1) × ( + 1) linear system for unknowns , ,. .,
must be solved and the coefficients in the linear system are of the form

∫ = then the matrix in the linear system is known as H ilbert

matrix .
E xample: Find the least square approximating polynomial of degree 2 for the function
( )= on the interval [0,1] .
Solution: The normal equations for ( )= + + are
∫ 1 + ∫ + ∫ = ∫
∫ + ∫ + ∫ = ∫
∫ + ∫ + ∫ = ∫
Performing the integration yields

⎡1 ⎤ ⎡ ⎤
⎢ ⎥ ⎢ ⎥
⎢ ⎥ = ⎢ ⎥
⎢ ⎥ ⎢ ⎥
⎣ ⎦ ⎣ ⎦
We get = ≈ −0.050465 = − = ≈ 4.12251
Consequently, the least square approximating polynomial of degree 2 for the function
( )= on the interval [0,1] is
( ) = −0.050465 + 4.12251 − 4.12251

Page 69

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

5.4. Weighted L east Square M ethod

a. Weighted L east Square for D istrict data points
I f one has more confidence in some data points than others, one can define a weighting
function to give more priority to those particular data points.
Find the best straight line fit for the data in the previous example.
0 1 2 3 4
1 1.8 3.3 4.5 6.3

N ow however, assume that we have five times the confidence in the accuracy of the
final two data points, relative to the other points.
D efine a square weighting matrix :
1 0 0 0 0
⎛0 1 0 0 0⎞
= ⎜0 0 1 0 0⎟
0 0 0 5 0
⎝0 0 0 0 5⎠
N ow we perform the following operations:
∗ =
∗ ∗ = ∗
( ∗ ) ∗( ∗ ) ∗ = ( ∗ ) ∗ ∗
= [( ∗ ) ∗ ( ∗ )] ∗ ( . ) ∗ ( ∗ )
= [ ∗ ] ∗ ∗ ----------------------------------- (5.13)
0.1365
=
1.5105
Then the best fit estimate is = 0.1365 + 1.5105
b. Weighted L east Square for continuous functions
D efinition: An integrable function ( ) is called a weight function on the interval if
( ) ≥ 0 for allin , but ( ) ≠ 0 on any subinterval of .
The purpose of a weighted function is to assign varying degrees of important to
approximations on a certain portions of the interval.
The function which are continuous on [ , ] and are given explicitly, we have
= ∫ ( )[ ( ) − ( )]
Where ( )= + + ⋯+ + = ∑
⇒ = ∫ ( ) ( )−∑ = minimum
The necessarily condition is to have a minimum value is that = 0, = 0,1,. .,
This gives a system of ( + 1) unknown , ,…, . These equations are called
normal equations. The normal equation become

= ∫ ( ) ( )−∑ = 0( , ,…, ) -------------------------------- (5.14)

Page 70

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample: Find the least square approximating polynomial of the polynomial

( ) = 1+ 2 on the interval [0, ] , with the weight function
i. ( )= 1
ii. ( )=
Solution: For = 2 and ( ) = 1
We have = ∫ [( 1 + 2 ) −( + + )]
The normal equations are
= −2 ∫ [( 1 + 2 ) −( + + )] = 0

= −2 ∫ [( 1 + 2 ) −( + + )] = 0

= −2 ∫ [( 1 + 2 ) −( + + )] = 0
This gives

⎧ + + =
⎪ 2 3
⎪
−
+ + =
⎨ 2 3 4 2
⎪
⎪ 2 −3
⎩3 + + =
4 5 6
Solving the system of equations
= 1.9549, = −0.6079, = 0.0000
( ) = 0.0000 − 0.6079 + 1.9549
For = 2 and ( ) =
We have = ∫ [( 1 + 2 )−( + + )]
The normal equations are
= −2 ∫ [( 1 + 2 ) −( + + )] = 0

= −2 ∫ [( 1 + 2 ) −( + + )] = 0

= −2 ∫ [( 1 + 2 ) −( + + )] = 0
This gives

⎧ + + =
⎪
+ + =
⎨
⎪
⎩ + + =
Solving the system of equations
= 3.4441, = −2.5040, = 0.5029
( ) = 0.5029 − 2.5040 + 3.4441

Page 71

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xercise:
1. Fit a straight line = + to the following data by method of least square
0 1 3 6 8
1 3 2 5 4

2. Fit the second degree parabola to the following data by the least square method
1 2 3 4 5
1090 1220 1390 1625 1915

3. Fit a least square geometric curve = to the following data

1 2 3 4 5
0.5 2 4.5 8 12.5

4. A pressure and volume of gas are related by the equation = ( and are
constants). Fit this equation to for the data given below
0.5 1.0 1.5 2.0 2.5 3.0
1.62 1.0 0.75 0.62 0.52 0.46

5. Find the normal equations that arise from filling by the least squares method, an
equation of the form = sin , to set of points ( 0,0) , , 1 , ( , 3) and
( ,2) solve for and .
6. Find the linear least squares polynomial approximation to ( ) on the indicated
interval if
a. ( )= + 3 + 2, [0,1] ;
b. ( )= , [0,2] ;
c. ( )= ln , [1,3] ;
7. Find the least squares polynomial approximation of degree two to the functions and
intervals in exercise 6.
8. Compute the error for the approximations in exercise 6 and 7.

Page 72

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

CH APT E R SI X
N U M E RI CAL D I FFE RE N T I AT I O N AN D I N T E GRAT I O N
6.1. N umerical D ifferentiation
The method of obtaining the derivatives of a function using a numerical technique is
known as N umerical D ifferentiation. N umerical D ifferentiation is required for these two
situations
1. The function values are known but the function is unknown, such functions are
called tabular functions.
2. The function to be differentiating is complicated and therefore it is difficult to
differentiate. The choice of formula
a. N ewton forward/ backward formula – if the derivative at a point near the
beginning/ end of a set of values respectively.
b. N ewton divided difference/ Lagrange I nterpolation formula – if the values of
are not equally spaced.
6.1.1. D erivative using N ewton’s forward interpolation formula:
Consider N ewton’s forward interpolation formula
( ) ( )( )
= + ∆ + ∆ + ∆ + ⋯------------------- (6.1)
! !
= ----------------------------------------------------------------------------- (6.2)
D ifferentiation (6.1) with respect to we get
= ∆ + ∆ + ∆ + ⋯------------------------------- (6.3)
D ifferentiation (6.2) with respect to we get
= -------------------------------------------------------------------------------- (6.4)

N ow from equation (6.3) and (6.4) =

⇒ = ∆ + ∆ + ∆ + ⋯ --------------------- (6.5)
Expression (6.5) becomes the value of at any either tabulated or not.
The formula (6.5) becomes simple for tabulated values of , in particular
when = and = 0.
Putting = 0 in (6.5) we get
= ∆ − ∆ + ∆ − ∆ + ⋯ ---------------------- (6.6)

D ifferentiation (6.5) with respect to

=

= ∆ + ( − 1) ∆ + ∆ + ⋯ ----------- (6.7)
Putting = 0 in (6.7) we have
= ∆ −∆ + ∆ − ⋯ ------------------------------- (6.8)

Page 73

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

6.1.2. D erivative using N ewton’s backward interpolation formula:

Consider N ewton’s backward interpolation formula
( ) ( )( )
= + ∇ + ∇ + ∇ + ⋯------------------- (6.9)
! !
= ---------------------------------------------------------------------------- (6.10)
D ifferentiation (6.9) with respect to we get
= ∇ + ∇ + ∇ + ⋯------------------------------ (6.11)
D ifferentiation (6.10) with respect to we get
= ------------------------------------------------------------------------------ (6.12)

N ow from equation (6.11) and (6.12) =

⇒ = ∇ + ∇ + ∇ + ⋯ -------------------- (6.13)
Expression (6.13) becomes the value of at any either tabulated or not.
The formula (6.13) becomes simple for tabulated values of , in particular
when = and = 0.
Putting = 0 in (6.13) we get
= = ∇ + ∇ + ∇ + ∇ + ⋯ --------------- (6.14)
D ifferentiation (6.13) with respect to
=

= ∇ + ( + 1) ∇ + ∇ + ⋯ --------------- (6.15)
Putting = 0 in (6.15) we get
= ∇ + ∇ + ∇ + ⋯ --------------------------- (6.16)

E xample: From the table below compute and for = 1, = 1.5

and = 2.
1 2 3 4 5 6
1 8 27 64 125 216

Page 74

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Solution: The difference table is

∆ ∆ ∆ ∆
1 1
7
2 8 12
19 6
3 27 18 0
37 6
4 64 24 0
61 6
5 125 30
91
6 216

We have = 1, ℎ = 1, and = 1 at the beginning of the table

= ∆ − ∆ + ∆ − ∆ +⋯

= 7 − ( 12) + ( 6) = 3

= ∆ −∆ + ∆ −⋯
= [ 12 − 6] = 6
We have = 1, ℎ = 1, and = 1.5 ⇒ = 0.5
= ∆ + ∆ + ∆ +⋯
( . ) ( . ) ( . )
= 7+ (12) + (6)
= [ 7 − 0.25] = 6.75
= ∆ − ∆ + ∆ − ∆ +⋯

= 19 − ( 18) + ( 6) = 12

= ∆ −∆ + ∆ −⋯
= [ 18 − 6] = 12
E xample 2: A rod is rotating in a plane about one of its ends. I f the
following table gives the angle radians through which the rod has turned
for different values of time seconds, find its angular velocity at = 0.7 .
, sec 0 0.2 0.4 0.6 0.8 1

, rad 0 0.12 0.48 1.1 2 3.2

Page 75

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Solution: The difference table is

∇ ∇ ∇ ∇
0 0
0.12
0.2 0.12 0.24
0.36 0.02
0.4 0.48 0.26 0
0.62 0.02
0.6 1.1 0.28 0
0.90 0.02
0.8 2 0.30
1.20
1 3.2

.
= = 1, ℎ = 0.2 and = = 0.7 where = = −1.5
.
From N ewton’s interpolation formula, we have
= ∇ + ∇ + ∇ +⋯
.
( . ) ( . ) ( . )
= 1.2 + ( 0.3) + ( 0.02) + ⋯
.
= [ 1.2 − 0.3 − 0.0008]
.
= 4.496 / This represents the angular velocity.
= [∇ + ( + 1) ∇ ]
.
= [ 0.3 + ( −1.5 + 1)( 0.02)]
( . )
= 7.25 / This represents the angular acceleration.

6.2. N umerical I ntegration

N umerical integration is used to obtain approximate answers for definite integrals
that cannot be solved analytically (has no explicit anti-derivative). N umerical
integration is a process of finding the numerical value of a definite integral
= ∫ ( ) ------------------------------------------------------------------------------ (6.17)
Where, a function = ( ) is not known explicitly. But we give only a set of values
of the function = ( ) corresponding to the same values of . To evaluate the
integral, we fit up a suitable interpolation polynomial to the given set of values of
( ) and then integrate it within the desired limits. H ere we integrate an
approximate interpolation formula instead of = ( ) . When this technique is
applied on a function of single variable, the process is called Q uadrature.

Page 76

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

6.2.1. N ewton cotes quadrature formula

Suppose we are required to evaluate the definite integral = ∫ ( )
firstly we have to approximate ( ) by a polynomial ( ) of suitable degree.
Then we integrate the limits with the interval [ , ] .
That is, ∫ ( ) ≈∫ ( ) and the difference

∫ ( ) −∫ ( ) is called the error of approximation.

Let ( ) take the values ( ) = , ( + ℎ) = , …, ( + ℎ) =
when = , = + ℎ, …, = + ℎ respectively.
To evaluate , we replace ( ) by a suitable interpolation formula ( ) . Let
the interval [ , ] be divided into sub-intervals with the division points
= < + ℎ< ⋯< + ℎ = , where ℎ is the width of each sub-
interval. Approximating ( ) by N ewton’s forward interpolation formula
we can write the integral (6.17) as
= ∫ ( )
( )
= ∫ + ∆ + ∆ + ⋯ -------------------------- (6.18)
! !
Since, = ⇒ = + ℎ ⇒ = ℎ
= , ⇒ = = = 0, ⇒ = 0
and when ( )
= + ℎ, ⇒ = = = , ⇒ =
Expression (6.18) can be written as
+ ∆ + ∆ + ∆
= ∫ ℎ
+ ∆ +⋯
Therefore,
∆ ∆
+ ∆ + − + − +
= ℎ --------- (6.19)
∆
+ − + −3 + ⋯
The equation (6.19) is called General Gauss L egendre Q uadrature
formula, for equidistant ordinates from which we can generate any
N umerical integration formula by assigning suitable positive integral value
to .
a. T rapezoidal rule
Substituting = 1 in the relation (6.19) and neglecting all differences
greater than the first we get
= ∫ ( ) = ℎ + ∆ = [2 + − ]

= [ + ] for the sub-interval [ , + ℎ]

Similarly, for the next sub intervals

Page 77

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

= ∫ ( ) = [ + ]

= ∫ ( ) = [ + ]
⋮
= ∫ ( )
( ) = [ + ]
N ow adding , ,…, we get
+ + ⋯+ =
∫ ( ) +∫ ( ) + ⋯+ ∫ ( )
( )

= [ + ]+ [ + ] + ⋯+ [ + ]
⇒ ∫ ( ) = [( + )+ ( + ) + ⋯+ ( + )]
= [( + ) + 2( + + ⋯+ )] -------------------------- (6.20)
The formula (6.20) is called T rapezoidal rule for numerical integration.
The error committed in this formula is given by
( )
≈− ( )= − ( ) .

b. Simpson’s ⁄ rule
Substituting = 2 in the General quadrature formula given by (6.19) and
neglecting the third and other higher order differences
We get
∆
= ∫ ( ) = ℎ 2 + 2∆ + −2

= ℎ 2 + 2( − )+ ( −2 + )
= [ +4 + ]
Similarly,
= ∫ ( ) = [ + 4 + ]

= ∫ ( ) = [ + 4 + ]
⋮
⁄ = ∫ ( )
( ) = [ + 4 + ]
Adding , , …, ⁄ we get,
= + + ⋯+ ⁄ =
∫ ( ) + ∫ ( ) + ⋯+ ∫ ( )
( )

= [ +4 + ]+ [ + 4 + ]+⋯
+ [ + 4 + ]

Page 78

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( + )
= + 4( + + ⋯+ ) -------------------------------------- (6.21)
+ 2( + + ⋯+ )
The above rule is known as Simpson’s one-third rule. The error
committed in Simpson’s one-third rule is given by
( )
≈− ( )( )= − ( )( ) .

c. Simpson’s ⁄ rule
We assume that within any three consecutive sub-intervals of width ℎ,
the interpolating polynomial ( ) approximating ( ) is of degree 3.
H ence substituting = 3, i.e., the General quadrature formula and
neglecting all the differences above ∆ , we get
= ∫ ( )
∆
= ℎ 3 + ∆ + 9− ∆ + − 27 + 9

3 + 9( − )+ ( −2 + )
= ℎ
+ ( −3 + 3 + )

= [ +3 +3 + ]
Similarly
= ∫ ( ) = [ + 3 +3 + ]
⋮
⁄ = ∫ ( )
( ) = [ + 3 +3 + ]
Adding , , …, ⁄ we get,
= + + ⋯+ ⁄ =
( + )+
= 3( + + + + ⋯+ + ) -------------- (6.22)
+ 2( + + ⋯+ )
 Simpson’s three-eighths rule can be applied when the range [ , ]
is divided into a number of subintervals, which must be a
multiple of 3.
 The error in Simpson’s three-eighths rule is
≈− ( )( ) .

Page 79

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

6.2.2. Gauss quadrature formula

Gaussian quadrature chooses the points for evaluation in an optimal, rather
than equally spaced way.
The nodes , , …, in the interval [ , ] and coefficients , ,…,
are chosen to minimize the expected error obtained in the approximation
∫ ( ) ≈∑ ( ) ---------------------------------------------------- (6.23)
For , , …, ,
,…, gives us 2 parameters to choose.
and
Class of polynomials at most 2 − 1
First for any finite intervals [ , ] transformed to [ −1,1]
⇒ ∫ ( ) = ∫ ( ) -------------------------------------------------- (6.24)
The mechanism to convert the limit of integration is use = +
= (1) +
⇒ ⇒ = and =
= ( −1) +
⇒ = + ⇒ =
⇒ ∫ ( ) = ∫ +
= ∫ + ---------------------------- (6.25)
a. O ne point Gaussian quadrature formula
For = 1, 2 − 1 ⇒ 2( 1) − 1 = 1 (degree one polynomial)
∫ ( ) ≈ ( ), ≤ ≤ , ( )= +

∫ ( + ) = + = ( − )+ - (6.26)
( )= [ + ]= + ( ) ------------------------- (6.27)
From (6.26) and (6.27) = − and =
( )( )
⇒ = ( − ) = = ⇒ =
⇒ ∫ ( ) ≈ ( )= ( − ) ( )
If ∫ ( ) ≈ ( 1 − ( −1) ) ( ) = 2 (0) ------------------------ (6.28)
.
E xample: Find the integral ∫ . 5
. . .
Solution: ∫ . 5 = ( 1.3 − 0.1)
= 1.2 ( 0.7)
( . )
= 1.2 5(0.7) = 1.0357

b. T wo point Gaussian quadrature formula

For = 2, 2( 2) − 1 = 3 (degree three polynomial)
∫ ( ) ≈ ( )+ ( ), ≤ ≤ and ≤ ≤

Page 80

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

And also choose the function ( ) = + + +

From integral calculus,
∫ ( + + + )
= (2) + (0) + ( )+ (0) --------------------------- (6.29)
And from

+ + + ( )
( )+ ( )=
+ ( + + + )
( + )+ ( + )
= ----------- (6.30)
+ ( + )+ ( + )
We have the four equations
+ = 2
⎫
+ = 0⎪
and can be determine the four
+ = ⎬
⎪
+ = 0⎭
unknowns , , , .
Solving simultaneously, we get = = 1 and
= ± and = ∓ .
√ √
Therefore, the two point Gauss quadrature is given by
∫ ( ) ≈ − + ---------------------------------------- (6.31)
√ √
.
E xample1: Evaluate the integral = ∫. 5 using Gauss
quadrature two point formula.
Solution: first transfer the interval [0.1,1.3] to [ −1,1] .
Let = +
1.3 = (1) +
⇒ ⇒ = 0.6 and = 0.7
0.1 = ( −1) +
= 0.6 + 0.7
( . . ) . .
Then 5 = 5(0.6 + 0.7) = (3 + 3.5) and
= 0.6
.
Therefore, ∫ . 5 = ∫ (0.6)(3 + 3.5) . .

= ∫ (1.8 + 2.1) . .
.

≈ − +
√ √
. . . .
= (1.8( − ) + 2.1) √ + (1.8( ) + 2.1) √ ,
√ √
= 0.52299 + 0.38719 = 0.91018
.
The exact solution of the integral ∫ . 5 is 0.8939.

Page 81

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

E xample2: Evaluate the integral = ∫ using Gauss quadrature two

point formula.
Solution: first transfer the interval [0,1] to [ −1,1] .
Let =
+
1 = (1) +
⇒ ⇒ = and =
0 = ( −1) +
= +
Then = = and =
( )

Therefore, ∫ = ∫ = ∫
≈ − + = +
√ √
√ √

= 0.41277 + 0.27954 = 0.69231

c. T hree point Gaussian quadrature formula
For = 3, 2( 3) − 1 = 5 (degree five polynomial)
∫ ( ) ≈ ( )+ ( )+ ( ), ≤ ≤ ,
≤ ≤ ≤and ≤
Choose the function ( ) = + + + + +
From integral calculus,
∫ ( + + + + + )
= ( 2) + ( 0) + + ( 0) + ( )+ (0) ------ (6.32)
And from
( )+ ( )+ ( )=
( + + + + + )
+ ( + + + + + )
( + + + + + )
( + + )+ ( + + )
= + ( + + )+ ( + + ) ----- (6.33)
+ ( + + )+ ( + + )
We have the four equations
+ + = 2
+ + = 0⎫
⎪
+ + = ⎪
and can be determine the six
+ + = 0⎬
+ + = ⎪
⎪
+ + = 0⎭
unknowns , , , , , .

Page 82

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Solving simultaneously, we get

= , = , = and = ± , = 0, = ∓ .

Therefore, the two point Gauss quadrature is given by

∫ ( ) ≈ ( )+ ( 0) + (− ) ----------------------- (6.34)

To evaluate = ∫ ≈ 5 ( ) + 8 ( 0) + 5 (− )

= ( 5) + ( 8) + (5) = 0.693122

Exact solution of the integral ∫ is ln 2 = 0.693147

Table of the roots and coefficients of the Gauss-Legendre quadrature

roots coefficients
1 0 2
2 0.57735 1
−0.57735 1
3 0.77460 0.55556
0.00000 0.88889
−0.77460 0.55556
4 0.86114 0.34785
0.33998 0.65214
−0.33998 0.65214
−0.86114 0.34785

E xample: Approximate ∫ cos using Gauss-Legendre quadrature

with = 3
Solution: Taking the entries from the table
0.55556 ∗ . ∗ cos( 0.77460)
∫ cos ≈ + 0.88889 ∗ ∗ cos( 0)
+ 0.55556 ∗ .
∗ cos( −0.77460)

= 0.86151 + 0.88889 + 0.18300 = 0.93340

Using integration by part, the exact solution of ∫ cos is
1.9334214 so, which have small error.

Page 83

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

6.2.3. D ouble integration

D ouble integral as an iterative integral over the region gives
= ∬ ( , ) = ∫ ∫ ( , )
= ∫ ∫ ( , ) ------------------------------- (6.35)
. .
E xample: Evaluate = ∫ ∫
Take
a. ℎ = 0.2 in -direction by Simpson’s 3⁄8 rule, and = 0.3 in
-direction by Simpson’s 1⁄3 rule
b. in both direction by Trapezoidal rule
Solution: Calculate the values of ( , ) =

4 4.2 4.4 4.6

2 0.125 0.1191 0.1131 0.1087

2.3 0.1087 0.1035 0.0988 0.0945
2.6 0.0962 0.0916 0.0874 0.0836
a. Applying Simpson’s 3⁄8 rule in -direction
 For = 2
= ( 0.2)[ 0.125 + 3( 0.1191 + 0.1136) + 0.1087] = 0.0698
 For = 2.3
= ( 0.2)[ 0.1087 + 3( 0.1035 + 0.0988) + 0.0945] = 0.0608
 For = 2.6
3
= ( 0.2)[ 0.0962 + 3( 0.0916 + 0.0874) + 0.0836] = 0.0538
8
Then,
2 2.3 2.6
( ) 0.0698 0.0608 0.0538

Applying Simpson’s 1⁄3 rule in -direction

.
 = [ 0.0698 + 4( 0.0608) + 0.0538] = 0.0367
b. Similarly, applying Trapezoidal rule in -direction
 For = 2
.
= [ 0.125 + 2( 0.1191 + 0.1136) + 0.1087] = 0.0699
 For = 2.3
.
= [ 0.1087 + 3( 0.1035 + 0.0988) + 0.0945] = 0.0608
 For = 2.6
.
= [ 0.0962 + 2( 0.0916 + 0.0874) + 0.0836] = 0.0538

Page 84

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

Then,
2 2.3 2.6
( ) 0.0698 0.0608 0.0538

Applying Trapezoidal rule in -direction

.
 = [ 0.0699 + 2( 0.0608) + 0.0538] = 0.0368

E xercise
1. find and at = 1, = 2

1 3 2 4 5 6
1 27 8 64 125 216
2. Find the first and second derivative for = 15, = 23

15 17 19 21 23 25
= √ 3.873 4.123 4.359 4.583 4.796 5.000
3. The table given below reveals the velocity ‘ ’ of a body during the time ‘ ’ specified. Find its
acceleration at = 1.1.
1.0 1.1 1.2 1.3 1.4
43.1 47.7 52.1 56.4 60.8
4. A rod is rotating a place the following table given the angle (radian) through with the rod has
turned for various values of the time seconds find
a. the angular velocity of the rod
b. its angular acceleration when = 1.0
0 0.2 0.4 0.6 0.8 1.0 1.2
0 0.122 0.493 1.123 2.022 3.200 4.666
5. Evaluate using Trapezoidal rule
i. ∫ ( )
⁄
ii. ∫ (
√
)
6. Evaluate the integral using Simpson’s rule dividing the range in to eight equal parts

log(1 + )
1+
7. Evaluate the problem #6 using two point and three point Gauss quadrature formula.
. .
8. Find the value of the double integral = ∫ ∫ (1 − )
a. using Simpson’s rule ℎ = 0.4, = 0.5
b. using Simpson 3⁄8 rule

9. Evaluate the double integral ∫ ∫ using

( )
a. Trapezoidal rule with ℎ = = 0.25
b. Simpson’s rule ℎ = = 0.25

Page 85

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

CH APT E R SE VE N

N U M E RI CAL SO L U T I O N S O F D I FFE RE N T I AL E Q U AT I O N S

1. T aylor’s Series M ethod

Consider the first order differential equation
= ( , ) with initial value ( )= ---------------------------------------------------------------- (7.1)
D ifferentiating (7.1) with respect to , we get
= + i.e., = + ----------------------------------------------------------------------- (7.2)
( ) ( )
D ifferentiating successively we obtain , , ,…
( )
Putting = and = we get , , , ,…
The Taylor series expansion of ( ) about = is given by
( − ) ( − )
( )= ( )+( − ) ( )+ ( )+ ( )+ ⋯
2! 3!
( ) ( )
= + ( − ) + + + ⋯ --------------------------------------------- (7.3)
! !
Substituting the values of , , , … we obtain ( ) for all values of for which (7.3)
converges.
Let = + ℎ and let ( + ℎ) = ( )= = + + + +⋯
! ! !
Once is known we compute , , ,… from (7.1) and (7.2) etc.
Then can expand in Taylor series about = and we have
( + ℎ) = ( )= = + + + + ⋯ ---------------------------------------- (7.4)
! ! !
Continuing in this way we find the solution of ( ) .
Equation (7.4) can also written as = + + + + 0( ℎ )
! ! !
Where 0( ℎ ) represent all the terms involving fourth and higher power of ℎ.
E xample: Using Taylor’s series method solve = 1+ with = 2
Find a) ( 0.1) b) ( 0.2) c) ( 0.3)
Solution:
a. The Taylor’s algorithm is = + + + +⋯
! ! !
H ere = 0 and ℎ = 0.1
Given = = 1+
Successively differentiating, we get
= +
= 2 +
N ow = ( , ) = 1+ = 1
= ( , )= + = 2
= ( , )= 2 + = 2
. ( . ) ( . )
we get = 2+ ( 1) + ( 2) + ( 2) = 2.1103
! ! !
b. The Taylor’s algorithm for the next approximation is

Page 86

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

= + + + +⋯
! ! !
N ow = ( , ) = 1+ = 1 + ( 0.1) 1103 = 1.21103
= ( , )= + = 2.1103 + ( 0.1) 21103 = 2.2314
= ( , )= 2 + = 2( 2.21103) + ( 0.1) 2.2314 = 2.6452
. ( . ) ( . )
Therefore, = 2.1103 + ( 1.21103) + ( 2.2314) + ( 2.6452) = 2.2430
! ! !
( 0.2) = 2.2430
c. The Taylor’s algorithm for the third approximation is
= + + + +⋯
! ! !
N ow = ( , ) = 1+ = 1.4486
= ( , )= + = 2.53272
= ( , )= 2 + = 3.4037
. ( . ) ( . )
Therefore, = 2.2430 + ( 1.4486) + ( 2.53272) + ( 3.4037) = 2.4011
! ! !
( 0.3) = 2.4011

2. E uler’s M ethod
The Taylor’s series method that we have discussed in the previous section yields solution of
differential equations in the form of a power series.
We now proceed to describe methods which give the solution in the form of table values at equally
spaced points.
Consider the first order differential equation
= ( , ) with initial condition ( )= ---------------------------------------------------------- (7.5)
Suppose that we want to solve (7.5) for at the points = + ℎ( , , ,…)
I ntegrating (7.5) between the points and we get
∫ = ∫ ( , )
H ence = + ∫ ( , ) ---------------------------------------------------------------------------- (7.6)
Assume that ( , ) = ( , ) in ≤ ≤ , we get
= +( − ) ( , )
∴ = + ℎ ( , ) -------------------------------------------------------------------------------------- (7.7)
Similarly if ≤ ≤ , we have

= + ( , )

Substituting ( , ) for ( , ) , we get

= + ℎ ( , ) ------------------------------------------------------------------------------------------ (7.8)
Proceedings like this we obtain the general formula
= +ℎ ( , ) = 0,1,2,…----------------------------------------------------------------- (7.9)
This is called the E uler’s algorithm.

3. M odified E uler’s method

I nstead of approximating ( , ) by ( , ) in equation (7.6),

Page 87

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

We approximate it the Trapezoidal rule by [ ( , )+ ( , )] ,

Which is the mean of the slopes of the tangents at the points corresponding to = and = .
( )
Using this value of we compute = +ℎ ( , ) , then = obtained from Euler’s
method
( ) ( )
Thus, we obtain = + ( , )+ , is the first modified value of .
( ) ( )
= + ( , )+ , is the second modified value of .
( ) ( )
= + ( , )+ , is the third modified value of .
We repeat the process till two consecutive values of agree.
Let be the final value obtained to the desired accuracy. Using this value of we compute
( )
= +ℎ ( + ℎ, ) , then = obtained from Euler’s method
( ) ( )
N ow = + ( , )+ , is the first modified value of .
( ) ( )
= + ( , )+ , is the second modified value of .
( ) ( )
= + ( , )+ , is the third modified value of .
Similarly, we repeat the process till two consecutive values of agree.
Then we proceed to calculate as above and we continue the process till we calculate .
E xample: Using Euler’s method solve = 1+ with ( 0) = 2.
Find ( 0.1) , ( 0.2) and ( 0.3) .
Also find the values by modified Euler’s method.
Solution: Given = ( , ) we know
= + ℎ ( , ) , = 0,1,2,…
H ere ( , ) = 1 + , = 0 , = 2 and ℎ = 0.1
To start putting = 0, we get
( 0.1) = = + ℎ ( , ) = 2 + ( 0.1) ( 0,2) = 2.1
N ow = + ℎ = 0.1
Putting = 1, we get ( 0.2) = = + ℎ ( , ) = 2.1 + ( 0.1)[ 1 + 0.1( 2.1) ] = 2.221
N ow = + ℎ = 0.2
Putting = 2, we get ( 0.3) = = +ℎ ( , )
= 2.221 + ( 0.1) [1 + 0.1( 2.221) ] = 2.3654
H ence ( 0.1) = 2.1 , ( 0.2) = 2.221 and ( 0.3) = 2.3654
( )
Using Modified Euler’s method, starting value = = 2.1 which is obtained by Euler’s method
( ) ( )
= + ( , )+ ,
.
= 2+ [1 + 1 + (0.1)(2.1)] = 2.1105
( ) ( )
= + ( , )+ ,
.
= 2+ [1 + 1 + (0.1)(2.1105)] = 2.1105
Final value = 2.1105
N ow starting value of , using Euler’s method

Page 88

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

( )
= = + ℎ ( + ℎ, )
= 2.1105 + 0.1[1 + ( 0.1)( 2.1105) ] = 2.2316
( ) ( )
= + ( , )+ ,
.
= 2.1105 + [1 + ( 0.1)( 2.1105) + 1 + ( 0.2) ( 2.2316) ] = 2.2434
( ) ( )
= + ( , )+ ,
.
= 2.1105 + [1 + ( 0.1)( 2.1105) + 1 + (0.2) (2.2434)] = 2.2435
( ) ( )
= + ( , )+ ,
.
= 2.1105 + [1 + ( 0.1)( 2.1105) + 1 + (0.2) (2.2435)] = 2.2435
Final value = 2.2435
( )
Starting value of = = + ℎ ( + 2ℎ, ) , using Euler’s method
= 2.2435 + ( 0.1) [1 + ( 0.2) ( 2.2435) ] = 2.3884
( ) ( )
= + ( , )+ ,
.
= 2.2435 + [1 + ( 0.2) ( 2.2435) + 1 + ( 0.3)( 2.3884) ] = 2.4018
( ) ( )
= + ( , )+ ,
.
= 2.2435 + [1 + ( 0.2) ( 2.2435) + 1 + ( 0.3)( 2.4018) ] = 2.4020
( ) ( )
= + ( , )+ ,
.
= 2.2435 + [1 + ( 0.2)( 2.2435) + 1 + ( 0.3) ( 2.4020) ] = 2.4020
Final value = 2.4020
H ence, = 2.1105, = 2.2435 and = 2.4020
4. Runge-Kutta M ethod
The Taylor’s series method of solving differential equation numerical is handicapped by the problem
of finding the higher order derivatives. Euler’s method is less efficient in practical problems since it
requires ℎ to be small for obtaining reasonable accuracy.
The Runge-K utta Methods do not require the calculations of higher order derivatives and they are
designed to give greater accuracy with the advantage of requiring only the function values at some
selected points on the sub-interval.
i. First order R-K method
Euler’s method is the R-K method of first order.
ii. Second order R-K method
The Modified Euler’s formula is
ℎ
= + ( , )+ + ℎ, +ℎ ( , )
2
= + ( + ) where = ℎ ( , ) and = ℎ ( + ℎ, + )
H ence the Modified Euler’s method is the R-K method of second order.
iii. T hird order R-K method

= + ( +4 + ) where = ℎ ( , ), = ℎ + ℎ 2, + 2
= ℎ ( + ℎ, + ,)[ ,
= ℎ ( + ℎ, + )]

Page 89

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

iv. Fourth order R-K method

This method is commonly used and is refereed as Runge-K utta Method.
The working rule for solving the initial value problem
= ( , ) , with initial value ( )= by 4 order R-K method is as follows:
= ℎ ( , )
⎧
= ℎ ( + ℎ⁄2 , + ⁄2)
Calculate successively:
⎨ = ℎ ( + ℎ⁄2 , + ⁄2 )
⎩ = ℎ ( + ℎ, + )
and ∆ = ( + 2 +2 + )
Then the required approximate value is given by = +∆
Similarly the value of in the second interval is obtained by replacing by and by
in the above set of formula and we obtain . I n general, to find substitute
( , ) in the expression of , , , .
E xample: Compute (0.1) and (0.2) by R-K method 4 order for differential
equation. = + , ( 0) = 1
Solution: The formula for the fourth order R-K method is
= ℎ ( , )
⎧
= ℎ ( + ℎ⁄2 , + ⁄2 )
and ∆ = ( +2 +2 + )
⎨ = ℎ ( + ℎ⁄2 , + ⁄2)
⎩ = ℎ ( + ℎ, + )
where ℎ is the interval of differencing and ( , ) is the interval value.
H ere ( , ) = + , = 0, = 1 and ℎ = 0.1
⎧ = ( 0.1)( 0 + 1) = 0.1
⎪ = ( 0.1)[ 0.05( 1.05) + ( 1.05) ] = 0.1155
N ow
⎨ = ( 0.1)[ 0.05( 1.05775) + ( 1.05775) ] = 0.1172
⎪ = ( 0.1)[ 0.1( 1.1172) + ( 1.1172) ] = 0.1360
⎩
∆ = ( 0.1 + 2( 0.1155) + 2( 0.1172) + 0.1360) = 0.1169
∴ = + ∆ = 1 + 0.1169 = 1.1169
∴ ( 0.1) = 1.1169
For the second approximation we have = 0.1 and = 1.1169
⎧ = ℎ ( , ) = 0.1[0.1( 1.1169) + ( 1.1169) ] = 0.1359
⎪ = ℎ ( + ℎ⁄2 , + ⁄2 ) = 0.1[0.15( 1.1849) + ( 1.1849) ] = 0.1582
⎨ = ℎ ( + ℎ⁄2 , + ⁄2) = 0.1[0.15( 1.196) + ( 1.196) ] = 0.1610
⎪
⎩ = ℎ ( + ℎ , + ) = 0.1[0.2( 1.2779) + ( 1.2779) ] = 0.1889
1 1
∴ ∆ = ( + 2 + 2 + ) = ( 0.1359 + 2( 0.1582) + 2( 0.1610) + 0.1889)
6 6
= 0.1605
∴ = + ∆ = 1.1169 + 0.1605 = 1.2774
∴ ( 0.2) = 1.2774

E xercise:

Page 90

Downloaded by Yemane Fsaha ([email protected])

lOMoARcPSD|52342570

Numerical Methods For Engineers

1. Use the Euler’s method to solve the following differential equation over the interval [0,1]
With the step size ℎ = 0.2
= − + + , ( 0) = 1
2. Solve = 1− , ( 0) = 0 using Euler’s method. Find at = 0.1 and = 0.2
Compare the result with result of exact solution.
3. Solve the equation in # 1 using Modified Euler method for the values of at = 0.2, = 0.4 and
= 0.6
4. Use the R-K 4 order method to find numerical solution at = 0.8
For = + , ( 0.4) = 0.41 and take ℎ = 0.2
5. Using R-K method compute the initial value problem = (2 − ) with ( 2) = 1 for =
2.2 = 2.4 (take ℎ = 0.1). Compare the result with result of exact solution.

Page 91

Downloaded by Yemane Fsaha ([email protected])