Numerical Methodsaaa
Numerical Methodsaaa
I am grateful to all those persons in the University of Gondar who have a believed on me for
this work. I am greatly thankful to my department council and those persons who have a
direct impact on these teaching materials.
I sincerely thank my Editors, [Link] A., Yohanis B., Sisay L. for their help and professionalism
on this work.
I would like to thank Mr. Abiyu E. Mr. Gosaye A. and [Link] A. for their support and
guidance on the work.
Finally, I express my heartfelt thanks to my wife Rahel Araya not only for her encouragement
and support but also for sharing all the pain, frustration, and fun of producing the martial.
I would appreciate being informed of errors, or receiving other comments and helpful
suggestions about the book.
The present Teaching material ‘Numerical Analysis I ' is planned to cover the syllabus of
Numerical Analysis I of different Ethiopian Universities in Natural science mainly to the
University of Gondar.
The subject matter is presented in a very systematic and logical manner. In each chapter, all
concepts, definitions and large number of examples in the best possible way have been
discussed in detail and simpler manner so that the students should feel no difficulty to
understand the subject.
In this teaching material, I tried to present algorithm and computer program in met lab, so that
students can see how computation are easily obtained through computer program. Brief
discussion on the contents included in each chapter shall be presented.
Thorough care has been taken to eradicate errors but perfection cannot be claimed. So that the
readers are requested for constructive suggestions to improve the book and it will be gratefully
accepted.
Haileyesus Tessema
June, 2015
haileyesust@[Link]
Table 7.1: Nodes and Weights for the Gauss-Legendre Integration Methods…………248
INTRODUCTION
Numerical technique is widely used by scientists and engineers to solve their problems.
Analytical methods have certain limitations in practical applications. Exact solutions are not
possible. A major advantage for numerical technique is that a numerical answer can be
obtained even when a problem has no analytical solution. However, the 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.
The main purpose of numerical analysis and scientific computing is to develop efficient and
accurate methods to compute approximations to quantities that are difficult or impossible to
obtain by analytic means. It has been convincingly argued that controlling rounding errors is
just a small part of this, and that the main business of computing is the development of
algorithms that converge rapidly. Even if we acknowledge the truth of this statement, it is still
necessary to be able to control different sources of errors, including round off errors, so that
these will not interfere with the computed results.
The errors associated with both calculations and measurements can be characterized with
regard to their accuracy and precision. Accuracy refers to how closely a computed or measured
value agrees with the true value. Accuracy is governed by the errors in the numerical
approximation. Precision refers to how closely individual computed or measured values agree
with each other. Precision is governed by the number of digits being carried in the numerical
calculations. Precision and accuracy are quantified by the errors in a numerical calculation.
These concepts can be illustrated graphically using an analogy from target practice. The
cross sign on each target in Fig.1.l can be thought of as the predictions of a numerical
technique, whereas the bull's-eye represents the truth.
Neither Precise NOR accurate: Precise, but NOT accurate:
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. After the computation is over, the result in the
machine form (with base b) is again converted to decimal form understandable to the users and
some more error may be introduced at this stage. In general, we can say that
Error = True value – Approximate value
Generally the errors may be divided into the following different types:
i) Inherent error: 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. Inherent
6 - 6? =
8@ 8B 8D
A! C! E!
+ + +…
a) Absolute error: 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 J? . Therefore
J? = | − ? |
The unit of exact or unit of approximate values expresses the absolute error.
b) Relative error: The relative error JK , defined by
JK =L L=
8M8N ON
8 8
Where ? is the approximate value of quantity . The relative error is independent of units.
c) Percentage error: The percentage error in ? which is the approximate value of x is given
by
JP = 100XJ? = 100XL L
8M8N
8
EXAMPLE 1.3: Find the absolute, percentage and relative errors if 6 is rounded-off to three
decimal digits. Given 6 = 0.005998.
SOLUTION: If 6 is rounded-off to three decimal places we get 6? = 0.006. Therefore
Error = True value – Approximate value
Error =0.005998 – 0.006 = – 0.000002
Absolute Error = Ea= |JQQRQ| =0.000002
EXAMPLE 1.4: Round-off the number 75462 to four significant digits and then calculate its
absolute error, relative error and percentage error.
SOLUTION: After rounded-off the number to four significant digits we get 75460.
Therefore,
Absolute Error J? = |6 − 6? | =|75462 − 75460| =2
JK =L L=
SMSN ON ;
S S
Relative Error = XCAE; = 0.0000265
L L<10M^_Y
SMSN
S
L L<]x10M^_Y
SMSN
S
L L<10M^
SMSN
S
EXAMPLE 1.5: If 6=8/9 and the exact decimal representation of X is 0.888. Find the absolute
error using an appropriate rule when it is rounded-off to three decimal digits.
SOLUTION: We have 6=8/9, n= 3
The decimal representation of X rounded-off to three decimal digits is X= 0.889
Then
Ea =L − 0.889L = L − L
V V VVU
U U YTTT
=L L
VTTTMVTTY
UTTT
=L L
MY
UTTT
≤ x10M>
]
[ by rule 1 ]
Therefore,
Ea = ≤ x10M>
]
………….+(δb );
dg
:
d8f
: ]+…………………… ……….(1.2)
Errors δY , δ; , … , δb all are small so that the terms containing (δY ); ,(δ; ); ,…,(δb );
and higher powers of δY , δ; , … , δb are being neglected.
Therefore
c6 = δY
dS dS dS
+δ; +…+δb
d8e d8: d8f
………(1.4)
Eq. (1.4) represents the general formula for Errors. If eq. (1.4) divided by X we get relative
error
δie jh δi: jh δik jh
S S d8e S d8: S d8f
δh
Er = = + + ……..+
Also from eq. (1.2), by taking modulus we get maximum absolute error.
EXAMPLE 1.5: If l= and the errors in , o, p be 0.001. Compute the relative maximum
Ai: m=
n@
error in l when = o = p = 1.
SOLUTION: We know cl = c + co + cp
dq dq dq
d8 dr dn
= =
dq Vim= dq Y;i: m: dq MYEi: m=
d8 n@ dr n@ dn nB
Since = , ,
Er =L L = where l = 4 at = o = p = 1
sq T.T>E
q A
=0.009
δl =δt dv =( 12V5-5) δt
dq
Now
δt6100
dq ( Y;wB MC)
q ;wD – Cw
x100 =
x0.056100
( Y;MC)
(;– C)
=
= − >x5 = -11.667%
X
δR =δ Q + δh
d| d|
dK d{
d| K d| MK : Y
dK { d{ ;{: ;
= and = +
~
K
The percentage error in r = X100 = 0.12
⇒ δr =
Y
YTT
(0.12)X 4.5= 0.0054
~
{
The percentage error in h = X100 = 0.505
⇒ δh =
Y
YTT
(0.505)X 5.5= 0.028
Now
= (0.0054)(0.8182)+ (0.028)(0.8347)
= 0.00441828+ 0.0233716
= 0.02778988
Therefore the percentage error in R is
~
|
Ep= X 100
T.T;XXVUVV
(@.B): B.B
_
= X 100
:(B.B) :
T.T;XXVUVV
CT.C/YY
= X 100
= 0.605 %
To find the error in the function X= f(Y , ; , … , b ) is to have a desired accuracy and to
evaluate errors δY , δ; , … , δb in Y , ; , … , b we have
dS dS dS
δX = δY d8 +δ; d8 +…+δb d8
e : f
Similarly, we get
and so on.
This formula is useful where error in the dependent variable is given and we are asked to find
errors in those independent variables of the respective relation.
δh = and δ =
~ ~
; ;
= 0.0025x3
= 0.0075
Now
~{ Y ~
{ { ;
a) Percentage Error in h = x100 = [ ]x100
Y [Link]
Y ; (YC)
= [ ] x100
= [0.00025 x100]
= 0.025
Y [Link]
= Y[ ; (E)
] x100
= [ 0.000625 x100]
= 0.0625
EXAMPLE 1.9: How accurately should the length and time of vibration of a pendulum should
be measured in order that the computed value of g is correct to 0.01%.
SOLUTION: Period of vibration T is given by T = 2/ , where l is the length of pendulum.
⇒
A : d A: d A :
: d
Therefore, g = = : and d =(−2) =
δ = and δ =
~ ~
; ;
……….(1.5)
⇒
~ ~
@:
x100 =0.001 x100 =0.001 ………..(1.6)
:
~ Y ~
;
a. Percentage Error in l = x100 = [ ]x100
Y ~
@:
= [
;( : )
] x100 [ from eq.(1.5) ]
Y ~
@:
= ;[
( : )
x100]
= 0.01= 0.005%
Y
;
[ from eq. (1.6) ]
= x100
Y ~
;
[ from eq.(1) ]
=A @: x100
Y ~
:
Y
A
= x 0.01 = 0.0025% [ from eq. (1.6) ]
The error committed in a series approximation can be evaluated by using the remainder after n
terms. Taylor’s series for () at = is given by,
() = () + ( − )′() + ′′() + ′′′() +…+ (bMY) () +b ()
(8M?): (8M?)= (8M?)fe
;! >! (bMY)!
This term b () is called remainder term and for a convergent series it tends to zero as n→∞.
Thus if we approximate () by the first terms of a series then maximum error committed in
this approximation is given by the b () and if accuracy required is already given then it is
possible to find the number of terms such that the finite series yields the required accuracy.
7 8 =1 + + 7 , 0 < <
8: 8= 8 (fe) 8f
;! >! (bMY)! b!
+ +……….+ +
7 8 =1 + + 7 ,0 < <
8: 8= 8 (fe) 8f
;! >! (bMY)! b!
SOLUTION: Given that + +…+ +
b () = 7
8f
b!
Ea(max) = b! 7 8
8f
Hence at = 1, Er(max) =
Y
b!
Hence we have 12 terms of the expansion in order that its sum is correct to 8 decimal places.
¤ MY = −
8= 8B 8 :fe
-…………+(−1)bMY
> C ;bMY
+ +-------
Find such that series determines ¤ MY (1) correct to eight significant digits.
that their sum gives the value of ¦R§ correct to five decimal places for all values of in
the range − ≤ ≤ . Find also the truncation error.
; ;
2. Two sides and included angle of a triangle are 9.6 cm, 7.8 cm and 45°respectively. Find the
possible error in the area of a triangle if the error in sides is correct to a millimeter and the
angle is measured correct to one degree.
3. The error in the measurement of area of a circle is not allowed to exceed 0.5%. How
accurately the radius should be measured.
4. In a ∆ABC, b = 9.5 cm, c = 8.5 cm and A = 45o, find allowable errors in b, c, and A such that
the area of ∆ABC may be determined nearest to a square centimeter.
5. Round off the following numbers to four significant figures:
a) 8.46235 c) 0.0022218
b) 0.70029 d) 9.235101
6. Round off the following numbers to two decimal places:
i) 48.21416
ii) 2.385
iii) 52.275
iv) 81.255
≈ 0.1
Y
YY
a)
≈ 0.44
A
U
b)
≈ 0.56
C
U
c)
d) 5007932≈ 5008
8. Prove that the relative error of a product of three non-zero numbers does not exceed the sum
of the relative errors of the given numbers.
9. If S = 42o3p–4, find the maximum absolute error and maximum relative errors in S. When
errors in = 1 , o = 2, p = 3 respectively are equal to 0.001, 0.002, 0.003.
10. Compute the relative maximum error in the function u = 7 : , when = o = p = 1 and
8©
8
errors in , o, p be 0.001.
11. Find the number of terms of the exponential series such that their sum yields the value of ex
correct to 8 decimal places at = 1.
INTRODUCTION
One of the most common problem encountered in engineering and science is that given a
function(), find the values of for which() = 0. The solution (values of) are known as the
roots of the equation () = 0,or the zeroes of the function ().
We have seen that expression of the form
() = T + Y + ; ; + > > + … … + b b
where ª 's are constant (b ≠0) and is a positive integer, is called a polynomial in of degree ,
and the equation ( )=0 is called an algebraic equation of degree n.
To solve the equation () = 0, using graphical method we usually write () = 0 as
Y () = ; () where Y () and ; () are both functions of . Then we draw the graphs of oY =
Y () and o; = ; () on the same axes. The real roots of the given equation are the intersections
of the two curves, since at these points oY = o; and so of the function o = ().
EXAMPLE 2.1: Solve the equation R«YT
8
= 1, graphically.
SOLUTION: The given equation R«YT
8
=1
Can be written it as
R«YT
8 Y
8
=
Where R«YT simpler than R«YT , constructing the curves y= R«YT and o = 8, we get
8 Y 8 8 Y
8
and
o=8
Y
1 2 3 4
y= R«YT
8
Direct methods give the exact value of the roots in a finite number of steps. We assume here
that there are no rounds off errors. Direct methods determine all the roots at the same time. For
example, the roots of the quadratic equation a 2+ + ¦ = 0, a ≠0, can be obtained using the
method.
= ;?[-b±√; − 4¦ ]
Y
There are direct methods for finding all the roots of cubic and fourth degree polynomials.
However, these methods are difficult to use. For higher order polynomial and non-polynomial
equations, it is difficult and in many cases impossible to get a closed form solution.
These methods are also called 'trial and error'. Indirect or iterative methods are based on the
concept of successive approximations. We start with one or two initial approximations to the root
and obtain a sequence of approximationsT , Y ,…, b until we obtain the desired degree of
accuracy. These methods generally give only one root at a time. The iterative methods are further
divided into two categories:
i. bracketing and
ii. open methods.
The bracketing methods require the limits between which the root lies, whereas the open
methods require the initial estimation of the solution. Bisection and False position methods are
two known examples of the bracketing methods. Among the open methods, the Newton-
Raphson and the method of fixed point iteration are most commonly used. The most popular
method for solving a non-linear equation is the Newton-Raphson method and this method has a
high rate of convergence to a solution.
In this chapter, we present the following indirect or iterative methods:
1. Bisection Method.
2. Method of False Position (Regular Falsi Method).
3. Secant method.
4. Iteration method (fixed point iteration).
5. Newton-Raphson Method (Newton’s method).
Root 3
Figure: 2.2 Shows the existence of a root of a continuous function b/n , .
2.3.1 Bisection method
It is also called as Bolzano method or Interval halving method. This is one of the simplest iterative
methods based on the intermediate value theorem '' If () is continuous on some interval ,
and () and () are of opposite signs (i.e. ()() < 0), then the equation () = 0 has at
least one real root or an odd number of real roots in the interval (, )''. For definiteness, let
() be negative and () be positive. Then there is a root of () = 0, lying between and .
Let the first approximation be
Y = ( + )
Y
;
(i.e., average of the ends of the range)
; Y
() < 0
; = ( + Y )
Y
;
If (; ) is negative, then the root lies in between ; and , otherwise it lies between ; and Y
and using the graph the third approximation to the root is
> = (Y + ; )
Y
;
and so on. The bisection is repeated until the interval has been reduced to a small value ε (where
0< ε <1), so that |b_Y − b |< ε
Procedure for Bisection Method to Find the Root of the Equation f(x) = 0
Step 1: Choose two initial guess values (approximation) a and b (say > ) such that
(). () < 0.
Step 2: Evaluate the midpoint Y of and given by Y = ( + ) and also evaluate (Y ).
Y
;
Step 3: If (). (Y ) < 0, then set = Y else set = Y . Then apply the formula of step 2
Step 4: Stop evaluation when the difference of two successive values of obtained from step
2, is numerically less than the prescribed accuracy.
EXAMPLE 2.2: Find the root of the equation 3– – 1 = 0 lying between 1 and 2 by
Bisection method.
and (Y ) = 0.875 >0. Therefore, the root will now lie between 1 and 1.5.
Second iteration: The Second approximation is given by
; = ( + Y )= (1+ 1.5 )= 1.25
Y Y
; ;
and (Y ) =– 0.297 < 0. So, the root will now lie between 1.25 and 1.5.
The rest results of the iteration are shown in Table:
v ª (ª ) Interval
a=1 -1
iteration b=2 5 (1,2)
1 1.5 0.875 (1,1.5)
2 1.25 – 0.297 (1.25,1.5)
3 1.375 0.2246 (1.25,1.375)
4 1.313 – 0.0494 (1.313,1.375)
5 1.344 0.0837 (1.313,1.344)
6 1.329 0.0183 (1.313,1.329)
7 1.321 – 0.0158 (1.321,1.329)
8 1.325
From above iterations, the root of the equation up to three places of decimals is 1.325, which is
of desired accuracy.
EXAMPLE 2.3: Find one root of 7 8 – 3 = 0 correct to two decimal places using the method
of Bisection.
SOLUTION: Let () = 7 8 – 3
(1.5) = –0.01831 and (1.6) =0.1530
Since (1.5) = –0.01831, which is negative and (1.6) = 0.1530, which is positive
Therefore, at least one real root will lie (1.5, 1.6).
and (Y ) = 0.06147 > 0. Therefore, the root will now lie between 1.5 and 1.55.
Second iteration: The Second approximation is given by
; = ;(a+ Y )= ;(1.5+ 1.55 )= 1.525
Y Y
and (; ) = 0.02014 > 0. So, the root will now lie between 1.5 and 1.525.
The rest results of the iteration are shown in Table:
ª (ª ) Interval
=1.5 0.01831
iteration =1.6 0.1530 (1.5,1.6)
1 1.55 0.06147 (1.5,1.55)
2 1.525 0.02014 (1.5,1.525)
3 1.5125 0.00056 (1.5,1.5125)
4 1.50625 –0.00896 (1.5125,1.50635)
5 1.50938 –0.00422 (1.5125,1.50938)
6 1.51094 –0.00184 (1.5125,1.51094)
EXERCISE 2.1
1. Find a root of () = 3– 4– 9 = 0 correct to three decimal places by using bisection
method.
2. Find the positive root of () = 3– 3 + 1.06 = 0, by method of bisection, correct to
three decimal places.
3. Compute one positive root of 2– 3§¬ – 5 = 0, by bisection method, correct to three
significant figures.
o − (T ) = ( − T )
¸(8e )M¸(8¹ )
8e M8¹
..……….(2.1)
The method consists in replacing the curve PQ by means of the line PQ and taking the point of
intersection of the line with x-axis as an approximation to next the root.
=¼ ½
8¹ ¸(8e )M8e ¸(8¹ )
¸(8e )M¸(8¹ )
x; = xT − ¼¾(i ½ f(xT )
ie Mi¹
e )M¾(i¹ )
=1 − ¼ ½ (−2)
;MY
;M(M;)
= 1.5
Now, (1.5) = (1.5)3– (1.5)2– 2
= 3.375 – 4.25 = – 0.875 < 0
Thus, the root lies between 1.5 and 2.
Second iteration: let xT = 1.5, xY = 2, (T ) = – 0.875, (Y ) = [Link]
x> = xT − ¼ ½ f(xT )
ie Mi¹
¾(ie )M¾(i¹ )
xA = xT − ¼ ½ f(xT )
ie Mi¹
¾(ie )M¾(i¹ )
=1.6522 − ¼;—T.;YUXÃ(−0.2197)
;[Link];;
= 1.6866
Now, (1.6866) = (1.6866)3– (1.6866)2– 2
=4.7977 – 4.8446 = – 0.0469 < 0
Thus, the root lies between 1.6866 and 2.
Fourth iteration: let xT = 1.6866, xY = 2, (T ) = – 0.0469, (Y ) = [Link]
= 1.6938
Now, (1.6938) = (1.6938)3– (1.6938)2– 2
= 4.8594 – 4.8690 = – 0.0096 < 0
Thus, the root lies between 1.6938 and 2.
Fifth iteration: let xT = 1.6938, xY = 2, (T ) = – 0.0096, (Y ) = [Link]
x> = xT − ¼ ½ f(xT )
ie Mi¹
¾(ie )M¾(i¹ )
= 1.6938 − º » (– 0.0096 )
;[Link]>V
;MÀ– [Link] Á
= 1.6953
Now, (1.6953) = (1.6953)3– (1.6953)2– 2
= 4.8724 – 4.8740 = – 0.0016 < 0
Therefore, the root lies between 1.6953 and 2.
Sixth iteration: let xT = 1.6953, xY = 2, (T ) = – 0.0016, (Y ) = [Link]
= 1.6953 + 0.0002
= 1.6955
Hence, the root is 1.695 correct to three places of decimal.
EXAMPLE 2.5: Use the Regula-Falsi methods to determine the root of the equation cos − 7 8 =0
which lies in the interval (0, 1).
SOLUTION: Let () = ¦R§ − 7 8 = 0
(0) = 1and (1) = −2.17797952
Since (0) and (1) are of opposite signs, therefore the root lies between 0 and 1, so taking
xT =0, xY =1, (T ) = 1 and (Y ) = −2.17797952, then by Regula-Falsi method, we get
First iteration:
x; = x T − ¼ ½ f(xT )
ie Mi¹
¾(ie )M¾(i¹ )
= xT − ¼ ½ (1)
YMT
M;.YXXUXUC;MY
= 0.3146653378
Now, f (0.3146653378) = cos (0.3146653378) – (0.3146653378)7 (T.>YAEEC>>XV)
=0.519871175 > 0
So that the root lies between 0.3146653378 and 1.
Second iteration: let xT =0.3146653378, xY =1, (T ) = 0.519871175 and (Y ) = -2.17797952
x> = xT − ¼ ½ f(xT ),
ie Mi¹
¾(ie )M¾(i¹ )
= 0.4467281466
Now, (0.4467281466) = ¦R§ (0.4467281466)– (0.4467281466)7 ([Link];VYAEE)
= 0.203544710 > 0
So that the root lies on the interval (0.4467281466, 1)
The rest results of the iteration are shown in Table:
=2−¼ ½ (−9)
>M;
YM(MU)
=2.9
Now, (2.9) =( 2.9)3– 9(2.9 ) + 1
=– 0.711 < 0
So that the root lies between 2.9and 3.
Second iteration: let T = 2.9, Y = 3, (T ) = −0.711, (Y ) = [Link]
= 2.9416
Now, (2.9416) = ( 2.9416)3– 9(2.9416) + 1
=– 0.0207< 0
So that the root lies between 2.9416 and 3.
The rest results of the iteration are shown in Table:
ª (ª ) Interval
=2 -9
iteration =3 1 (2,3)
1 2.9 – 0.711 (2.9,3)
2 2.9416 – 0.0207 (2.9416,3)
3 2.9428 – 0.0003 (2.9428,3)
4 2.942817
EXERCISE 2.2
1. Compute the root of the equation ; + 4 + 4 = 0, by regular-Falsi method, correct to
three decimal places.
2. Compute the root of the equation > − 4 − 9 = 0, by Regular –Falsi method, correct
to two decimal places.
3. Compute the root of the equation §¬ + ¦R§ = 1, by Regula-Falsi method, correct to
four decimal places.
Geometrically, in Secant method we replace the function () by a straight line passing through
the points (bMY , bMY ) and (b , b ) and take the point of intersection of the straight line
with the x-axis as the next approximation to the root.
; ! Y I ¼ ½ Y
8e M8¹
¸8e M¸8¹
=1I¼ ½ I3
YMT
M>MY
= 0.25
Hence, ; ! 0.25and ; ! – 0.23437.
Second iteration:
> ! ; I ¼ ½ ;
8: M8e
¸8: M¸8e
A = > − ¼ ½ (> )
8= M8:
¸(8= )M¸(8: )
= 0.186441 − º » (0.074276)
[Link].;C
[Link];XEMÀ– T.;>A>XÁ
=0.201736
Hence, A = 0.201736 and (A ) = – 0.000470.
Fourth iteration:
C = A − ¼ ½ (A )
8@ M8=
¸(8@ )M¸(8= )
= 0.201640
EXAMPLE 2.8: Find the root of the equation () = 4 §¬ + 2= 0 by Secant method.
SOLUTION: In this method we neglect the condition (b ). (bMY ) < 0. Initially take T = −1
and Y = −2 and (T ) = – 2.36588 ,(Y ) = 0.36281, the next approximation to the root by
Secant method is given by
First iteration:
; = Y − ¼ ½ (Y )
8e M8¹
¸(8e )M¸(8¹ )
= −1 − º » (0.36281)
M;M(MY)
T.>E;VYMÀ– ;.>ECVV Á
= – 2 + 0.13296
=– 1.86704
Hence, ; =– 1.86704 and (; ) = – 0.33992.
Second iteration: Now, the next approximation > can be obtained by using the values of Y and
; in Secant method, which is given by
> = ; − ¼ ½ (; )
8: M8e
¸(8: )M¸(8e )
=– 1.93384
Hence, A =– 1.93384 and (A ) = 0.00045.
Fourth iteration:
= – 1.93384 – º » (0.00045)
(– Y.U>>VA) –(– Y.U>Y>C)
[Link]À– [Link];EUÁ
= – 1.93375
Hence, C =– 1.93375 and (C ) = – 0.00002.
Fifth iteration:
= – 1.93375 – º » (– 0.00002)
À– Y.U>>XCÁ–(– Y.U>>VA)
À– [Link];Á[Link]
= – 1.93375
Thus, the approximation value to the root is –1.93375, correct up to five decimals.
EXAMPLE 2.9: Find the root of the equation 3– 2 – 5 by Secant method.
SOLUTION: Here () = 3– 2 – 5 then (2) = – 1, (3) = 16 and (2) (3) < 0
Therefore initial approximations are taken as T = 2, Y = 3, (T ) =-1 and (Y ) = 16
Then by Secant method, the next approximation is given by
First iteration:
=3−¼ ½ (16)
>M;
YEM(MY)
= 2.0588
Hence, ; = 2.0588and (; ) = – 0.3911.
Second iteration:
= 2.0813
Hence, > = 2.0813and (> ) = – 0.1468.
Third iteration:
A = > − ¼ ½ (> )
8= M8:
¸(8= )M¸(8: )
= 2.0813 − º » (– 0.1468)
;.TVY> M;.TCVV
– [Link]À– T.>UYYÁ
= 2.0948
Hence, A = 2.0948 and (A ) = – 0.0028.
Fourth iteration:
E = C − ¼ ½ (C )
8B M8@
¸(8B )M¸(8@ )
= , = ( 5 − 1)= , =
e
8 = _Y C8MY
C 8
, etc. .............. (2.4)
Now, finding a root of () = 0 is same as finding a number α such that = Æ(), that is, a
fixed point of Æ(). A fixed point of a function Æ is a point α such that = Æ(). This result is
also called the fixed point theorem.
Using Eq. (2.3), the iteration method is written as
b_Y = Æ(b ), = 0, 1, 2, … ……………(2.5)
The function Æ() is called the iteration function. Starting with the initial approximationT , we
compute the next approximations as
Y = Æ(T )
; = Æ(Y )
……………….
b_Y = Æ(b )
The stopping criterion is same as used earlier. Since, there are many ways of writing () = 0 as
= Æ(), it is important to know whether all or at least one of these iteration methods converges.
REMARK: Convergence of an iteration method b_Y = Æ(b ), n= 0, 1, 2, ..., depends on the
choice of the iteration function Æ(), and a suitable initial approximation x0, to the root.
Consider again, the iteration methods given in eq. (2.4), for finding a root of the equation f () =
3– 5 + 1 = 0 The positive root lies in the interval (0, 1).
b =
8f = _Y
C
i. , n= 0, 1, 2, ….
b = ( 5b − 1)= , n= 0, 1, 2, …
e
ii.
With T = 1, we get the sequence of approximations as
Y = 1.5874, ; = 1.9072, > = 2.0437, A = 2.0968, ….
b =
C8f MY
8f
iii.
Æ() = , Æ′() = , and |Æ′()|= < 1 for all in 0 < <1. Hence, the method
8 = _Y >8 : >8 :
C C C
i.
Æ() = ( 5 − 1)= , Æ′() = . Now |Æ′()|< 1, when is close to 1 and| Æ′() | > 1
e
C
:
>( C8MY)=
ii.
| Æ′() | > 1 in the other part of the interval. Convergence is not guaranteed.
Procedure for Iteration Method to Find the Root of the Equation f(x) = 0
Step 1: Take an initial approximation as T .
Step 2: Find the iteration function which satisfies the condition |Æ′()|<1 and the first
approximation Y by using the formula Y = Æ(T )
Step 3: Follow the above procedure to find the successive approximations b_Y using the
formula b_Y = Æ(b )when =1, 2,…..
Step 4: Stop the process when where ε is the prescribed accuracy.
|b_Y − b |< ε
EXAMPLE 2.10: Find a real root of 3– 2– 3 = 0, correct to three decimal places using fixed
point iteration method.
SOLUTION: Here () = 3– 2– 3 = 0
Also (0) = −3 < 0, (1) = −4 < 0, and (2) = 1 > 0. Therefore, root of the equation
lies between 1 and 2. Since (1) < (2), we can take the initial approximation T = 1
now we have to rewrite the given equation in the following forms, such that we will take for the
successive approximation the one which satisfies the condition |Æ′()| < 1 at the initial
approximation.
=
8 = M>
;
i.
= 8 :M;
>
ii.
Case 1: let = Æ()= , then Æ′() = and |Æ′(T )| = |Æ′(1)| =1.5 > 1. So it is failed.
8 = M> >8 :
; ;
Case 2: let = Æ()= , then Æ′() = and |Æ′(T )| = |Æ′(1)| =6 > [Link] failed.
> >(;8)
8 : M; (8 : M;):
Case 3: let = Æ()=(2 + 3)= , then Æ′() = and |Æ′(T )| = |Æ′(1)|=0.24000 <[Link]
e
;
:
>(;8_>)=
Second iteration:
Third iteration:
Fourth iteration:
Fifth iteration:
Sixth iteration:
Seventh iteration:
Eighth iteration:
Hence, from the 7th and 8th iterations, we can conclude that one of the real roots of () = 0 is
1.8933 correct to four decimal places.
EXAMPLE 2.11: Find a real root of ¦R§ = 3 – 1, correct to three decimal places using fixed
point iteration method.
SOLUTION: Here, we have () = ¦R§ – 3 + 1.
= (¦R§ + 1)
Y
>
ii.
Now let's check that which satisfies the condition |Æ′()|<1 at the initial approximation.
|Æ′(T )| = L L =L L |§¬|< 1.
M˪b8 Y
> >
Second iteration:
Third iteration:
Fourth iteration:
Fifth iteration:
Sixth iteration:
Therefore, from the 5 and 6th iterations, we can conclude that the real roots of () = 0 is
th
We obtain
Æ Ç ()=
A(U8 : _;8
>8 = _8 : _Y;
We find | φ′(x) |< 1 for all x in the interval (– 1, 0). Hence, the iteration converges.
The iteration method is written as
b_Y = − >(8
A
f)
= _(8
f)
: _Y;
Let T = – 0.25be the initial approximation to the root. We obtain the following results.
First iteration:
Y = Æ(T )= −
A
>(8¹ )= _(8¹ ): _Y;
= – 0.33290
Second iteration:
; = Æ(Y )= −
A
>(8e )= _(8e ): _Y;
= – 0.33333
Third iteration:
> = Æ(; )= −
A
>(8: )= _(8 :)
: _Y;
= – 0.33333
3. Use the method of Iteration to find a positive root between 0 and 1 of the equation
§¬ =
8_Y
8MY
.
Since ℎ is small, we can neglect second and higher degree terms in hand therefore, we get
(T ) + ℎ′(T ) = 0 …………(2.12)
From which we have,
Hence, if T be the initial approximation, then next (or first) approximation Y is given by
Y = T +ℎ = T −
¸(8¹ )
¸Ç(8¹ )
; = Y − ¸Ç(8e )
¸(8 )
e
In general
b_Y = b −
¸(8f )
¸Ç(8f )
Step 3: Follow the above procedure to find the successive approximations b_Y using the
formula
Step 4: Stop the process when where ε is the prescribed accuracy |b_Y − b |< ε.
b_Y = b − = b −
¸(8f ) 8f
@ – YY8 _ V
f
¸Ç(8f ) = – YY
A8f
8f ÀA8f
= – YYÁM(8 @ – YY8 _ V)
f f
A8f
= – YY
=
A8f
@ – YYi M8 @ _YY8 M V)
k f f
= – YY
A8f
=
>8f
@M V
A8f
= – YY
=
First iteration:
Y =
>8¹@ M V
A8¹= – YY
= 1.90476
Second iteration:
; =
>8e@ M V
A8e= – YY
= 1.89209
Third iteration:
> =
>8:@ M V
A8:= – YY
= 1.89188
Fourth iteration:
A =
>8=@ M V
A8== – YY
= 1.89188
Therefore, the 3th and 4th iteration values are the same, hence the real roots of () = 0 is
1.89188correct to five decimal places.
EXAMPLE 2.14: Find the smallest positive real root of the equation = 7 M8 using the Newton-
Raphson’s method.
SOLUTION: Let () = 7 8 − 1 then ′() = 7 8 +7 8 = (1 + ) 7 8 and (0) = −1 < 0
,(1) =1.718282>0. Then let's take T = 1. Now f(b ) =b 7 8f – 1 and ′(b ) = (1+b )7 8f , then
using the Newton-Raphson formula, we have
b_Y = b − = b -
¸(8f ) ik Ì f – Y
¸Ç(8f ) (Y_ik )¢£k
=
8f ¢£k _(ik): ¢£k Mik Ì f _ Y)
(Y_ik )¢£k
=
(ik ): ¢£k _ Y)
(Y_ik)¢£k
First iteration:
Y =
(i¹ ): ¢£¹ _ Y)
(Y_i¹ )¢£¹
= 0.6839397
Second iteration:
; =
(ie ): ¢£e _ Y)
(Y_ie )¢£e
= 0.5774545
Third iteration:
> =
(i: ): ¢£: _ Y)
(Y_i: )¢£:
= 0.5672297
Fourth iteration:
A =
(i= ): ¢£= _ Y)
(Y_i= )¢£=
= 0.5671433
b_Y = b − = b -
¸(8f ) 8f @ M8f MYT
¸Ç(8f ) A8f = MY
=
8f (A8f = MY))M(8f @ M8f MYT)
A8f = MY
=
A8f @ M8f M8f @ _8f _YT)
A8f = MY
>8f @ _YT)
A8f = MY
=
First iteration:
Y =
>8¹ @ _YT)
A8¹ = MY
= 1.8710
; =
>8e @ _YT)
A8e = MY
= 1.8558
Third iteration:
> =
>8: @ _YT)
A8: = MY
= 1.8556
formula
b_Y = b −
¸(8f )
¸Ç(8f )
e
( M?)
= b –
f
e
M
f :
= b (2−b )
ii. Square root: The quantity √ can be considered as a root of the equation 2− = 0. Here
() = 2− ⇒ ′() = 2 .From this, by Newton's formula, we get a recursion formula
b_Y = b −
¸(8f )
¸Ç(8f )
= b − (b +
8f : M? Y ?
;8f ; 8f
= )
Inverse square root: the inverse square root is the root of the equation8: − = 0. From
Y
iii.
iv. Formula for pth root and reciprocal pthroot: For computing pth root of a, we can solve the
equation P – =0. Here () = P – ′() = Í PMY .
Hence, by Newton's formula, we obtain the recursion formula
b_Y = b −
¸(8f )
¸Ç(8f )
=b − P8f
8 Î M?
Îe
f
(PMY)8f Î _?
P8f Îe
= for p=2, 3, 4 ,…..
b_Y =
;8f = _? Y ?
= (2b + )
>8f : > 8f
Similarly, the reciprocal of Íth root of a can be obtained by solving the equation 8Î - a =0.
Y
Then
= b −
8f Î M?
MP8f (Îe)
= b [ ] for Í=3, 4, 5, ….
(PMY)_?8f Î
P
b_Y = b [
;_?8f =
>
]
EXAMPLE 2.16: Evaluate √29 to five decimal places by Newton-Raphson iterative method.
SOLUTION: Let = √29 , then ; − 29 = 0
We consider () = ; − 29 and ′() = 2
The Newton-Raphson iteration formula gives
b_Y = b − = b -
¸(8f ) 8f : M;U
¸Ç(8f ) ;8f
:
8f ( ;8f )M(8f M;U)
;8f
=
:
= ;( b + 8 )
;8f : M8f _;U) 8f : _;U Y ;U
;8f ;8f
= =
f
Now (5) = – 4 < 0 and (6) = 7 > 0. Hence, a root of () = 0 lies between 5 and 6.
Taking T = 5.3, as the initial approximation and using the above simplified formula, we have
First iteration:
Y = ( T +
Y ;U
; 8¹
) = 5.38585
Second iteration:
; = ;( Y + 8 ) = 5.38516
Y ;U
e
Third iteration:
> = ( ; +
Y ;U
; 8:
) = 5.38516
√8<√17< √27
= = =
[taking cubic root to each side]
2 < √17< 3
=
b_Y =
;8f = _? Y ?
>8f :
= >(2b +8 )
f
Y = (2T +
Y ?
> 8¹
) = 2.5733
Second iteration:
; = (2Y +
Y ?
> 8e
) = 2.5713
Third iteration:
Fourth iteration:
A = (2> +
Y ?
> 8=
) = 2.57138
Since > and A are accurate to four decimal places hence the required root is 2.5713.
EXERCISE 2.6
1. A root of () = 7 8 – 22 lies in the interval (1, 2). Determine this root with the Newton-
Raphson method to five decimal accuracy.
2. Find the cube root of 10.
limª→Ò L ÌÓ¥e
L ≤A
Ì
Ó
Where A, is a non-zero finite number called asymptotic error constant and it depends on
derivative of () at an approximate root 7ª and 7ª_Y are the errors in successive approximation.
Hence 7ª and 7ª_Y are the errors in ¬ th and (¬ + 1)th iterations respectively. Comparing the above
equation with
limª→Ò L L ≤A
ÌÓ¥e
ÌÓ
We get k= 1 and A= 0.5. Thus the Bisection Method is first order convergent or linearly
convergent.
Let α denotes the exact value of the root of () = 0, and let ª andª_Y , be two successive
approximations to the actual root α. If 7ª and 7ª_Y are the corresponding errors, we have
ª = + 7ª and ª_Y = + 7ª_Y
by Newton-Raphson’s iterative formula
+ 7ª_Y = + 7ª −
¸(Ô_ÌÓ )
¸Õ (Ô_ÌÓ )
7ª_Y − 7ª = −
¸(Ô_ÌÓ )
¸Õ (Ô_ÌÓ )
,
7ª_Y = 7ª − ¸Õ (Ô_ÌÓ )
¸(Ô_Ì )
Ó
Ö :
ÌÓ ¸Õ (Ô)_ Ó ¸ÕÕ (Ô)_⋯
= 7ª − Ö
:!
:
¸Õ (Ô)_ÌÓ ¸ÕÕ (Ô)_ Ó ¸ÕÕÕ (Ô)_⋯
[ b/c f(α) = 0 ]
:!
ÌÓ : ¸ÕÕ (Ô)
; ¸Õ (Ô){Y_ Ì ÙÕÕ (Ú)}
=
Ó Õ Ù (Ú)
+ 7ª
ÌÓ : ¸ÕÕ (Ô) ¸ÕÕ (Ô) -1
; ¸Õ (Ô) ¸Õ (Ô)
= [1 ]
− 7ª
ÌÓ : ¸ÕÕ (Ô) ¸ÕÕ (Ô)
; ¸Õ (Ô) ¸Õ (Ô)
= {1 +…..}
− { } +……….
ÌÓ : ¸ÕÕ (Ô) ÌÓ = ¸ÕÕ (Ô) ;
; ¸Õ (Ô) ; ¸Õ (Ô)
=
= − { } +……..
ÌÓ¥e Y ¸ÕÕ (Ô) ÌÓ ¸ÕÕ (Ô) ;
ÌÓ : ; ¸Õ (Ô) ; ¸Õ (Ô)
Where ª_Y , ª , ªMY are successive approximations to the required root of () = 0.
Let α be the actual (true) root of () = 0, i.e., () = 0. If 7ª_Y , 7ª , 7ªMY are the successive
errors in the (¬+1)th, (¬)th, (¬ −1)th iterations respectively, then
ª_Y = + 7ª_Y , ª = + 7ª and ªMY = + 7ªMY ………………(2.14)
Using these in eq. (2.14), we obtain
+ 7ª_Y = + 7ª −
(¢Ü M¢Üe )¾(Ý_¢Ü )
¾(Ý_¢Ü )–¾(Ý_¢Üe )
7ª_Y =7ª −
(ÌÓ MÌÓe )¸(Ô_ÌÓ )
¾(Ý_¢Ü )–¸(Ô_ÌÓe )
…………….(2.15)
Ö :
(ÌÓ MÌÓe ){¸(Ô)_ÌÓ ¸Õ (Ô)_ Ó ¸ÕÕ (Ô)}
= 7ª − :!
Ö :ÖÓe : ÕÕ
à(ÌÓ MÌÓe )¸Õ (Ô)_{ Ó }¸ (Ô)á
, [ on ignoring the higher order terms ]
:!
Ö :
{ ¸(Ô)_ÌÓ ¸Õ (Ô)_ Ó ¸ÕÕ (Ô) }
= 7ª − Ö Ö
:!
â¸Ç(Ô)_{ Ó Óe }¸ÕÕ (Ô)ã
:
Ö :
{ ÌÓ ¸Õ (Ô)_ Ó ¸ÕÕ (Ô) }
7ª − Ö Ö
:!
â¸Ç(Ô)_{ Ó Óe }¸ÕÕ (Ô)ã
=
:
Ö : ÙÕÕ (Ú)
{ ÌÓ _ Ó Õ }
= 7ª −
:! Ù (Ú)
Ö Ö ÙÕÕ (Ú)
àY_{ Ó Óe } Õ á
: Ù (Ú)
= 7ª − { 7ª + }[1 +{ } ¸Õ (Ô)
ÌÓ : ¸ÕÕ (Ô) ÌÓ MÌÓe ¸ÕÕ (Ô) -1
;! ¸Õ (Ô) ;
= 7ª − { 7ª + }[1 − â ã ¸Õ (Ô) + ⋯
ÌÓ : ¸ÕÕ (Ô) ÌÓ MÌÓe ¸ÕÕ (Ô)
;! ¸Õ (Ô) ;
= 7ª − + − ¸Õ (Ô)
ÌÓ (ÌÓ MÌÓe) ¸ÕÕ (Ô) ÌÓ : ¸ÕÕ (Ô) ÌÓ : (ÌÓ MÌÓe) ¸ÕÕ (Ô) ;
; ¸Õ (Ô) ;! ¸Õ (Ô) A
Where M =
¸ÕÕ (Ô)
;¸Õ (Ô)
and would be a constant.
In order to find the order of convergence, it is necessary to find a formula of the type
7ª_Y =å7ª , with an appropriate value of k …………….(2.18)
With the help ofeq. (5), we can write
7ª =å7ªMY or 7ªMY =( Ó )
Ì e
æ
å7ª =7ª ( æÓ ) .M
Ì e
Therefore the order (or rate) of convergence of false position method is 1.618. Hence regula-falsi
is supper linear convergence.
A multiple root corresponds to a point where a function is tangent to the x-axis. For example, a
double root results from
() = ( − 3)( − 1)( − 1)
or, multiplying terms, () = 3−5 2+7 − 3. The equation has a double root because one value
of makes two terms equal to zero. Graphically, this corresponds to the curve touching the x-axis
tangentially at the double root.
A triple root corresponds to the case where one x value makes three terms in an equation equal
to zero, as in
() = ( − 3)( − 1)( − 1)( − 1)
or, multiplying terms, () = 4−6 3+122−10 + 3. Notice that the graphical depiction again
indicates that the function is tangent to the axis at the root, but that for this case the axis is
crossed. In general, odd multiple roots cross the axis, whereas even ones do not. For example,
the quadruple root in does not cross the axis.
Let α be a root of the equation
() = 0 ……………….(2.21)
which is repeated m times, then () =0, ′() = 0, ….., (°MY) () =0 and (°) (α) ≠0. Such
equation can be written as
() = ( − )° «() = 0
MULTIPLE ROOTS BY NEWTON'S METHOD: Let be a root of the polynomial equation () = 0
which is repeated m times. If T , Y , …., b_Y be its successive approximations, then on the lines
of Newton's iterative method , we have
T − ¯ , T − (¯ − 1) , b − (¯ − 2)
¸(8¹ ) ¸Ç(8¹ ) ¸ÕÕ (8f )
¸Ç(8¹ ) ¸ÇÇ(8¹ ) ¸(=) (8¹ )
, ………
EXAMPLE 2.18: Find the double root of the equation 3− 2− + 1 = 0
SOLUTION: Let () = 3− 2− + 1. Then ′() = 32−2 − 1 and ′′() = 6 − [Link]
(1) = 0 and (−1) = [Link] this = 1 è = −1 are the roots of . Since is a
polynomial of degree 3, one of this root is a double root.
Starting with T =0.9, we have
T − 2
¸(8¹ ) [Link]
¸Ç(8¹ ) MT.>X
= 0.9 – 2( ) = 1.003
And
T −
¸Ç(8¹ ) MT.>X
¸ÇÇ(8¹ ) >.A
= 0.9 - = 1.009
The closeness of these values implies that there is a double root near = 1
Therefore choosing Y =1.01 for the next approximation, we get
Y − 2
¸(8e ) [Link];
¸Ç(8e ) [Link]>
= 1.01 – 2( ) = 1.001
And
Y −
¸Ç(8e ) [Link]>
¸ÇÇ(8e ) [Link]
= 1.01 – = 1.001
This shows that there is a double root at =1.001, which is quite near the actual root = 1.
Here o = Í() is approximated by a second degree parabola passing through these three points
Figure 2.6: Muller’s method
To find the values of A and B, we will solve simultaneously. To find value of A multiply the top of
eq. (2.25) by (ªM; − ª ) and the bottom by (ªMY − ª ). Then we obtain
(8Ó: M8Ó )(rÓe MrÓ )M(8Óe M8Ó )(rÓ: MrÓ )
(8Óe M8Ó: )(8Óe M8Ó )(8Ó: M8Ó )
A=
Similarly to find the value of B, multiply the top of eq. (2.25) by (ªM; − ª ); and the bottom by
(ªMY − ª ); . Then we obtain
(8Ó: M8Ó ): (rÓe MrÓ )M(8Óe M8Ó ): (rÓ: MrÓ )
(8Ó: M8Óe )(8Óe M8Ó )(8Ó: M8Ó )
B=
From A and B, the quadratic equation Í() = A(x − ª ) ; + B(x-ª ) + oª = 0 give the next
approximation
ª_Y − ª =
Më±Êë : MAærÓ
;æ
…………….(2.26)
A direct solution from eq.(2.26) leads to inaccurate results. And therefore it is usually written in
the form
ª_Y − ª =−
;rÓ
ë±Êë : MAærÓ
……………..(2.27)
Note that:
1. If B>0, we use the positive sign with square root of equation and if B < 0, we use negative
sign with square root of the equation.
2. Once ª_Y is determined, the process is repeated. This brings up the issue of which point
is discarded. Two general strategies are typically used:
(i) If only real roots are being located, we choose the two original points that are nearest to
the new root estimate, ª_Y .
(ii) If both real and complex roots are being evaluated, a sequential approach is employed.
That is, just like the secant method, ªMY , ª , and ª_Y take the place of ªM; , ªMY , and ª .
= 2
(8Ó: M8Ó ): (rÓe MrÓ )M(8Óe M8Ó ): (rÓ: MrÓ )
(8Ó: M8Óe )(8Óe M8Ó )(8Ó: M8Ó )
B=
=5
The next approximation to the desired root is
ª_Y = ª −
;rÓ
ë±Êë : MAærÓ
,
=2−
;Y
C±Ê;CMA(;)(Y)
=2−
;Y
C_Ê;CMA(;)(Y)
[taking +ve sign, b/c B > 0 ]
= 1.780776
The procedure can now be repeated with three approximations as 1, 2, and 1.780776.
Let ªM; = 1, ªMY = 2, ª = 1.780776
Then oªM; = −2, oªMY = 1, oª = −0.304808
Now, using the formula like in the above we obtain
A=3.780773 and B= 5.123098
The next approximation to the desired root is
ª_Y = ª −
;rÓ
ë±Êë : MAærÓ
= 1.780776 −
;(MT.>TAVTV)
C.Y;>TUV_Ê(C.Y;>TUV): MA(>.XVTXX>)(MT.>TAVTV)
= 1.837867
The procedure can repeat with three approximations as 2, 1.780776, and 1.837867.
Let ªM; = 2, ªMY = 1.780776, ª = 1.837867
ª_Y = ª −
;rÓ
ë±Êë: MAærÓ
= 1.837867 −
;([Link])
[Link];;C_Ê([Link];;C): MA([Link];A )([Link])
= 1.839284
The procedure can now be repeated with three approximations as 1.780776, 1.837867, and
1.839284.
Let ªM; = 1.780776, ªMY = 1.837867, ª = 1.839284
Then oªM; = −0.304808, oªMY = −0.007757, oª =– 0.000015
Now, using the formula like in the above we obtain
A=4.20000 and B= 5.20000
The next approximation to the desired root is
ª_Y = ª −
;rÓ
ë±Êë : MAærÓ
,
= 1.839284 −
;(– [Link])
= 1.839287
Hence the required root is 1.839287.
EXAMPLE 2.20: Using Muller’s method, find the root of the equation () = 3– 2 – 5 = 0,
which lies between 2 and 3.
SOLUTION: Let ªM; = 1.9, ªMY = 2, ª = 2.1
Then oªM; =– 1.941, oªMY = −1, oª = 0.061
Now, we obtain
(8Ó: M8Ó )(rÓe MrÓ )M(8Óe M8Ó )(rÓ: MrÓ )
(8Óe M8Ó: )(8Óe M8Ó )(8Ó: M8Ó )
A=
([Link];.Y )([Link])M(;M;.Y)([Link] )
(;MY.U)(;M;.Y)([Link];.Y)
=
=6
(8Ó: M8Ó ): (rÓe MrÓ )M(8Óe M8Ó ): (rÓ: MrÓ )
(8Ó: M8Óe )(8Óe M8Ó )(8Ó: M8Ó )
B=
= 11.21
The next approximation to the desired root is
ª_Y = ª −
;rÓ
ë±Êë : MAærÓ
= 2.1 −
;([Link])
YY.;Y±Ê(YY.;Y): MA(E)([Link])
= 2.094542
The procedure can now be repeated with three approximations as 2, 2.1 and 2.094542.
Let ªM; = 2, ªMY = 2.1, ª = 2.094542
Then oªM; = −1, oªMY = 0.061, oª =– 0.0001058
Now, using the formula like in the above we obtain
A=6.194492 and B= 11.161799
The next approximation to the desired root is
ª_Y = ª −
;rÓ
ë±Êë : MAærÓ
= 2.094542 −
;(– [Link])
= 2.094551
Hence the required root is 2.0945 correct up to 4 decimal places.
2. Determine the solution of the equation 8 – (– §¬ ) = 0 by using the Bisection method
U
;
accurate to five decimal places in the interval (2, 3). Use the Bisection method to compute
the root of 7 8 – 3 = 0 correct to three decimal places in the interval (1.5, 1.6).
3. Find the root of R« = ¦R§ correct to two decimal places using Bisection method.
4. Use the Bisection method to find a root of the equation 3– 4– 9 = 0 in the interval
(2, 3), accurate to four decimal places.
5. Use Bisection method to find the smallest positive root of ¦R§ =(1/2) + §¬ in the
interval (0.41, 043). Use an error tolerance of = 0.0001.
7. Use the method of False Position to find solution accurate to within 10–4 for the function
() = – ¦R§ in the interval (0, π/2).
8. Use the method of False Position to find solution accurate to within 10–4 for the function
() = – 0.8 – 0.2 §¬ = 0 in the interval (0, π/2).
9. Use the method of False Position to solve the equation ¤ + 1 =0 accurate to three
decimal places starting with 2.5 and 3.0 as the initial approximations to the root.
10. Use the method of False Position to solve the equation 7 8 – ¦R§ =0 correct to four
decimal places in the interval (0, 1).
11. Use the method of False Position to find a root correct to three decimal places of the
function () = 3– 4– 9 = 0. ).
12. Find a root of the equation ¤ – ¤ℎ = 0 using the secant method with initial
approximations: T = 7 and Y = 7.5.
13. Find a root of the equation 6–– 1 = 0 using the secant method approximations: T =2 and
Y = 1.0.
14. Use the method of fixed point iteration to determine a solution accurate to within 10–2 for
4– 3 2– 3 =0 on [1, 2]. Use T =1.
15. Use the method of fixed point iteration to find a root of the equation 7 8 – 3 =0 in the
interval (0, 1) accurate to four decimal places.
Chapter objectives
At the end of this chapter you will be able to:
Describe the system of linear equations.
Identify if the system of linear equation has a unique solution.
Define direct and iterative methods.
State some of direct methods used for solving system of linear
equations.
Understand how we use the different direct methods.
Understand the relation and rate of convergence of each method.
State the difference between direct and iterative methods.
State the sufficient condition for iterative methods.
Compare the rate of convergence of those iterative methods.
Find the nth approximation formula for iterative methods.
INTRODUCTION
Consider a system of n linear algebraic equations in n unknowns Y , ; , …, b :
YY Y + Y; ; + ⋯ + Yb b = Y
;Y Y + ;; ; + ⋯ + ;b b = ;
⋮ ⋮ ⋯ ⋱ ⋮
………….(3.1)
bY Y + b; ; + ⋯ + bb b = b
In matrix notation we write the system as
AX =b …………(3.2)
Where
ða ó ð x2 ó ð b2 ó
a a .... a x b
å=ï
11 12 1n 1 1
ï M M ò ï M ò ï M ò
21 a 22
..... a2n
, X= and b=
M O
î an1 a n2
.... a nn ñ î xn ñ î bn ñ
The matrix [A | b], obtained by appending the column b to the matrix A is called the
Augmented matrix. That is
(A: b ) =ï
ï ò
.... M ò
21 a 22
..... 2n 2
.... ....
î a n1 a n2
.... a :b
nn n ñ
Linear systems of equations are associated with many problems in engineering and science, as
well as with applications of mathematics to the social sciences and quantitative study of business
and economic problems.
If the column vector = 0, then the system in eq. (3.1) is called homogeneous. If b≠ 0, then the
system in eq (3.1) is called non-homogeneous. i.e at least one b is non-zero.
A system of linear equations in n unknowns has a unique solution, provided that the
determinant of the coefficient matrix is non-singular i.e., if |A| ≠0. The rows and columns of a non-
singular matrix are linearly independent in the sense that no row (or column) is a linear
combination of the other rows (or columns).
Definition 3.1( Elementary row operations): Reducing matrices to nice and simple forms if any of
the following elementary row operations is performed.
i) Interchanging two rows:ª ↔ õ
ii) Multiplying a row by a non-zero number: ª → ª
iii) Add a multiple of another row to a given row:ª + õ → ª
Definition 3.2: A matrix is said to be in echelon form (also called a row echelon matrix) if it
satisfies the following three conditions:
a) All zero rows (rows consisting entirely of zeros) are at the bottom.
b) Each leading non- zero entry of a row in a column is to the right of the leading entry of
a row above it.
c) All entries in a column below a leading entry are zero.
EXAMPLE 3.2.: The following matrices are in row reduced echelon form
ð0 0 1 0ó 1 2 0 0
1 2 0 0
A=ï , ø = ö0 0 1 0÷
0 0 0 1ò 0 0 0 1
î 0 0 0 0 ñ
And the following are not in row reduced echelon form
1 0 5 1 0 5 1 0 0
ù = ö 0 1 3÷ ú = ö0 2 1÷ J = ö0 0 0÷
0 0 1 0 0 0 0 0 1
Definition 3.4 (Rank) The number of non-zero rows in the echelon form of a matrix is called the
rank of a matrix A , it is denoted by rank(å) or û(å).
1 2 1
EXAMPLE 3.3: Compute the echelon form reduced echelon form and rank of
B=ö2 3 1÷
1 1 2
SOLUTION:
1 2 1 ||: M;| e →|: 1 2 1
= M|e →|=
B=ö2 3 1÷ üýýýýýýýþ ö0 −1 −1÷
1 1 2 0 −1 1
|= _|: →|= 1 2 1
üýýýýýýþ ö0 −1 −1÷.
0 0 2
This is the echelon form. Therefore rank (B) =3.
Y
| →|=
; =
M|: →|:
1 2 1
üýýýýþ ö0 −1 −1÷
0 0 2
|: _|= →|:
|e M|= →|e
1 2 1
üýýýýýýþ ö0 1 −1÷
0 0 1
|e M;|: →|e 1 2 0
üýýýýýýýþ ö0 1 0÷
0 0 1
Prepared by Haileyesus T., UOG Page 66
1 0 0
ö0 1 0÷
0 0 1
This is the reduced echelon form of matrix B.
The solution of the system in eq. (3.1) depends on the consistent of the system and this consistent
is checked using rank.
The system of eq. (3.1) is consistent (has at least one solution), if
Rank (A) =Rank [A | b] = r
If r = n, then the system has unique solution.
If r < n, then the system has (n – r) parameter family of infinite number of solutions.
a) The system of eq. (3.1) is inconsistent(has no solution) if
Rank (A) ≠Rank [A | b].
Linear algebraic equations occur in almost all branches of engineering. Their most important
application in engineering is in the analysis of linear systems (any system whose response is
proportional to the input is deemed to be linear). Linear systems include structures, elastic solids,
heat flow, seepage of fluids, electromagnetic fields and electric circuits i.e., most topics taught in
an engineering curriculum. If the system is discrete, such as a truss or an electric circuit, then its
analysis leads directly to linear algebraic equations.
Here are two classes of methods for solving system of linear equations:
direct methods and
Iterative methods.
Direct methods: The common characteristics of these methods are that they transform the original
equation into equivalent equations (equations that have the same solution) using elementary row
operations that can be solved more easily. These methods take finite steps and give an exact
solution to the system.
Some of direct Methods:
a) matrix inversion method
b) Cramer's rule
c) Gauss Elimination Method
d) Gauss-Jordan Method
e) Lu-decomposition method
f) Cholesky’s method
g) partition method
If the system of equations has some special forms, then the solution is obtained directly.
We consider three such special forms.
Let A be a diagonal matrix, A =D. That is, we consider the system of equations in eq.(2)
become Dx = b .i,e
YY Y = Y
;; ; = ;
⋱ ⋮
bb b = b
This system is called a diagonal system of equations. Solving directly, we obtain
Let A be a lower triangular matrix, A = L. That is, we consider the system of equations in
eq.(2) become Lx =b . i.e
YY Y = Y
;Y Y + ;; ; = ;
⋮ ⋮ ⋯ ⋱ ⋮
bY Y + b; ; + ⋯ + bb b = b
This system is called a lower triangular system of equations. Solving for the unknowns in the order
x1,x2,….,xn ,we get
; =
: M?:e 8e
?::
.
……………………
f M∑fe
e ?f 8
b =
?ff
.
The unknowns are obtained by forward substitution and this procedure is called the forward
substitution method.
Let A be an upper triangular matrix, A = U. That is, we consider the system of equations in
eq.(2) become Ux =b . i.e
YY Y + Y; ; + ⋯ + Yb b = Y
;; ; + ⋯ + ;b b = ;
⋱ ⋮
bb b = b
This system is called an upper triangular system of equations. Solving for the unknowns in the
order xn,xn-1,….,x1 ,we get
b =
f
?ff
,
fe M?fe 8f
bMY =
?fefe
,
…………………………
e M∑f
: ?e 8
b =
?ee
,
The unknowns are obtained by back substitution and this procedure is called the back substitution
method.
Mat lab cod for upper triangular system
%%M-file: Usolve.m
function X= Usolve(U,b)
%input: upper triangular matrix U, right hand side vector b
%output: solution X
n=max(size(U)); %determines the dimension of the matrix U
X(n) = b(n)/U(n,n);
for k=n-1:-1:1 %Backward for loop
X(k) = ( b(k) - sum( U(k, k+1:n).*X(k+1:n) ))/U(k,k);
end
Ax=b …………(3.4)
Ax= b
1 1 1 1
A=ö1 2 3÷, x= o and b=ö6÷
1 3 4 p 6
where
Ax=b
A=ö ;Y ;; ;> ÷ , x = ö; ÷ and b= ö; ÷
>Y >; >> > >
If the coefficient matrix of the system is nonsingular, then the system has a unique solution. And
the solution of the system is given as follows:
1 3 6 2
SOLUTION: The coefficient matrix A and the column vector b is given as
A= ö3 −1 4 ÷ 9÷
1 −4 2 7
and b=ö
|å|= 1L−1 4 3 4 3 −1
Then let us find the determinant of A and others
L-3L L+ 6L L= -58
−4 2 1 2 1 −4
2 3 6 1 2 6 1 3 2
∆= 9 −1 4 =-116, ∆o= 3 9 4 =58, ∆p= 3 −1 9 =-29
7 −4 2 1 7 2 1 −4 7
From this, the solution of the system is
∆8 ∆r ∆n
= |æ| = 2, o = |æ| = −1 and p = |æ| = 0.5
Step 2: To eliminate x2 from the third equations using the second equation, following operations
are performed:
ª − ; → ª for i=3
?Ó:
?::
After these two steps, we will get an equivalent system of linear equation given as:
YY Y + Y; ; + Y> > = Y .
;; ′; + ;> ′> = ; ′
>> ′′> = > ′′
Step 3: To findY , ; and > , we apply back substitution
EXAMPLE 3.6: Solve the system of equation by Gauss-Elimination method.
2x + 3 y − z = 5
4 x + 4 y − 3z = 3
2x − 3 y + 2z = 2
3 −1 x
SOLUTION: We can write the system as
2 5
ö4 4 −3÷ = ö3÷
y
2 −3 2 z 2
Now let us apply elementary row operations on the augmented matrix
−3 2 ∶ 2
|= M|e →|=
2
2 3 −1 ∶5
> + 3; → > ö0 2 1 ∶ 7÷
0 −6 3 ∶ 3
2 3 −1 ∶ 5
ö0 2 1 ∶ 7÷
0 0 6 ∶ 18
Now the corresponding equivalent system of linear equation is
2x + 3y − z= 5
2y + z= 7
6 z = 18
3 3 4 ∶ 20
SOLUTION: The augmented matrix is given b
ö2 1 3 ∶ 13÷
1 1 3 ∶6
3 3 4 ∶ 20
Now let us apply the an elementary row operation on to the augmented matrix
→ Y ö2 1 3 ∶ 13÷
Y
> Y
1 1 3 ∶6
1 1 4/3 ∶20/3
ö2 1 3 ∶ 13 ÷
|: M;|e →|:
1 1 3 ∶ 6
|= M|e →|=
1 1 4/3 ∶ 20/3
ö0 −1 1/3 ∶ −1/3÷
0 0 5/3 ∶−2/3
Back substitution gives the solution as
EXAMPLE 3.8: Solve the system of equations using the Gauss elimination method
x + 10 y − z = 3
2 x + 3 y + 20 z = 7
9 x + 22 y + 79 = 45
1 10 −1 ∶ 3
SOLUTION: We have an augmented matrix
ö2 3 20 ∶ 7 ÷,
9 22 79 ∶45
10 −1 ∶ 3
Applying the elementary row operation on to the augmented matrix, we have
1
ö23 20 ∶ 7 ÷
|: M;|e →|:
22 79 ∶45
|= MU|e →|=
9
1 10 −1 ∶ 3
> − 4; → > ö0 −17 22 ∶ 1 ÷
0 −68 88 ∶18
1 10 −1 ∶ 3
ö0 −17 22 ∶ 1 ÷
0 0 0 ∶ 14
But we obtain that the third row elements in the coefficient matrix are all zero this leads us that
the rank of the coefficient matrix and augmented matrix have different value( i.e.
R(A)=2≠3=R(AG)), consequently the system is inconsistent and hence has no solution.
Mat lab Code 0f Gaussian elimination
function [x] = GaussianEliminate(A, b)
%work out the number of equations
N = length(b)
%Gaussian elimination
for column=1:(N-1)
%work on all the rows below the diagonal element
for row = (column+1):N
%work out the value of d
d = A(row,column)/A(column,column);
%do the row operation
A(row,:) = A(row,:)-d*A(column,:)
b(row) = b(row)-d*b(column)
end%loop through rows
end %loop through columns
−1 ∶3
SOLUTION: We have the augmented matrix
1 10
ö2 3 20 ∶ 7÷
10 −1 2 ∶4
Since we must apply partial pivoting we perform the following elementary row transformations
10 −1 2 ∶ 4
and do the eliminations
Y ↔ > ö 2 3 20 ∶ 7÷
1 10 −1 ∶3
Then apply ; − → ; and > − → > ,we get
|e |e
C YT
10 −1 2 ∶ 4
ö0 3.2 19.6 ∶ 6.2÷
0 10.1 −1.2 ∶2.6
10 −1 2 ∶ 4
Since we must apply partial pivoting her also, we interchange R2 and R 3 Then we get
3.2;
Next we apply an elementary operation to eliminate y from the third row .Apply
> − → >
10.5
10 −1 2 ∶ 4
Then we get
REMARK: Gauss elimination method is a direct method. Therefore, it is possible to count the
total number of operations, that is, additions, subtractions, divisions and multiplications. Without
going into details, in the system of linear equations the total number of divisions and
multiplications (division and multiplication take the same amount of computer time) is n (n2+ 3n–
x+ 2y − 3z =4 3x + 5 y + 2 z =4
c) − 3y + − 4 y + 7x = 5 d) 8 y + 2 z = −7
6 z + 8x − 9 y = 1 6 x + 2 y + 8 z = 26
1 1 1 ∶1
SOLUTION: We have the augmented matrix a
ö4 3 −1 ∶6÷
3 5 3 ∶4
Now let us apply an elementary operation on the augmented matrix
; − 4Y → ; 1 1 1 ∶ 1
Apply
ö4 3 −1 ∶6÷
> − 3Y → >
3 5 3 ∶4
ö0 −1 −5 ∶2÷
Y + ; → Y
0 2 0 ∶1
1 0 −4 ∶ 3
Apply
−1
> → > ö0 −1 −5 ∶ 2÷
10
0 0 −10 ∶5
5> → ;
ö0 −1 −5 ∶ 2 ÷
Y + 4> → Y
0 0 1 ∶ −1/2
0 ∶ 1
Apply
1 0
−; → ; ö0 −1 0 ∶ −1/2÷
0 0 1 ∶−1/2
1 0 0 ∶ 1
ö0 1 0 ∶ 1/2 ÷ , b 6=
0 0 1 ∶−1/2
[ w/h is linear system of the form b' ]
1 3 2 ∶ 17
SOLUTION: We have the augmented matrix a
ö1 2 3 ∶ 16÷
2 −1 4 ∶13
Now let us apply an elementary operation on the augmented matrix
2 ∶ 17
Apply
; − Y → ; 1 3
ö1 2 3 ∶ 16÷
> − 2Y → >
2 −1 4 ∶13
ö0 −1 1 ∶ −1 ÷
Y + 3; → Y
0 −7 0 ∶−21
5 ∶ 14
Apply
−1 1 0
> → > ö0 −1 1 ∶ −1 ÷
7
0 0 −7 ∶ −14
5 ∶ 14
Apply
; − > → ; 1 0
ö0 −1 1 ∶ −1÷
Y − 5> → Y
0 0 1 ∶2
Apply
1 2 2
EXAMPLE 3.13: Find the inverse of the matrix using the Gauss-Jordan method.
ö2 3 2 ÷
3 2 1
1 2 2 ∶ 1 0 0
SOLUTION: Consider the augmented matrix
ö2 3 2 ∶ 0 1 0 ÷
3 2 1 ∶ 0 0 1
We perform the following elementary row operations and do the eliminations
2 ∶ 1 0
Apply
; − 2Y → ; 1 2 0
ö2 3 2 ∶ 0 1 0 ÷
> − 3Y → >
3 2 1 ∶ 0 0 1
ö0 −1 −2 ∶ −2 1 0 ÷
Y + 2; → Y
0 −4 −5 ∶ −3 0 1
1 0 −2 ∶ −3 2 0
Apply
1
→ > ö0 −1 −2 ∶−2 1 0 ÷
3 >
0 0 3 ∶ 5 −4 1
0 −2 ∶ −3 2 0
Apply
1
−; → ; ö0 1 2 ∶ 2 −1 0 ÷
0 0 1 ∶ 5/3 −4/3 1/3
ö0 1 2 ∶ 2 −1 0 ÷
|: M;|= →|:
1 2 1
EXAMPLE 3.14: Find the inverse of the matrix using the Gauss-Jordan method.
ö2 3 −1 ÷
2 −1 3
1 2 1 ∶ 1 0 0
SOLUTION: Consider the augmented matrix
ö2 3 −1 ∶ 0 1 0 ÷
2 −1 3 ∶ 0 0 1
We perform the following elementary row operations and do the eliminations
; − 2Y → ; 1 2 1 ∶ 1 0 0
Apply
ö2 3 −1 ∶ 0 1 0 ÷
> − 2Y → >
2 −1 3 ∶ 0 0 1
ö0 −1 −3 ∶−2 1 0 ÷
Y + 2; → Y
0 −5 1 ∶ −2 0 1
1 0 −4 ∶ −3 2
Apply
1 0
> → > ö0 −1 −2 ∶−2 1 0 ÷
16
0 0 16 ∶ 8 −5 1
; + 2> → > 1 0 −4 ∶ −3 2 0
Apply
ö0 −1 −2 ∶ −2 1 0 ÷
Y + 2> → Y
0 0 1 ∶ 1/2 −5/16 1/16
Apply
1 1 1 1 2 3
3. Find the inverse of the matrix using the Gauss-Jordan method.
ö4 3 −1 ÷ b) ö 2 3 1 ÷
3 5 3 −1 0 2
a)
Y;
|YY | ≠ 0 ,LYY ;; ;> ≠ 0
;; L ≠ 0 , ;Y
;Y
>Y >; >>
It is possible to show that any square matrix A can be expressed as a product of a lower triangular
matrix L and an upper triangular matrix U.
A=LU
2 1 4 12
SOLUTION: We have
2
UX=K
1 1/2 6
ö0 1 2 ÷ o 4÷
0 0 1 p 1
=ö
2 1 1 7
SOLUTION: We have
Hence
2 0 0 1 1/2 1/2
=ö1 3/2 0 ÷ and = ö0 1 1/3÷
1 1/3 4/3 0 0 1
Now
LK=b
2 0 0 Y 7
ö1 3/2 0 ÷ ö; ÷= ö8÷
1 1/3 4/3 > 9
From this, by forward substitution, we get
Y = 3.5,; = 3,> = 3
And then, to get the solution we use
1/2
UX=K
1 1/2 3.5
ö0 1 1/3÷ o =ö 3 ÷
0 0 1 p 3
By backward substitution, the solution is given as:
= 1, o = 2, p = 3
EXAMPLE 3.17: Solve the following set of equations by using the LU-Decomposition method
x+ y−z =2
2 x + 2 y + 5 z = −3
3x + 2 y − 3z = 6
SOLUTION: We have
1 1 −1 2
A = ö2 2 5 ÷, X=o and b = ö−3÷
3 2 −3 p 6
Let A = LU , where
YY 0 0 1 lY; Y>
L=ö;Y ;; 0÷ U=ö0 1 ;> ÷
>Y >; >> 0 0 1
Then,
3 0 0
Hence
1 2/3 −1
= ö2 3/2 0 ÷ and = ö0 1 21/2÷
1 1/3 −7/3 0 0 1
Now
LK=b
2
UX=K
1 2/3 −1
ö0 1 21/2÷ o =ö−14/3÷
0 0 1 p 2/3
By backward substitution, the solution is given as:
= 82/9, o = −35/3, p = 2/3
EXERCIS 3.4
Solve the following set of equations by using the LU-Decomposition method
2x + 2 y + z = 1 2x + y + z + w = 2
a) 4 x + 3 y + 3z = 2 4 x + 0. y + 2 z + w = 3
x+ y+z =3 3x + 2 y + 2 z + 0 w = −1
b) x + 3 y + 2 z + 6w = 2
A =
and positive definite, then the matrix A can be decomposed as
………. (3.11)
YY 0 0
Where = ö;Y ;; 0 ÷ is a lower triangular matrix ( i.e L = (lij), lij =0,i<j.
>Y >; >>
A =
Alternatively, A may be decomposed as
SOLUTION: We have
1 2 3 5
A = ö2 8 22÷, X=o and b = ö 6 ÷
3 22 82 p −10
Let A = , where
YY 0 0 YY Y; Y>
L=ö;Y ;; 0÷ and =ö 0
;; ;> ÷
>Y >; >> 0 0 >>
Then
1 2 3 YY 0 0 YY Y; Y>
ö2 8 22÷= ö;Y ;; 0 ÷ö 0 ;; ;> ÷
3 22 82 >Y >; >> 0 0 >>
YY ; YY Y; YY Y>
= ;Y YY ;Y ; + ;; ; ;Y Y> + ;; ;>
>Y YY >Y Y; + >; ;; >Y ; + >; ; + >> ;
YY = 1, ;Y YY = 2 → ;Y= 2, >Y YY= 3→ >Y = 3
Y; ; + ;; ; = 8 → ;; =8−4 = 2,
>Y Y; + >; ;; =22→ >; =8
>Y ; + >; ; + >> ; =82 → >> =3
Hence, we get
=
Then using eq. (3.14)
1 0 0 Y 5
ö2 2 0÷ ö; ÷ = ö 6 ÷
3 8 3 > −10
From this, solving by forward substitution, we get Y = 5, ; = -2 and > = -3. Finally using
6 =
1 2 3 5
ö0 2 8÷ = ö−2÷
o
0 0 3 p −3
Solving this by backward substitution, we obtain the solution
= 2, o = 3, p = −1
EXAMPLE 3.19: Solve the following set of equations by using the Cholesky's method.
2x + y + z = 7
x + 2y + z = 8
x + y + 2z = 9
SOLUTION: To make our work simple, let us rewriting the example as
x + 2y + z = 8
2x + 5 y + z = 7
x + y + 6z = 9
We have
1 2 1 8
A = ö2 5 1÷, o ö 7÷
1 1 6 p 9
X= and b =
where B, C, E, D, are of orders QQ, Q§, §Q and §§ respectively, with Q + § = . Similarly,
we partition åMY as
where X, Y, Z and V are of the same orders as B, C, E and D respectively. Using the identity
We obtain
V= (D − EBMY C)MY , Y= −B MY CV,
Z=-VEB MY , X=−B MY −BMY CZ
where we have assumed that !M exists. If !M does not exist but "M exists then the equationscan
be modified suitably. This procedure requires finding the inverse of two lower order matrices,
!M and (# − $!M %)M .
1 1 1
EXAMPLE 3.21: Determine the inverse of the matrix
ö4 3 −1÷
3 5 3
using the partition method. Hence, find the solution of the system of equations
x + y + z =1
4x + 3 y − z = 6
3x + 5 y + 3 z = 4
SOLUTION: Let the matrix Abe partitioned as
Now,
3 −1 MY 1 4
Y= −ø MY CV = − (− ) = (− )
Y Y
−4 1 −1 YT YT −5
3 −1 MY
Z=-VEø MY = (− ) (3 5) = − (-11 2)
Y Y
YT −4 1 YT
X =−ø MY −ø MY CZ
−3 1 −3 1 1
= − (− ) (-11 2)
Y
4 −1 YT 4 −1 −1
1.4 0.2
=
−1.5 0
The solution of the given system of equations is
1.4 0.2 −0.4 1 1
X=o = ö−1.5 0 0.5 ÷ ö6÷= ö 0.5 ÷
p 1.1 −0.2 −0.1 4 −0.5
So far, we have discussed only direct methods of solution. The common characteristic of these
methods is that they compute the solution with a finite number of operations. Moreover, if the
computer were capable of infinite precision (no round off errors), the solution would be exact.
Iterative, or indirect methods, start with an initial guess of the solution x and then repeatedly
improve the solution until the change in x becomes negligible. Since the required number of
iterations can be large, the indirect methods are, in general, slower than their direct counter parts.
However, iterative methods do have the following advantages that make them attractive for certain
problems:
a) Possible to deal with very large matrices that are sparse, but not necessarily
banded.
b) Iterative procedures are self-correcting, meaning that round off errors (or
even arithmetic mistakes) in one iterative cycle are corrected in subsequent
cycles.
Iterative methods are based on the idea of successive approximations. We start with an initial
approximation to the solution vector = (T) , to solve the system of equations Ax = b, and
obtain a sequence of approximate vectors (Y) , (;) ,…, (b) . We stop the iteration procedure when
the magnitudes of the differences between the two successive iterates is significantly small.
The sufficient condition for the use of these methods is that the diagonal elements of the
coefficient matrix must be dominant. That means if ªª , i=1,2,3,…n are diagonal elements then
|ªª | ≥ ∑bõ)1(ªõ ( i=1,2,…,n
ª*õ
Diagonal dominant condition is a sufficient condition for convergence but it is also a necessary
condition. Under the category of iterative method, we shall describe the following two methods:
i) Jacobi’s method
ii) Gauss-Seidel method.
REMARK: A serious drawback of iterative methods is that they do not always converge to the
solution and it can be shown that convergence is guaranteed only if the coefficient matrix is
diagonally dominant. In the absence of any better estimate we take the initial approximation
(T , oT , pT ) = (0,0,0) for the iterative methods
… … … … … … … … … … … … … … … … … … … … …. ……………..(3.16)
b = b − bY Y − b; ; − ⋯ − bbMY bMY
Y
?ff
…………………………………………………
b = ? b − bY Y − b; ; − ⋯ − bbMY bMY
(Y) Y (T) (T) (T)
ff
… … … … … … … … … … … … … … … … … … …,
b = b − bY Y − b; ; − ⋯ − bbMY bMY
(;) Y (Y) (Y) (Y)
?ff
… … … … … … … … … … … … … … … … … … …,
EXAMPLE 3.22: Solve the following system of equation using Jacobi’s method
5 x − y + z = 10
2 x + 4 y = 12
x + y + 5 z = −1
Start with the solution (2, 3, 0).
SOLUTION: The coefficient matrix of the given system is diagonally dominant hence we can
write the system as:
= C 10 +o − p
Y
o = A 12 − 2
Y
p = C −1 − − o
Y
Second approximation:
(;) = 10 +o (Y) − p (Y) =2.64
Y
C
Third approximation:
(>) = C 10 +o (;) − p (;) = 2.564
Y
Fifth approximation:
(C) = C 10 +o (A) − p (A) = 2.5534
Y
Sixth approximation:
(E) = C 10 +o (C) − p (C) = 2.5557
Y
Seventh approximation:
(X) = 10 +o (E) − p (E) = 2.5558
Y
C
Hence, approximating solution after having some other approximations up to 3 decimal places is
= 2.556, o = 1.722, p = – 1.055.
EXAMPLE 3.22: Solve the system
3 x + 4 y + 15 z = 54.8
x + 12 y + 3 z = 39.66
10 x + y − 2 z = 7.74
using Gauss-Jacobi's method.
SOLUTION: The coefficient matrix of the given system is not diagonally dominant. Hence, we
rearrange the equations such that the elements in the coefficient matrix are diagonally dominant .so
that the system can be written as:
1
Now to apply the Gauss-Jacobi's method let us write this system in the form:
= 7.74 − o + 2p
10
1
o= 39.66 − − 3p
12
1
p= 54.8 − 3 − 4o
15
Now, since there is no initial approximation given, we start with initial the
approximation,( (T) , o (T) , p (T) ) = (0, 0, 0):
First approximation:
(Y) = YT 7.74 − o (T) + 2p (T) =0.774
Y
Second approximation:
(;) = YT 7.74 − o (Y) + 2p (Y) =1.1742
Y
Third approximation:
(>) = 7.74 − o (;) + 2p (;) =1.0647
Y
YT
Fourth approximation:
(A) = 7.74 − o (>) + 2p (>) =1.0879
Y
YT
Sixth approximation:
(E) = YT 7.74 − o (C) + 2p (C) =1.0751
Y
Seventh approximation:
(X) = YT 7.74 − o (E) + 2p (E) =1.0748
Y
Therefore from the 6th and 7th approximation, the solution of the system is = 1.075, o =
2.524, p = 2.765 approximated up to three decimal places.
EXAMPLE 3.23: solve the system
26 x + 2 y + 2 z = 12.6
3 x + 27 y + z = −14.3
2 x + 3 y + 17 z = 6
using Gauss-Jacobi's method.
SOLUTION: The coefficient matrix of the given system is diagonally dominant
1
Now to apply the Gauss-Jacobi's method let us write this system in the form:
= 12.6 − 2o − 2p
26
1
o= −14.3 − 3 − p
12
1
p= 6 − 2 − 3o
15
Now, since there is no initial approximation given, we start with initial the approximation
,( (T) , o (T) , p (T) ) = (0, 0, 0):
Second app
Second approxima
roximation:
pproximation:
(;) = 12.6 − 2 (Y) − 2p (Y) = 0.49821
Y
;E
Third app
Third approxima
roximation:
pproximation:
(>) = ;E 12.6 − 2o (;) − 2p (;) = 0.50006
Y
Fourth app
Fourt approxima
roximation:
pproximation:
(A) = ;E 12.6 − 2o (>) − 2p (>) = 0.50000
Y
Fifth app
Fift approxima
roximation:
pproximation:
(C) = ;E 12.6 − 2o (A) − 2p (A) = 0.50001
Y
We find that
( (C) − (A) ( =0.00001, (o (C) − o (A) ( =0.00001, (p (C) − p (A) (= 0.00011
Since, all the errors in magnitude are less than 0.0005, the required solution is
x = 0.5, y = -0.6, z = 0.4
Remark: What is the disadvantage of the Gauss-Jacobi method?. At any iteration step, the value
of the first variable Y is obtained using the values of the previous iteration. The value of the
… … … … … … … … … … … … … … … … … … … … …. ……………..(3.18)
b = b − bY Y − b; ; − ⋯ − bbMY bMY
Y
?ff
1
b = ¼ − bY Y − b; ; − ⋯ − bbMY bMY ½
(Y) (Y) (Y) (Y)
bb b
In general, ( + 1)th iteration can be written as:
1
Y = − Y; ; − Y> > − ⋯ − Yb b
(b_Y) (b) (b) (b)
YY Y
; = ; − ;Y Y − ;> > − ⋯ − ;b b
(b_Y) Y (b_Y) (b) (b)
?::
… … … … … … … … … … … … … … ………………………
This process is continued until the value of Y , ; , … , b are obtained to the desired degree of
accuracy. That is when the difference between ( + 1)th and ()th iterations is significantly small.
EXAMPLE 3.24: Solve the system of equations by Gauss-Seidel iteration method
8 x − 3 y + 2 z = 20
6 x + 3 y + 12 z = 35
4 x + 11y − z = 33
SOLUTION: The coefficient matrix of the given system is not diagonally dominant. Hence, we
rearrange the equations such that the elements in the coefficient matrix are diagonally dominant.
So that the system can be written as:
8 x − 3 y + 2 z = 20
4 x + 11y − z = 33
6 x + 3 y + 12 z = 35
1
Now to apply the Gauss-Jacobi's method let us write this system in the form:
= 20 + 3o − 2p
8
1
o= 33 − 4 + p
11
1
p= 35 − 6 − 3o
12
Second app
Second approxima
roximation:
pproximation:
(;) = 20 + 3o (Y) − 2p (Y) =2.9981061
Y
V
Third app
Third approxima
roximation:
pproximation:
(>) = 20 + 3o (;) − 2p (;) =3.0266228
Y
V
Fourth app
Fourt approxima
roximation:
pproximation:
(A) = 20 + 3o (>) − 2p (>) = 3.0165121
Y
V
Fifth app
Fift approxima
roximation:
pproximation:
(C) = 20 + 3o (A) − 2p (A) = 3.0166005
Y
V
Sixth app
Sixt approxima
roximation:
pproximation:
(E) = V 20 + 3o (C) − 2p (C) = 3.0167568
Y
Since at the sixth and seventh approximations, the values of , o, p are the same, correct to four
decimal places, we can stop the iteration process and the solution is = 3.0167, o = 1.9858, p =
0.9118.
EXAMPLE 3.25: Find the solution of the system of equations
45 x + 2 y + 3 z = 58
− 3x + 22 y + 2 z = 47
5 x + y + 20 z = 67
Correct to three decimal places, using the Gauss-Seidel iteration method.
SOLUTION: The given system of equations is strongly diagonally dominant. Hence, we can
expect fast convergence. Gauss-Seidel method gives the iteration
= AC 58 − 2o − 3p
Y
o = ;; 47 + 3 − 2p
Y
p= 67 − 5 − o
Y
;T
Starting with initial the approximation (T) =o (T) =p (T)= 0, we get the following results:
First app
First approxima
roximation:
pproximation:
(Y) = 58 − 2o (T) − 3p (T) =1.28889
Y
AC
Second app
Second approxima
roximation:
pproximation:
(;) = 58 − 2o (Y) − 3p (Y) = 0.99198
Y
AC
Third app
Third approxima
roximation:
pproximation:
(>) = 58 − 2o (;) − 3p (;) = 0.99958
Y
AC
Fourth app
Fourt approxima
roximation:
pproximation:
(A) = 58 − 2o (>) − 3p (>) = 1.00000
Y
AC
Fifth app
Fift approxima
roximation:
pproximation:
C = 58 − 2oA − 3pA = 1.00000
Y
AC
Since at the 4th and 5th approximations, the values of , o, p are the same, correct to four decimal
places, we can stop the iteration process and the solution is = 1, o = 2, p = 3.
EXAMPLE 3.26: Solve the system of equations by Gauss-Seidel iteration method
3x − 6 y + 2 z = 23
− 4 x + y − z = −8
x − 3 y + 7 z = 17
SOLUTION: The coefficient matrix of the given system is not diagonally dominant. Hence, we
rearrange the equations such that the elements in the coefficient matrix are diagonally dominant.
So that the system can be written as:
4 x + y − z = −8
3x − 6 y + 2 z = 23
x − 3 y + 7 z = 17
Now to apply the Gauss − Jacobi′s method let us write this system in the form:
1
= −8 − o + p
4
1
o= 23 − 3 − 2p
−6
1
p = 17 − + 3o
7
Now, since there is initial approximation given, we start with initial approximation
(T) =o (T) =p (T) = 0.
Second app
Second approxima
roximation:
pproximation:
(;) = −8 − o (Y) + p (Y) = 0.9953
Y
A
Third app
Third approxima
roximation:
pproximation:
(>) = A −8 − o (;) + p (;) =0.9988
Y
Fourth app
Fourt approxima
roximation:
pproximation:
(A) = A −8 − o (>) + p (>) =0.9998
Y
We find that
( (A) − (>) ( =0.0010, (o (A) − o (>) ( =0.0012,(p (A) − p (>) (= 0.0003.
Since, all the errors in magnitude are less than 0.005, the required solution is
= 0.9998, o =– 3.0002, p = 0.9999
Rounding to two decimal places, we get = 1. o = −3, p = 1
Mat lab cod of Gauss Seidel method
% A PROGRAM OF GAUSS SEIDEL METHOD
function x = GSsolve(A, b, tolerance)
n = size(A,1); % Get matrix size
x = zeros(n,1); % Result vector for iteration n
xn= zeros(n,1); % Result vector for iteration n+1
done = false; % Flag to exit iteration loop
ni = 0; % Number of iterations
while done==false
% Count the number of iterations
ni = ni+1;
% Do one iteration for each element of the unknown vector
While solving simultaneous equations, we also come across two forms of instabilities
Inherent: Inherent instabtity of the system is the property of the given problem and occurs
due to the problem being ill-conditioned. It can be avoided by reformations of the problem
suitably.
Induced: Induced instability occurs because of the incorrect choice of the method.
How do we know the linear system is ill-conditioned?
There are several ways to check a matrix A for ill-conditioning. If the magnitude of the
determinant of the matrix is small, the matrix is ill-conditioned. However, this is not a foolproof
test. The inverse matrix åMY can be calculated, and AåMY can be computed and compared to I.
Similary, (åMY )MY can be computed and compared to A. A close comparison in either case
suggests that matrix A is well-conditioned. A poor comparison suggests that the matrix is ill-
conditioned. Some of the elements of A and/or b can be changed slightly, and the solution
The 1-norm of a square matrix is the maximum of the absolute column sums.
n
The Infinity-Norm : A = max (∑ aij )
∞ j i =1
The infinity-norm of a square matrix is the maximum of the absolute row sums.
The Euclidean norm of a square matrix is the square root of the sum of all the squares of the
elements.
EXAMPLE3.30: Calculate the norm of
5 −4 2
B = ö−1 2 3÷
−2 1 0
SOLUTION: Let us see the norm in three different norms
= max (8, 7, 5)
=8
= max (11, 6, 3)
= 11
= √25 + 16 + 4 + 1 + 4 + 9 + 4 + 1 + 0
=√64 =8
CONDITION NUMBER
Sometimes, one comes across a system of equations which are very sensitive to round off errors.
That is, one gets different solutions when the elements are rounded to different number of digits. In
such cases, the system is called an ill-conditioned system of equations. The measure of the ill-
conditionedness is given by the value of the condition number of the matrix A.
Definition: Condition number of a system is a measure of the sensitivity of the system to small
changes in any of its elements.
To investigate the condition number of a linear algebraic system, we perturb the equation and then
examine the size of resulting change in the solution. First, consider the effect of change in the
Or
‖s8‖ ‖s‖
‖8‖
≤ (å) ‖‖
1 = I = A −1A ≤ A −1 ⋅ A
III) = Cond(A)
EXAMPLE 3.31: show that the system is ill-condition.
1 2 x 2
2 3.999 y = 4
SOLUTION: For
[A] =
2
1
2 3.999
A ∞ = 5.999
A −1 = 5999
∞
Cond ( A) = A ∞ A−1
∞
= 5.999 × 5999.4
= 35990
Which is so large, implies that the system is ill-conditioned.
EXAMPLE 3.32: Show that the system is ill-condition.
1 2 x 4
2 3 y = 7
SOLUTION: For
[A] =
1 2
2 3
It can be shown
−3
[A]−1 =
2
2 − 1
Then
A ∞ = 5,
A−1 = 5.
∞
Cond (A) = A ∞ A −1
∞
= 5× 5
= 25
Which is so large, implies that matrix A is ill-conditioned.
is given by x=2.8, y= 1, z= 1.8. Using the iterative method, improve this solution.
SOLUTION: Substituting the approximate value, we get
2 x − y + 3z = 4
b) x + 9 y + 3 z = −8
4 x − 8 y + 11z = 15
2. Solve the following system of linear equations using the method of Gaussian elimination.
2x + y + z = 4
a) 0 x + 3 y − 3z = 0
0x − y − 2z = 1
x + 2 y + 3 z + 4w = 8
2 x − 2 y − z − w = −3
b)
x − 3 y + 4 z − 4w = 8
2 x + 3 y − 3z + 4w = −2
3. Solve the following set of simultaneous linear equations by the Gauss-Jordan method.
4 x − 3 y + 5 z = 34
a) 2x − y − z = 6
x + y + 4 z = 15
2 x − y + z = −1
b) 3x + 3 y + 9 z = 0
3x + 3 y + 5 z = 4
4. Solve the following set of simultaneous linear equations using the LU decomposition
3x + 2 y + 7 z = 4
a) 2x + 3y + z = 5
3x − 4 y + z = 7
x − 4 y + 7z = 9
b) 2 x − 3 y + 4 z = 13
3x + y + 5 z = 40
5. Solve the following set of simultaneous linear equations using the Cholesky's method.
7. Solve the following system of simultaneous linear equations using the Gauss-Seidal method.
4 x − 3 y + 5 z = 34
a) 2x − y − z = 6
x + y + 4 z = 15
2 x − y + 5 z = 15
b) 2x + y + z = 7
x + 3 y + z = 10
8. An approximate solution of the system
x + 4 y + 7z = 5
2 x + 5 y + 8z = 7
3 x + 6 y + 9.1z = 9.1
is given by x=1.8, y= -1.2, z= 1.8. Using the iterative method, improve this solution
INTRODUCTION
The calculus of finite differences deals with the changes that take place in the value of the function
(dependent variable) due to finite changes in the independent variable. The study of finite
difference calculus has become very important due to its wide variety of application in everyday
life. It has been originated by Sir Issac Newton. It has been of great use for Mathematicians as
well as Computer Scientists for solution of the Scientific, Business and Engineering problems.
There, it helps in reducing complex mathematical expressions like trigonometric functions in terms
of simple arithmetic operations.
Numerical methods are very important tools to provide practical methods for calculating the
solution of problems to applied mathematics for a desired degree of accuracy.
If o = () is a function for ≤ ≤ , this means that one or more values of o = () exist
corresponding to every value of in the given range. However if the function is not known, the
value of o can be obtained, when a set of values of is given. The method to find out such values
is based on principle of finite differences provided the function is continuous.
If a function o = () assumes the values (), ( + ℎ), ( + 2ℎ), . .. corresponding to the
values of ( T = , Y = + ℎ, ; = + 2ℎ, … ) then each value of is called Argument and
Entry.
its corresponding values of o is called Entry
FORWARD DIFFERENCE
4.1 FORW
The forward difference operator is denoted by ∆. Let us see some of the first order forward
differences
∆oT = oY – oT
∆oY = o; – oY
…………………..
∆ob = ob_Y – ob
or ∆oª = oª_Y – oª ¬ = 0, 1, 2, … . , – 1
∆; oT = ∆oY – ∆oT
……………
x y ∆y ∆; o ∆> y ∆A o
T oT
∆oT
Y oY ∆; oT
∆oY ∆> oT
; o; ∆; oY ∆A oT
∆o; ∆> oY
> o> ∆; o;
∆o>
A oA
EXAMPLE 4.1: Construct a forward difference table for the following values
0 5 10 15 20 25
() 7 11 14 18 24 32
SOLUTION: For = 0, 1, 2, 3, 4, 5, we get the values of o are 1, 1, 9, 31, 73, 141. Therefore,
difference table for these data is as:
o ∆y ∆; o ∆> y ∆A o
0 1
0
1 1 8
8 6
2 9 14 0
22 6
3 31 20 0
42 6
4 73 26 0
8 6
5 141 32
100
6 241
⇒ oE − 141=100 ⇒ oE = 241
Let's confirm that for given function y = > + ; – 2 + 1, at = 6, 63+ 62– 2(6) + 1= 241.
EXAMPLE 4.3: Find (6) given the table below and the third differences being constant.
0 1 2 3
() -3 6 8 12
o ∆y ∆; o ∆> y
0 -3
9
1 6 -7
2 9
2 8 2
4
3 12
= –3 + 6 ×9 + 15 ×(–7) + 20 ×9
The backward difference operator is denoted by the symbol ∇. Let us see some of the first order
backward differences
∇oY = oY – oT
∇o; = o; – oY
………………..
∇ob = ob – obMY
or ∇oª = oª – oªMY ¬ = , – 1, − 2, … , 1
∇; o; = ∇o; – ∇oY
……………
∇; ob = ∇ob – ∇obMY .
x y ∇y ∇; o ∇> y
T oT
∇oY
Y oY ∇; o;
∇o; ∇> o>
; o; ∇; o>
∇o>
> o>
x 1 2 3 4 5 6 7 8
y 1 8 27 64 125 216 343 512
⇒ ∇; o>C = – 0.0039
⇒∇o>T = 0.0775
⇒ 0.500−o;C = 0.0775
⇒ −o;C = 0.4225.
EXAMPLE 4.6: Construct the backward difference table for the data and find ∇2f(0)
-1 0 1 2
() -8 3 1 12
o ∇o ∇; o ∇> y
-1 -8
11
0 3 -13
-2 26
1 1 13
11
2 12
The central difference operator is denoted by the symbol δ. Let us see some of the first order
central differences
δoe = oY − oT
:
δo= = o; − oY ………………(4.7)
:
…………..
δobMe = ob − obMY
:
δ; oª = δoª_ − δoªM
: :
δ; oY = δo= − δoe
: :
δ; o; = δoB − δo=
: :
…..……….
δ; ob = δob_e − δobMe
: :
o δy δ; o δ> y
T oT
δoY/;
Y oY δ; oY
δo>/; δ> o>/;
; o; δ; o;
δoC/;
> o>
EXAMPLE 4.7: Construct a central difference table and obtain from the following data.
x 1 2 3 4 5
y 2 5 10 17 26
Find δ; oY and δoC/; .
x y δy δ; o δ> y δA o
1 2
3
2 5 2
5 0
3 10 2 0
7 0
4 17 2
9
5 26
The divided difference operator is denoted by the symbol ϴ . Let (T ), (Y ), … , (b ) be
( + 1) entries of a function o = (), where () is a polynomial function of degree n
corresponding to the arguments T , Y , ; ,..., b (which need not be equally spaced) . Then the
first order divided differences of () for the arguments T , Y is defined as
ϴyT = ϴf(xT ) =
¸(8e )M¸(8¹ )
8e M8¹
Like wise
⊖; oT = , ⊖; oY =
ϴme Mϴm¹ ϴm: Mϴme
8: M8¹ 8= M8e
and …………..
The third order divided differences for the arguments T , Y , ; , > is defined as
⊖> oT = ⊖> oY =
⊖: re M⊖: r¹ ⊖: r: M⊖: re
8= M8¹ 8@ M8e
, and ………..
o = () ϴy ⊖; o ⊖> o
T oT
re Mr¹
ϴyT =
8e M8¹
Y oY ⊖; oT =
ϴme Mϴm¹
8: M8¹
⊖> oT =
r: Mre ⊖: re M⊖: r¹
ϴyY = 8= M8¹
; o;
8: M8e
⊖; oY =
r= Mr: ϴm: Mϴme
ϴy; = 8= M8e
8= M8:
⊖> oY =
⊖: r: M⊖: re
> o> 8@ M8e
⊖; o; =
r@ Mr= ϴm= Mϴm:
ϴy> = 8@ M8:
8@ M8=
A oA
1 6 10 11
() 3 31 1011 1343
o ϴo ⊖; o ⊖> o
1 3
14
3 31 10
64
6 223 19 1
179
10 1011 27 1
332
11 1343
EXAMPLE 4.9: Obtain the divided difference table for the data
-1 0 2 3
() -8 3 1 12
SOLUTION: We have the following divided difference table for the data.
o ϴo ⊖; o ⊖> o
-1 -8
11
0 3 −4
−1 2
2 1 4
11
3 12
…………………………………………
J MY () = ( − ℎ)
…………..
II. Differential operator: The differential operator for a function o = () is defined by:
D() = ä8 ()
ä
III. Averaging operator: The operator µ is averaging operator and is defined by:
IV. Unit operator: The unit operator 1 has a property that 1. () = (). It is also called
identity operator.
Remark: 1 = J T = ∆T = ∇T = c T = μT = úT
BETWEEN OPERATORS
4.6 RELATION BETW
There are few relations defined between these operators. Some of them are:
1. ∇= 1 – E–1 or E= (1 – ∇)–1
E∇= ∇E= ∆
2. ∆= E– 1 or E= 1 + ∆
3.
4. E= ehD= 1 + ∆, where D is the differential operator.
5. δ= E1/2– E–1/2
Y
;
6. µ= (E1/2+ E–1/2)
1/2
7. δE =∆
PROOF:
4. J() = ( + ℎ)
= () + ℎD() +
({5): ({5)=
;! >!
f(x) + f(x) + …………
({5): ({5)=
;! >!
= (1+ℎD + + +……………)()
= 7 {5 ()
J () = 7 {5 () or E = 7 {5
Since, E = 1 + ∆, therefore ∆= 7 {5 – 1.
=J Y/;
() –J MY/;
(x) = (JY/; – J MY/; ) ()
Therefore, δJY/; = ∆
c. ∆; =δ; +µδ> +
δ@
A
d. ∆ o = ∇ o_b = δb o_f
b b
:
PROOF:
δ: Y ( e/: – e/: ):
; ; ;
a. µδ + = (E1/2+ E–1/2) ( E1/2– E–1/2) +
OMO e OM;_ Oe
; ;
= +
;OM;
;
=
= E-1
=∆
δ: Y
b. 1 + A = 1+ A(E1/2- E–1/2)2
= 1 + A(E1/2− Oe/:)2
Y Y
e
A_( Y/;M e/: );
6
A
=
6e
A_( e/: );
6
A
=
(7e):
A_
7
A
=
AO_ : M; _Y
AO
=
: _; _Y
AO
=
_;_Y/
A
=
Y
A
= (E1/2+E–1/2)2
= µ;
δ:
c. Since µδ + ; = ∆, squaring both sides
∆; = (µδ + ); = µ ; δ; + µδ> +
δ: δ@
; A
:
δ@ δ:
= (1 + )δ; + µδ> + [b/c µ ; = 1 + ]
A A A
δ
=( ) o_b
MY b
=( ) o_f
MY b
e/:
:
= (J − 1) J Mb/; o_f
b
[b/c J Mb/; o_f =o_f ]
: : :
=(J − 1)b o
= ∆b o [b/c E-1=∆ ]
EXAMPLE 4.11: Prove that
a. ∆∇ = ∇∆ = ∆ − ∇ =δ;
b. ∆ + ∇ = −
∆ ∇
∇ ∆
= ∆; E MY ()
]
…………..(i)
∇∆() = (∆E MY )∆()
= ∆; E MY ()
(∆ − ∇)() = (∆ − ∆E MY )()
…………(ii)
= ∆(1 − E MY )()
=∆(E − 1)E MY () [b/c ∆= E − 1
= ∆; E MY ()
]
……………..(iii)
Also
δ; () = (∆E MY/; ); () [b/c δ =∆E MY/;
= ∆; E MY ()
]
……………..(iv)
Therefore,
= (E − E MY )()
O : MY
O
=( )()
[b/c ∆= E − 1 ]
O_Y
O
= (E-1) ( )()
= ∆(1+E MY )()
= ∆(1+E MY )()
= (∆+∆E MY )() [ b/c ∇=∆E MY ]
= (∆+∇)()
Therefore,
∆+∇=∇−∆
∆ ∇
EXAMPLE 4.12: Evaluate the following function with the interval of differencing being h.
1. ∆; (37 8 )
2. ∆; ( ; + )
3. ∆ (¤ MY )
SOLUTION:
1. We have∆(3ex) = 3(7 8_{ − 7 8 )
= 3(7 8 7 { − 7 8 )
= 37 8 (7 { − 1)
Therefore,
∆; (37 8 ) =∆(∆(37 8 )) =∆(37 8 (7 { − 1))
= 3(7 { − 1)[7 8_{ − 7 8 ]
= 37 8 (7 { − 1);
2. ∆ ( ; + ) =( + ℎ); − ; + + ℎ– ]
= ; − 2ℎ + ℎ; − ; + + ℎ −
= ℎ; − 2ℎ + ℎ
∆; (x2+x) = ∆ (∆ (x2+x))
Hence,
= ∆ (ℎ; − 2ℎ + ℎ)
= ℎ; − 2ℎ( + ℎ − ) + ℎ
=¤ MY Y_(8_{)8
(8_{)M8
= ¤ MY
{
Y_8{_8 :
∆:
EXAMPLE 4.13: Show that 7 8 =( ) 7 8 . ∆: ,with interval of differencing being h.
¢£
7 8 =∆; J MY ei = ∆; eiM
∆:
= ei
= 1 − + − +………...]lT
Y ∆ ∆: ∆=
; ; A V
x y ∆y ∆; o
T oT
∆oT
Y oY ∆; oT
∆oY
; o; ∆; oY + c
∆o; + c
> o> +c ∆; o; − 2c
∆o> − c
A oA ∆; o> + c
∆oA
C oC ∆; oA
∆oC
E oE
x 1 3 3 4 5 6 7
y 2 5 10 18 26 37 50
x Y ∆y ∆; o ∆> o
1 2
3
2 5 2
5 1
3 10 3
8 -3
4 18 0
8 3
5 26 3
11 -1
6 37 2
13
7 50
From the seventh properties above the maximum error lies in the two middle terms -3 and 3 and the
incorrect value of y lies between these two middle terms as shown by the horizontal arrow. Here the
sum of all the third differences is zero and the adjacent values –3, 3 are equal in magnitude.
Hence horizontal line between –3 and 3 points out the incorrect functional value 18. Therefore
coefficient of first middle term on the expansion of (1 – p)3= –3
⇒ –3e = –3 ⇒ e= 1
Therefore, Correct functional value = 18 – 1 = 17.
In the above table, the fourth order differences column have algebraic sum of all the values equal to 0.
The middle term of this difference column is [Link],
6e= 6 or e= 1
Hence, correct value is given by 16 – 1 = 15.
EXAMPLE 4.17: Locate the error in following entries and correct it.
1.203, 1.424, 1.681, 1.992, 2.379, 2.848, 3.429, 4.136
SOLUTION: Difference table for given data is as follows:
In eq. (4.11), except missing term, each term is known and hence from this way missing term can
be obtained. If two values of () are missing then in that case only (– 1) values of () can be
given by a polynomial of degree (– 2). That is
∆bMY () = 0 or (E − 1)bMY ( ) = 0
This gives for = T , (the first tabulated value) and for = Y , (second tabulated value) and by
solving these two we get the two missing values for given function (). Similarly method
proceeds to find three and more missing terms in given function ().
EXAMPLE 4.18: Estimate the missing term in the following table:
0 1 2 3 4
o = () 1 3 9 - 81
oT oY o; o> oA
∆A ( ) = 0 ,Ʉ
SOLUTION: Since we have given 4 values, therefore
∆A oT = 0 ,Ʉ
Or let's do in the following way. Since we have given four values of the arguments, we have
∆> oY − ∆> oT = 0
∆; o; − ∆; oY − (∆; oY − ∆; oT ) = 0
∆; o; − 2∆; oY + ∆; oT ) = 0
∆o> − ∆o; − 2(∆o; − ∆oY ) + ∆oY − ∆oT = 0
∆o> − 3∆o; + 3∆oY − ∆oT = 0
oA − o> − 3(o> − o; ) + 3(o; − oY ) − (oY − oT ) = 0
oA − 4o> + 6o; − 4oY + oT = 0
Now substituting the values from the table, we have
81 − 4o> + 6(9) − 4(3) + 1 = 0
81 − 4o> + 54 − 12 + 1 = 0
4o> = 124, hence o> = 31
EXAMPLE 4.19: Find the missing value of the data:
1 2 3 4 5
() 7 - 13 21 37
oT oY o; o> oA
SOLUTION: Since 4 values are known, let us assume the fourth order differences being zero.
Also since one value is unknown, we assume
∆A oT = 0 ,Ʉ
Or let's do in the following way. Since we have given four values of the arguments, we have
∆> oY − ∆> oT = 0
∆; o; − ∆; oY − (∆; oY − ∆; oT ) = 0
∆; o; − 2∆; oY + ∆; oT ) = 0
∆o> − ∆o; − 2(∆o; − ∆oY ) + ∆oY − ∆oT = 0
∆o> − 3∆o; + 3∆oY − ∆oT = 0
oA − o> − 3(o> − o; ) + 3(o; − oY ) − (oY − oT ) = 0
oA − 4o> + 6o; − 4oY + oT = 0
Now substituting the values from the table, we have
37 − 4(21) + 6(13) − 4oY + 7 = 0
37 − 84 + 78 − 4oY + 7 = 0,
4oY = 38,
oY = 9.5.
EXAMPLE 4.20: Find the missing values in the table:
45 50 55 60 56
() 3 - 2 - -2.4
As only three entries y0, y2, y4 are given, the function y can be represented by a second degree
c. ∆; ( 8), ℎ = 1
Y
INTRODUCTION
Interpolation is the technique of estimating the value of a function for any intermediate value of
the independent variable. The process of computing or finding the value of a function for any value
of the independent variable outside the given range is called extrapolation. Here, interpolation
denotes the method of computing the value of the function o = () for any given value of the
independent variable x when a set of values of o = () for certain values of are known or
given. It is required to find a simpler function, sayÆ(), such that () and Æ() agree at the set
of tabulated points. If Æ() is a polynomial, then the process is called Polynomial interpolation
and Æ() is called the interpolating polynomial. Similarly, different types of interpolation arise
depending on whether Æ() is a finite trigonometric series, series of Bessel functions, etc. In this
chapter we shall be concerned with Polynomial interpolation only.
The study of interpolation is based on the assumption that there are no sudden jumps in the values
of the dependent variable for the period under consideration. It is also assumed that the rate of
change of figures from one period to another is uniform.
Numerical interpolation actually approximates the actual curve by numerical polynomials using
the discrete data. Interpolation is based on calculus of finite differences.
If the discrete data is unequally spaced, we use the following numerical interpolation methods
Newton
wton’
5.1.1 Newton’s forw
forward
ard difference Int
Interp
erpolat
olation
Statement: Let o = () be a function which takes the values (T , oT ) (Y , oY ) (; , o; ) ,…, (b ,
ob ) for ( + 1) equally spaced values of such that ª = T + ¬ℎ for ¬ = 0,1,2, … , . One form of
the unique nth-degree polynomial that passes through the n + 1 points is given by
Where l =
8M8¹
{
…………….. (5.1)
Proof: Since we have ( + 1) data points there exist a function polynomial () of nth degree,
given by
T =oT …………….(5.3)
Y = = …………….(5.4)
re Mr¹ ∆r¹
{ {
∆: r¹
; = …………… (5.5)
;!{:
Similarly,
∆= r¹ ∆f r¹
> = ,…., b =
>!{= b!{f
Let = T + lℎ or − T = lℎ
− Y = − T + T − Y = lℎ − ℎ = (l − 1)ℎ
− ; = − T + T − ; = lℎ − 2ℎ = (l − 2)ℎ
……………………………………………………… ….………(5.6)
− b = − T + T − b =lℎ – ℎ = (l − )ℎ
Then substituting the values in eq. (5.6) in to in the above equation, we get
∆: r¹ ∆= r¹ ∆@ r¹ ∆f r¹
() = oT +u∆oT +u(u-1) ;!
+ u(u-1)(u-2) >!
+u(u-1)(u-2)(u-3) A!
+……+u(u-1)(u-2)…(u-n+1) b!
8M8¹
{
Where u=
Which is required Newton's forward difference interpolation formula. This formula is useful
for interpolating the values of () near the starting of the set of data values given.
EXAMPLE 5.1: The values of §¬ are given below for different values of . Find the value of sin
32°.
30T 35T 40T 45T 50T
o = §¬ 0.500 0.5736 0.6428 0.7071 0.7660
SOLUTION: = 32° is very near to the starting value T = 30°. We compute sin 32° by using
Newton’s forward interpolation formula.
Therefore,
l=
8M8¹
= 0.4
{
∆: r¹ ∆= r¹ ∆@ r¹ ∆f r¹
() = oT +u∆oT +u(u-1)
;! >! A! b!
+ u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……+u(u-1)(u-2)…(u-n+1)
= 0.529936
EXAMPLE 5.2: Find the number of students from the following data who secured marks
Therefore,
l=
8M8¹ Y
=
{ ;
∆: r¹ ∆= r¹ ∆@ r¹ ∆f r¹
() = oT +u∆oT +u(u-1)
;! >! A! b!
+ u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……+u(u-1)(u-2)…(u-n+1)
Y Y Y ;; Y Y Y (MC;) Y Y Y Y EA
= 35 + ;(48) +;(; -1) ;! + ;(; -1)(; -2) >!
+;(; -1)(; -2)(; -3) A!
= 59 – 8.5
= 50.5
≅51
Therefore the Number of students who secured not more than 45 marks are 51.
EXAMPLE 5.3: Find the cubic polynomial which takes the following data:
0 1 2 3
() 1 0 1 10
∆: r¹ ∆= r¹ ∆@ r¹ ∆f r¹
() = oT +u∆oT +u(u-1)
;! >! A! b!
+ u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……+u(u-1)(u-2)…(u-n+1)
Now = T +uh which implies u= = , where ℎ = 1, T =0. Hence the formula becomes
8M8¹ 8
{ Y
=
∆: r¹ ∆= r¹
() = oT +u∆oT + u(u-1)
;! >!
+ u(u-1)(u-2)
EXERCISE 5.1
(1) Use Newton's forward difference interpolation formula to find o at = 23 given the
following table of Values:
21 25 29 33 37
o 18.4708 17.8144 17.1070 16.3432 15.5154
(2) From the following table estimate the number of students who obtained marks in computer
programming between 45 and 50.
Statement: Let o = () be a function which takes the values (T , oT ), (Y , oY ), (; , o; ),…, (b ,
ob ) for ( + 1) equally spaced values of such that ª = b + ¬ℎ for ¬ = 0,1,2, … , . One form of
the unique nth-degree polynomial that passes through the n + 1 points is given by
8M8f
{
Where u= ……………..(5.7)
Proof: Since we have ( + 1) data points, there exists a polynomial function () of degree n,
given by
Put = b T = ob ……….(5.9)
Y = =
rfe Mrf ∇rf
{ {
……….(5.10)
∇; oT = ; (2ℎ; )
∇: rf
;!{:
a2 = …… (5.11)
Similarly,
> = ,…., b =
∇= r¹ ∇f rf
>!{= b!{f
()=ob +( – b )∇ob +(– b )( –bMY ) +(– b )( –bMY )( –bM; ) +…+( –b )( –bMY )…( –Y )
∇: rf ∇= rf ∇f rf
;! >! b!
I bMY ! I b 9 b I bMY ! lz 9 z ! l 9 1 z
I bM; ! I b 9 b I bM; ! lz 9 2z ! l 9 2 z
……………………………………………………………………….. …………(5.12)
I T ! I b 9 b I T ! lz 9 z ! l 9 z
Then substituting the values in eq. (5.12) in to the above equation, we get
Which is required Newton's backward difference interpolation formula. This formula is useful
for interpolating the values of near the end of the set of data values given.
EXAMPLE 5.4: Using Newton’s backward difference formula's find the value of 7 MY.U from the
following table of value of 7 M8
SOLUTION: The backward difference table for the given data as follows
∇: rf ∇= rf ∇@ rf
() =ob +u∇ob +u(u+1)
;! >! A!
+u(u+1)(u+2) +u(u+1)(u+2)(u+3) +……….
[Link] [Link]>>
;! >!
=0.1353 +(-0.4)(-0.0385) +(-0.4)(-0.4 +1) +(-0.4)(-0.4+1)(-0.4+2)
= 0.149
EXAMPLE 5.5: The areas A of a circle of diameter d is given for the following values and fined A
for 105.
d 80 85 90 95 100
A 5026 5674 6362 7088 7854
∇: rf ∇= rf ∇@ rf
() =ob +u∇ob +u(u+1)
;! >! A!
+u(u+1)(u+2) +u(u+1)(u+2)(u+3) +……….
= 8666
EXAMPLE 5.6: The following data gives the melting point of an alloy of lead and zinc, where ¤
is the temperature in degrees c and Í is the percentage of lead in the alloy.
Í 40 50 60 70 80 90
¤ 184 204 226 250 276 304
Find the melting point of the alloy containing 84 per cent lead.
SOLUTION: The value of 84 is near the end of the table, therefore we use the Newton’s backward
difference interpolation formula.
¤ Í ∆p ∆; Í ∆> Í ∆A Í
40 184
20
50 204 2
22 0
60 226 2 0
24 0
70 250 2 0
26 0
80 276 2
28
90 304
∇: Pf ∇= Pf ∇@ Pf
Í() =Íb +u∇Íb +u(u+1)
;! >! A!
+u(u+1)(u+2) +u(u+1)(u+2)(u+3) +……….
= 304- 16.8 - 0. 24
=286.96
EXERCISE 5.2
1. Using the Newton’s backward difference formula, find the polynomial () satisfying the
following data. Hence, find (2).
x 0 5 10 15
() 14 379 1444 3584
2. The profits of a company (in thousands of Birr) are given below: Calculate the total profits
between 1999–2002.
Let Æ() denote the interpolating polynomial. Also, let the function () be continuous and
possess continuous derivatives within the interval (T , b ). Now, defining the auxiliary function
;(¤) as
(t−x )(t−x )…(t−x )
;(¤) = (¤) − Æ(¤) − {( ) − Æ()} (x−x 0)(x−x1 )…(x−xn )
0 1 n
……………….. (5.13)
The expression(t − xT )(t − x; ) … (t − x^) is a polynomial of degree (n+ 1) in ¤ and the coefficient
of ¤ = 1. Hence, the (n + 1)th derivative polynomial is (n+ 1)!. That is
(^_Y)!
Hence
Since − T =lℎ or − Y =(l − 1)ℎ, … , − b =(l − )ℎ where ℎ is the interval of differencing,
we have
ú = {∆
Y
We have
úb_Y = ∆b_Y
Y
{f¥e
(b_Y) () =
∆f¥e ¸(8¹ )
b_Y
Or ………….(5.16)
∆f¥e r¹
() = l(l − 1)(l − 2) … (l − )
(b_Y)!
……….(5.17)
(t − x^ )(t − x^MY ) … (t − xT )
; (¤) = (¤) − Æ(¤) − {() − Æ()}
(x − x^ )(x − x^MY ) … (x − xT )
and proceeding as above, we obtain the error in the Newton’s backward interpolation formula as
∇f¥e r¹
() = (b_Y)!
l(l + 1)(l + 2) … (l + )
Where l =
8M8f
{
EXAMPLE 5.7: Using Newton’s forward interpolation formula find the value of ¦R§ 52º from the
following data and estimate the error.
l= = = 1.4
8M8¹ C;MAC
{ C
∆: r¹ ∆= r¹ ∆@ r¹ ∆f r¹
() = oT +u∆oT +u(u-1)
;! >! A! b!
+ u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……+u(u-1)(u-2)…(u-n+1)
−0.0049 0.0005
o = (52) = 0.7071 + 1.4(−0.0643) + 1.4(1.4 − 1) + 1.4(1.4 − 1)(1.4 − 2)
2! 3!
= 0.615680
∆f¥e r¹
Error = (b_Y)!
l(l − 1)(l − 2) … (l − ) when = 2
∆: r¹
= l(l − 1)(l − 2) >!
= 1.4(1.4 − 1)(1.4 − 2)
[Link]
E
= –0.000028
In the previous sections we have seen interpolation formulas that are fundamental and are
applicable to nearly at the beginning and the end of set of data points. But they do not converge as
rapidly as desired if we use them to interpolate near the middle points. So we need methods which
will be used for interpolation at the central. The central difference interpolation formula answers
our question. These formulas are used for interpolation near the middle of an argument values.
Before looking the central difference methods, let's see basic formulas for the derivation of them.
Statement: Suppose o = () is a function which takes the values …(M; , oM; ),
(MY , oMY ), (T , oT ), (Y , oY ) (; , o; ),…, for ( + 1) equally spaced values of such that ª =
T + ¬ℎ for ¬ = 0,1,2, … , . One form of the unique nth-degree polynomial that passes through the
n + 1 points is given by
Proof:
∆: r¹ ∆= r¹ ∆@ r¹ ∆f r¹
() = oT +u∆oT +u(u-1)
;! >! A! b!
+u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……+u(u-1)(u-2)…(un+1)
Let us assume a polynomial equation by using the arrow marks shown in the table below:
oP =oT +?Y ∆oT +?; ∆; yMY +?> ∆> yMY +?A ∆A yM; +……….. …….(5.19)
Where?Y , G; , … , ?b , are coefficients to be determined. From the Newton’s forward interpolation
formula eq.(5.18), we have
∆: r¹ ∆= r¹ ∆@ r¹
oP =oP_T = J @ oT = oT +u∆oT +u(u-1) ;!
+ u(u-1)(u-2) >!
+u(u-1)(u-2)(u-3) A!
+……… …….(5.20)
Similarly,
and so on.
x y ∆ ∆] ∆ ∆A ∆B ∆C
x-3 y-3
∆yM>
x-2 y-2 ∆; yM>
∆yM; ∆> yM>
x-1 y-1 ∆; yM; ∆A yM>
∆yMY ∆> yM; ∆C yM>
x0 y0 ∆; yMY ∆A yM; ∆E yM>
∆yT ∆> yMY ∆C yM;
x1 y1 ∆; yT ∆A yMY
∆yY ∆> yT
x2 y2 ∆; yY
∆y;
x3 y3
oP = oT +?Y ∆yT +?; (∆; yT -∆> yT +∆A yT - ∆C yT +… )+ ?> (∆> yT -∆A yT + ∆C yT -∆E yT + … )
?Y = u
?; =
q(qMY)
;!
?> =
(q_Y)q(qMY)
>!
Here we haveT = 1.4, =1.45 and ℎ=0.2. Then using the formula
8M8¹ [Link].A
{ T.;
u= = = 0.25
(0.256) +
T.;C(T.;CMY) (T.;C_Y)T.;C(T.;CMY)
;! >!
= -0.016 +0.25(0.352)+ (0.048)+….
= 0.047875
SOLUTION: Let us take the origin at T = 29, =30 and ℎ =4. Then using the formula
8M8¹ >TM;U
{ A
u= = = 0.25
(−0.0564 ) + ( −0.0054)+….
T.;C(T.;CMY) (T.;C_Y);C(T.;CMY)
;! >!
= 17.1070 + 0.25(-0.7638) +
= 16.9216
Proof: This formula is also solved by using Newton’s forward difference formula.
∆: r¹ ∆= r¹ ∆@ r¹
o = () =oT +u∆oT + u(u-1)
;! >! A!
+ u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……… ….(5.22)
Let us assume a polynomial equation by using the arrow marks shown in the table below:
oP = yo+?YÇ ∆yMY +?;Ç ∆; yMY +?>Ç ∆> yM; +?AÇ ∆A yM; +……….. …….(5.23)
Where ?YÇ , ?;Ç , ….,?bÇ are coefficients to be determined. From the Newton’s forward interpolation
formula eq. (5.22), we have
∆: r¹ ∆= r¹ ∆@ r¹
oP = J P oT = oT +u∆oT + u(u-1)
;! >! A!
+ u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……… …….(5.24)
Similarly,
and so on.
x y ∆ ∆] ∆ ∆A ∆B ∆C
x-3 y-3
∆yM>
x-2 y-2 ∆; yM>
∆yM; ∆> yM>
x-1 y-1 ∆; yM; ∆A yM>
∆yMY ∆> yM; ∆C yM>
x0 y0 ∆; yMY ∆A yM; ∆E yM>
∆yT ∆> yMY ∆C yM;
x1 y1 ∆; yT ∆A yMY
∆yY ∆> yT
x2 y2 ∆; yY
∆y;
x3 y3
oP =yT +?YÇ (∆yT − ∆; yT + ∆> yT − ∆A yT + ⋯ . )+?;Ç (∆; yT -∆> yT +∆A yT -∆C yT +…)
?YÇ = u
?;Ç =
q(q_Y)
;!
?>Ç =
(q_Y)q(qMY)
>!
?AÇ =
(q_Y)(q_;)q(qMY)
A!
, etc.
EXAMPLE 5.10: Using Gauss backward interpolation formula, find the population for the year
[Link] that
(1) + (−4)
MT.C(MT.C_Y) (MT.C_Y)MT.C([Link])
;! >!
= 39 +(-0.5)(12) +
= 32.625
EXAMPLE 5.11: If () is a polynomial of degree four find the value of (5.8) using Gauss’s
backward formula from the following data
⇒(6)=682+648
⇒(6)=1330
⇒ (E−1)3(4) =132
⇒(7) = 2448
Take T = 6, ℎ = 1, = 5.8
8M8¹ [Link]
{ Y
u= = = -0.2
= 1162.944
Proof: Stirling’s Formula will be obtained by taking the average of Gauss forward difference
formula and Gauss Backward difference formula.
Therefore,
Remark: The method is useful when |l| < ½ RQ − 1/2 < l < 1/2. It gives best estimation
when −1/4 < l < 1/4.
2 4 6 8 10
o -2 1 3 8 20
u o ∆y ∆; o ∆> y ∆A o
-2 2 -2
3
-1 4 1 -1
2 4
0 6 3 3 0
5 4
1 8 8 7
12
2 10 20
(3) + ¼ ½+0
;_C q: (q_Y)q(qMY) A_A
; ;! >! ;
= 3+u( )+
= 3+ u + l ; + (l > − l)
X > ;
; ; >
= 3 + l > + l ;+ l
; > YX
> ; E
l > + ; l ; + E l+3
; > YX
>
=
= ( ) + ( ) + ( )+3
; 8ME > > 8ME ; YX 8ME
> ; ; ; E ;
u o ∆y ∆; o ∆> y ∆A o
-2 20 49225
-909
-1 25 48316 -171
-1080 -59
0 30 47236 -230 -21
-1310 -80
1 35 845926 -310
-1620
2 40 44306
By Stirling’s formula,
= 47691.8256.
∆ yMY +
E(EM) ∆] Fµ _∆] F (EM )E(EM) ∆A F M∆A F]
oP = oT +u∆yT + ¼ ½+ +⋯
] (q_Y)q(qMY)(qM;)
]! ] ! A! ]
Similarly
Now Substituting eq. (5.27) and eq. (5.28) in to eq. (5.26), we get
∆> yMY +
q(qMY) ∆: me ∆: me ∆@ me ∆@ me
oP =oT +u∆yT+ + )+ ( + )+……
(q_Y)q(qMY) (q_Y)q(qMY)(qM;)
;! ; ; >! A! ; ;
∆ yMY +
E(EM) ∆] Fµ _∆] F (EM )E(EM) ∆A F M∆A F]
oP = oT +u∆yT + ¼ ½+ +⋯
] (q_Y)q(qMY)(qM;)
]! ] ! A! ]
Remark: This formula is very useful when l = 1/2 and gives best result when¼ < l < 3/4.
EXAMPLE 5.14: The pressure Í of wind corresponding to velocity I is given by following data.
Estimate pressure whenI = 25.
µ.BM µ.B(µ.BM)
¼ ½+ (−0.4)
T.C([Link]) [Link].Y ]
;! ; !
= 2+0.5(2.4) +
= ½(6.4)+0-0.16250+0
= 3.2-0.16250+0
= 3.0303750
EXAMPLE 5.15: Find the value of (25) from the following data using Bessel’s formula. Data
being (20) = 2854, (24) = 3162, (28) = 3544, (32) = 3992
∆> yMY +
q(qMY) ∆: m¹ ∓me (qM )q(qMY) ∆@ me M∆@ m:
oP = oT +u∆yT + ¼ ½+ +⋯
: (q_Y)q(qMY)(qM;)
;! ; >! A! ;
µ.BM µ.B(µ.BM)
¼ ½+ ( −8)
µ.]B(µ.]BM) MB_CC ]
]! ] !
= 3162 + 0.25(382) +
= 3250.875
EXAMPLE 5.16: Using Laplace Everett’s formula, find (30), if (20) = 2854, (28) =
3162, (36) = 7088, (44) = 7984
l=
8M8¹ >TM;V
{ V
= = 0.25
By Everett’s formula,
= 4064
Hence = 4064.
l=
8M8¹ ;>.XCM;>
{ Y
= = 0.75
By Everett’s formula,
(−0.0008) + (−0.0001)+….}
µ.]BÀ(µ.]B)] MÁ µ.]BÀ(µ.]B)] MÁÀ(µ.]B)] MAÁ
! B!
={ 0.25(1.3617+
= 1.375874
Now
EXERCISE 5.3
4. Apply Everett’s formula for the following data, to find the value of (31) and (32),if
(20) = 3010, (25) = 3979, (30) = 4771, (35) = 5441, (40) = 6021, (45) =
6532.
Statement:
Let (T ), (Y ), … , (b ) be (n+ 1) entries of a function o = (), corresponding to the
arguments T , Y ,..., b which may not be equally spaced. One form of the unique nth-degree
polynomial that passes through the n + 1 points is given by
¾(i¹)
For = T ,(T ) = T ( − Y )( − ; )…( − b )⇒ T =
(i¹ Mie )(i¹ Mi: )…(i¹ Mik )
..…(5.31)
¾(ie )
For = Y ,(Y ) = Y ( − T )( − ; )…( − b )⇒ Y =
(ie Mi¹ )(ie Mi:)…(ieMik )
……(5.32)
Similarly
This is called Lagrange’s interpolation formula. In eq. (5.34), dividing both sides by
¾(i: ) Y ¾(ik ) Y
(i: Mi¹ )(i: MY)…(i: Mik ) (iMi: ) (ik Mi¹ )(ik Mie )…(ik Mike ) (iMik )
+ +…+
Error of interpolation
We assume that () has continuous derivatives of order up to n + 1 for all ∈(a, b). Since, ()
is approximated by oP (), the results contain errors. We define the error of interpolation or
truncation error as
Without giving the derivation, we write the expression for the error of interpolation as
Since, ∈ is an unknown, it is difficult to find the value of the error. However, we can find a bound
of the error. The bound of the error is obtained a
Note that in eq. (5.37), we compute the maximum absolute value of (x − xT )(x − xY ) … (x − x^ ),
that is max |(x − xT )(x − xY ) … (x − x^ )|and not the maximum of (x − xT )(x − xY ) … (x − x^ ).
EXAMPLE 5.18: Using the data §¬(0.1) = 0.09983 and §¬ (0.2) = 0.19867, find an
approximate value of §¬ (0.15) by Lagrange interpolation. Obtain a bound on the error at
= 0.15.
SOLUTION: We have two data values. The Lagrange linear polynomial is given by
() =
(iMie) (iMi¹ )
(0 ) + (1 )
(i¹ Mie ) (ie Mi¹)
(0.15−0.2) ([Link].Y)
= 0.09983+ 0.19867
(0.1−0.2) (T.;MT.Y)
(0.15−0.1)(0.15−0.2)
. J = (−sin (∈)) = 0.00125 sin (∈)
;!
EXAMPLE 5.19: Use Lagrange’s formula, to find (2) that takes the values
0 1 3
o 0 1 0
EXAMPLE 5.20: Find the unique polynomial O() of degree 2 such that (1) = 1 , (3) =
() = 8 2– 19 + 12
(>CM>T)(>CMAT)(>CMCT) (>CM;C)(>CMAT)(>CMCT)
(;CM>T)(;CMAT)(;CMCT) (>TM;C)(>TMAT)(>TMCT)
= (52) + (67.3)
(>CM;C)(>CM>T)(>CMCT) (>CM;C)(>CM>T)(>CMAT)
(ATM;C)(ATM>T)(ATMCT) (CTM;C)(CTM>T)(CTMAT)
+ (84.1) +…+ (94.4)
The only difference of this formula from Lagrange’s method is that is assumed to be expressible
as a polynomial in o. So on interchanging and o in the Lagrange’s formula we have,
(mMm¹)(mMm:)…(mMmk )
= (m T + 1
(mMme )(mMm: )…(mMmk)
¹ Mme )(m¹ Mm: )…(m¹ Mmk ) (me Mm¹ )(meMm: )…(meMmk )
EXAMPLE 5.22: Find the value of correct to one decimal place for which o = 7, given
1 3 4
o 4 12 19
(mMm¹)(mMm:)…(mMmk )
= T + Y
(mMme )(mMm: )…(mMmk )
(m¹ Mme )(m¹ Mm: )…(m¹ Mmk ) (me Mm¹ )(me Mm: )…(me Mmk )
(XMA)(XMYU) (XMA)(XMY;)
= (AMY;)(AMYU) (1) +
(XMY;)(XMYU)
(Y;MA)(Y;MYU) (YUMA)(YUMY;)
(3) + (4)
EXAMPLE 5.23: Tabulate o = 3 for = 2, 3, 4, 5 and calculate the cube root of 10 correct to
three decimal places.
(YTMV)(YTMEA) (YTMV)(YTM;X)
(2) +
(YTM;X)(YTMEA)
(VM;X)(VMEA) (;XMV)(;XMEA) (EAMV)(EAM;X)
= (3) + (4)
=1.726+0.501−0.07144
= 2.15556
Remark: The main advantage of Lagrange's Interpolation formulas is, it can also be used in case
of equal intervals .
EXERCIS 5.4
1. Using Lagrange’s interpolation formula, find y(10) from the following table
x 5 6 9 11
y 12 13 14 16
2. If o(1) =– 3, o(3) = 9, o(4) = 30 and o(6) = 132, find the four point Lagrange’s
interpolation polynomial which takes the same values as the function y at the given points.
3. Given that (10) = 1754, (15) = 2648, (20) = 3564, find the value of for () =
3000 by using, iterative method of inverse interpolation.
Proof: We drive the Newton's divided difference interpolation formula by induction. Let we have
the following data.
T
o = () oT o
Then
rMr¹
ϴoT =
8M8¹
Y T
o = () oY oT o
Then
ϴyM ϴo0
⊖; oT =
8M8e
⊖2 oM⊖2 o0
⊖> oT =
8M8:
Also
This gives
EXAMPLE 5.24: Apply Newton’s divided difference formula to find the value of (8) if
o ⊖y ⊖; o ⊖> y ⊖A o
1 3
14
3 31 10
64 1
6 223 19 0
197 1
10 1011 17
332
11 1343
= 3 + 98 + 350 + 70
= 521.
EXAMPLE 5.25: Calculate the value of (6) from the following data
1 2 7 8
() 1 5 5 4
o ⊖y ⊖; o ⊖> y
1 1
4
2 5 -2/3
0 1/14
7 5 -1/6
-1
8 4
−
MAT YT
> X
= 1+20+
= 6.2381
0 1 2 4 5 6
() 1 14 15 5 6 19
SOLUTION: To use the Newton divided difference interpolation formula, we first construct the
divided difference table:
o ⊖y ⊖; o ⊖> y ⊖A o ⊖C o
0 1
-13
1 14 -6
1 1
2 15 -2 0
0
-5 1
4 5 2 0
1 1
5 6 6
13
6 19
=10
Remark: The main advantage of Newton's divided difference interpolation formulas is, it can
also be used in case of equal intervals.
1. Find when () = 14 for the following data using Lagrange’s inverse interpolation formula
0 5 10 15
() 16.35 14.88 13.59 12.46
2. By means of Newton’s divided difference formula, Find the value of (8) and (15) from the
following table:
4 5 7 10 11 13
5 7 11 13 17
() 150 392 1452 2366 5202
2. The following table give the marks secured by 100 students in Mathematics:
Range of marks() 30-40 40-50 50-50 60-70 70-80
No. of students 25 35 22 11 7
4 6 8 10
o 1 3 8 10
-1 2 4 5
() -5 13 255 625
6. Using Newton’s divided difference formula, determine f(3) for the data
0 1 2 4 5
() 1 14 15 5 6
0 4 8 12 16
o 15 25 34 37 42
9. Apply Gauss’s backward interpolation formula to find the population of a city in 1986 based
on the following data:
10. Apply Bessel’s interpolation formula to obtain (25), given that (20) = 515, (24) =
438, (28) = 348 è (32) = 249.
11. Use Stirling’s formula to find (35), data being: (20) = 512, (30) = 439, (40) = 346,
and f(5) = 243.
12. Apply Everett’s formula to find the value of of (31), (32), if (20) = 3010, (25) =
3979, (30) = 4771, (35) = 5441, (40) = 6021, (45) = 6532.
INTRODUCTION
The method of obtaining the derivatives of a function using a numerical technique is known as
numerical differentiation.
Approximation to the derivatives can be obtained numerically using the following two
approaches
The choice of the formula is the same as discussed for interpolation if the derivative at a point
near the beginning of a set of values given by a table is required then we use Newton forward
formula, and if the same is required at a point near the end of the set of given tabula values, then
we use Newton’s backward interpolation formula.
If the values of are not equally spaced, we use Newton’s divided difference interpolation
formula or Lagrange’s interpolation formula to get the required value of the derivative.
∆: r¹ ∆= r¹ ∆@ r¹ ∆f r¹
() = oT +u∆oT +u(u-1)
;! >! A! b!
+u(u-1)(u-2) +u(u-1)(u-2)(u-3) +……+u(u-1)(u-2)…(u-n+1)
Where l =
8M8¹
{
…….(6.1)
är ∆: r¹ ∆= r¹ ∆@ r¹
=∆oT +(2u-1)
äq ;! >! A!
+(3u2-6u+2) + (4u3-18u2+22u-6) +…… …...(6.2)
är är äq Y är äq Y
ä8 äq ä8 { äq ä8 {
Now = = where =
Therefore
∆: r¹ ∆= r¹ ∆@ r¹
= ∆yo+ (2u-1)
är Y
ä8 { ;! >! A!
+( 3u2-6u+2) +(4u3-18u2+22u-6) +…… ] …...(6.3)
For tabular values of , we can easily find the derivative by substituting = T .i.e u=0 in eq.
(6.3), we get
∆: r¹ ∆= r¹ ∆@ r¹
= ∆yo− −
är Y
ä8 { ; > A!
+ +….]
= ( ) ( )
ä: r ä är äq Y ä är
ä8 : äq ä8 ä8 { äq ä8
=
o ∆y ∆; o ∆> y ∆A o ∆C o
1 2.7183
1.2
0.6018
3.3201 0.1333
1.4
0.7351 0.0294
4.0552 0.1627 0.0067
1.6
0.8978 0.0361 0.0013
4.9530 0.1988 0.0080
1.0966 0.0441 0.0014
1.8 6.0496 0.2429 0.0094
1.3395 0.0535
2.0 7.3891 0.2964
1.6359
2.2 9.0250
The given point = 1.2 is a tabular value and found at the beginning, we use the derivative of
Newton's forward interpolation formula, we have
∆: r¹ ∆= r¹ ∆@ r¹
= ∆yo− −
är Y
ä8 { ; > A
+ +….]
0.7351 − −
Y [Link];X T.T>EY [Link] [Link]
T.; ; > A C
= + + ….]
= 3.32015
0.1627 − 0.0361 +
Y YY
(T.;): Y;
= (0.0080 ) ]
= 3.34825
o ∆y ∆; o ∆> y ∆A o
1 0
0.1280
1.2 0.1280 0.288
0.4160 0.048
1.4 0.5449 0.336 0
0.7520 0.048
1.6 1.2960 0.384 0
1.1360 0.048
1.8 2.4320 0.432
1.5680
2.0 4.000
Here we have to find the derivatives at = 1.1 which lies between given arguments 1.0 and 1.2.
So apply derivative of Newton’s forward formula, we have
8M8¹ 8MY
{
u= = u = T.; = 5(-1) = 5(1.1-1) = 0.5
∆: r¹ ∆= r¹ ∆@ r¹
= ∆yo+ (2u-1)
är Y
ä8 { ;! >! A!
+( 3u2-6u+2) +(4u3-18u2+22u-6) +…… ]
0.288 +(0.5-1)0.048+ ]
Y
(T.;):
=
∆: r¹ ∆= r¹ ∆@ r¹
= ∆yo− −
är Y
ä8 { ; > A!
+ +….]
0.4160 −
Y T.>>E [Link]
T.; ; >
= + -….]
0.336 − 0.048 ]
Y
(T.;):
=
0.288] = 7.2
Y
[Link]
=
EXAMPLE 3: The table given below reveals the velocity v of a body during the time t. Find its,
acceleration at ¤ = 1.1.
¤ I ∆v ∆; I ∆> v ∆A I
1 43.1
4.6
1.1 47.7 -0.2
4.4 0.1
1.2 52.1 -0.1 0.1
4.3 0.2
1.3 56.4 -0.1
4.4
1.4 60.8
The given point is a tabular value, we have, ¤T =1.1, IT =4.77 and ℎ = 0.1
= ∆vo− −
äJ Y ∆: J¹ ∆= J¹ ∆@ J¹
ä { ; > A!
+ +….]
= T.; 4.4–
Y (MT.Y ) T.;
;
+>]
EXERCISE 6.1
0 1 2 3 4
o 6.9897 7.4036 7.7815 8.1291 8.4510
2. A particle is moving along a straight line. The displacement at some time instance tare given
below. Find the velocity and acceleration of the particle at ¤ = 4.
¤ 1 3 5 7 9 11
0.1405 0.7676 3.5135 9.9351 21.5892 40.0324
Given (n+1) equally spaced data points (T , oT ), (Y , oY ), … (b , ob ).The Newton’s backward
difference interpolation formula is given by
∇: rf ∇= rf ∇@ rf ∇f rf
() = ob +u∇ob + u(u+1) ;!
+ u(u+1)(u+2) >!
+u(u+1)(u+2)(u+3) A!
+…+u(u+1)…(u+n-1) b!
8M8f
{
Where u= ………..(6.6)
= [ ∇yn + (2u+1)
är Y ∇: rf ∇= rf ∇@ rf
ä8 { ;! >! A!
+ (3u2+6u+2) +(4u3+18u2+22u+6) +……] ……...(6.7)
For tabular values of , we can easily find the derivative by substituting = T .i.e. l = 0 in eq.
(6.7), we get
= [ ∇yn +
är Y ∇: rf ∇= rf ∇@ rf
ä8 { ; > A
+ + +……] ……….(6.8)
∇= rf ∇@ rf
[∇; ob + (6u+6)
ä: r Y
ä8 : {: >! ;A
= +(12u2+36u+22) +……]
For tabular values of , we can easily find the derivative by substituting = b .i.e l = 0 in the
above formula, we get
[∇; ob + ∇> ob + 11
ä: r Y ∇@ rf
ä8 : {: Y;
= +……] ………..(6.10)
EXAMPLE 6.4: Find the first and second derivative value at = 2.03 and = 2.02
Since = 2.03 is not found in the table and the value is found at the end of the data. Sothat we
use the Newton's backward interpolation formula
= [ ∇yn + (2u+1)
är Y ∇: rf ∇= rf ∇@ rf
ä8 { ;! >! A!
+ (3u2+6u+2) +(4u3+18u2+22u+6) +……]
Y ([Link];) ([Link];)
T.T; ;! >!
= [−0.0090+ [2(-0.5) +1] + [3(-0.5)2+6(-0.5)+2]
([Link])
A!
+ (4(-0.5)3+18(-0.5)2+22(-0.5) + 6) ]
= -0.44875
=-1.0415
And since = 2.02, b = 2.02 is in the table ,we use the formula
= [ ∇yn +
är Y ∇: rf ∇= rf ∇@ rf
ä8 { ; > A
+ + +……]
Y
T.T;
= [-0.0088 +0.5(0)+0.333(0.0002)]
= -0.43667
[∇; ob + ∇> ob + 11
ä: r Y ∇@ rf
ä8 : {: Y;
= +……]
EXAMPLE 6.5: Find the values of o′(4) and o′′(4) from the following table:
o
0 1 2 3 4
5 8 12 17 26
o ∆y ∆; o ∆> y ∆A o
0 5
1
3
8 1
2
4 0
12 1 3
3
5 3
17 4
9
4 26
Since = 4 is found in the table and the value is found at the end of the data, so that we use the
derivative of Newton's backward formula
= [ ∇yn +
är Y ∇: rf ∇= rf ∇@ rf
ä8 { ; > A
+ + +……]
= [9 +0.5(4)+0.333(3)+(0.25)(3)]
= 12.749
: [∇ ob + ∇ ob + 11
ä: r Y ; > ∇@ rf
ä8 : = { Y;
+……]
= [4 + 3 + Y; (3)] = 9.75
YY
o ∇y ∇; o ∇> y ∇A o
30T 0.500
35
0.0736
T 0.5736 -0.0044
40
0.0692 -0.005
T 0.6428 -0.0049 0
45
0.0643 -0.005
T 0.7071 -0.0054
50
0.0589
T 0.7660
Since = 48 is not found in the table and the value is found at the end of the data and we use
the derivative of Newton's backward formula.
= [ ∇yn + (2u+1)
∇: rf ∇= rf ∇@ rf
¦R§ 48° =
är Y
ä8 { ;! >! A!
+ (3u2+6u+2) +(4u3+18u2+22u+6) +……]
Y ([Link] ) [Link]
= YT {0.0589+ [2(-0.2)+1] ;!
+ [3(-0.2)2+6(-0.2)+2] >!
}
Y
YT
= {0.0589+ 0.000162-0.0015994}
= 0.0060661
EXAMPLE 6.7: The following data gives the velocity of a particle for 10 seconds at an interval of
2 seconds. Find the acceleration when the time ¤=8sec.
¤(sec) 0 2 4 6 8 10
I(m/sec) 4 6 16 34 60 94
¤ I ∇v ∇; I ∇> v
0 4
2
2
6 8
4
10 0
16 8
6
18 0
34 8
26 0
8 60 8
34
10 94
Since ¤ = 8 is found in the table, l = 0. And is found at the end of the data and we use the
derivative of Newton's backward difference interpolation formula . The first derivative at ¤ = 8
= { [ ∇Ib +
äJ Y ∇: Jf ∇= Jf ∇@ Jf
ä ; > A
+ + +…]
= [26 + (;)
Y V
;
= 15
EXERCISE 6.2
1. Find the values of o′(4) and o′′(4) from the following table:
x 0 1 2 3 4
5 8 12 17 26
y
2. The elevations above a datum line of seven points of roads 300 units apart are 135, 149, 157,
183, 201, 205, 193 unit. Find the gradient of the road at the end point.
Given (n+1) equally spaced data points (T , oT ), (Y , oY ), … (b , ob ).If we want to determine the
values of the derivatives of the function near the middle of the given set of arguments. We may
apply any central difference formula. Therefore using Stirling’s formula,
where u= ………….(6.11)
8M8¹
{
∆A yM; +………..
∆me _∆m¹ (>q: MY) ∆= m: _∆= me
) + u∆; yMY + ¼ ½+
är (Aq= M;q)
äq ; >! ; A!
=(
Now
är är äq Y är äq Y
ä8 äq ä8 { äq ä8 {
= = where =
= [(
är Y
∆A yM; +…]
∆me _∆m¹ >q: MY ∆= m: _∆= me
) + u∆; yMY + ¼ ½+
Aq= M;q
ä8 { ; >! ; A!
……..(6.12)
= [(
är Y ∆me _∆m¹ Y ∆= m: _∆= me
) − ¼ ½+…]
ä8 { ; E ;
…….(6.13)
∆A yM; +….
∆= m: _∆= me
= {:[∆; yMY + u ¼ ½+
ä: r Y Y;q: M;
ä8 : ; A!
………(6.14)
At the tabular values, = T , i.e u =0. Therefore, putting l = 0 in eq. (6.14), we get
EXAMPLE 6.8: Find ′(0.6) è ′′(0.6) from the following table:
SOLUTION: Here, the derivatives are required at the central point x= 0.6, so we use Stirling’s
SOLUTION:
= (
är Y ∆me _∆m¹ Y ∆= m: _∆= me
) − ¼ ½+…]
ä8 { ; E ;
(
Y
= ) −E¼ ½+…]
T.;AEV_ T.;V>> Y [Link]>C_ [Link]>X
T.Y ; ;
= 10[0.26505- 0.0006]
’(0.6) = 2.6445
= 0.0365 − (0.0002)+….
Y Y
{: Y;
=1000.0365- 0.000016
= 3.6484
x 60 75 90 105 120
f(x) 28.2 38.2 43.2 40.9 37.7
l o ∇y ∇; o ∇> y
-2 60 28.2
10
-1 75 38.2 -5
5 -2.3
u= = = =0.2
8M8¹ U>MUT Y
{ YC C
Since =93 is not found in the table and the value is found at the middle of the data, so that we
use the Striling's formula. Here ℎ = 0.1, = 0.6
= (
är Y
∆A yM; +…]
∆me _∆m¹ >q: MY ∆= m: _∆= me
) + u∆; yMY + ¼ ½+
Aq= M;q
ä8 { ; >! ; A!
(
Y
= ) + (0.2)(−7.3 )+ ¼ ½+ (8.7)+…]
C_;.> >(T.;): MY M;.>_E.A A(T.;)= M;(T.;)
YC ; >! ; A!
Y [Link]>
YC ;
= [1.35- 1.46- -0.1334]
=-0.03627
1. Compute the values of o′(3) and o′′(3) from the following table:
x 1 2 3 4 5
y 0 1.4 1.65 5.673 8.0978
2. Compute the values of o′(3) and o′′(3) from the following table:
1 2 3 4 5
3. The table given below reveals the velocity ‘I’ of a body during the time ‘¤’ specified. Find the
acceleration at ¤ = 1.3.
Suppose we have (n+1) equally spaced data points (T , oT ), (Y , oY ), … (b , ob ).The Newton’s
divided difference formula is given by
EXAMPLE 6.10: Find ′(5) and ′′(5) from the following table:
1 2 4 8 10
() 0 1 5 21 27
SOLUTION: Here the data are unequally spaced. So we use Newton’s divide difference formula.
o ⊖y ⊖; o ⊖> y ⊖A y
1 0
2
1
1 1/3
4
2 0
5 1/3 -1/144
8
4 -1/16
21 -1/6
3
10 27
Using the the first derivative formula 0f Newton's divided difference, we have
′(5) = 1 +(10- 1- 2) + [(5 -2)(5- 4) +((5- 1)(5 -4)+ (5 -1)(5 -2)](0)+[(5 -2)(5- 4)(5- 8)
Y
>
MY
+ (5 -1)(5- 4)(5- 8)- (5- 1)(5- 2)(5 -8)+ (5 -1)(5- 2)(5- 4)]YAA
2 4 9 13 16 21 29
() 57 1345 66340 492052 1118209 4287844 21242820
SOLUTION: Here, the arguments are not equally spaced and therefore we shall apply Newton’s
divided difference formula. The divided difference table:
o = () ⊖y ⊖; o ⊖> o ⊖A o ⊖C o ⊖E o
2 57
644
4 1345 1765
12999 556
6 66340 7881 45
83928 1186 1
13 402052 22113 64 0
238719 2274 1
16 1118209 49401 89
633927 4054
21 4287844 114265
2119372
29 21242820
+( − 2)( − 4)( − 9)( − 13)(45) + ( − 2)( − 4)( − 9)( − 13)( − 16)(1)
+45( A -28 > + 257 ; − 878x + 936) + C −44 A +705 > − 4990 ; −14984x-14976
Then,
= 644 +1765(2x- 6) 556(3 ; − 30x+ 62) +45(4 > − 84 ; + 514- 878)
är
ä8
= 3530+ 556(6x- 30) +45(12 ; -168x+ 514) +20 > − 528 ; + 4230x- 9980
ä: r
ä8 :
= 60 ; + 24 + 6
When = 5;
=1626
EXAMPLE 6.12: Find the first and second derivatives at = 1.6, for the function represented by
o
1 1.5 2.0 3.0
0.000 0.40547 0.69315 1.09861
SOLUTION: The data is not equally paced. We use the divided difference formulas to find the
1 0.000
0.81094
1.5 0.40547 -0.23558
0.57536 0.061157
2.0 0.69315 -0.113267
0.40546
3 1.09861
′(1.6) = 0.81094 + [(1.6 – 1.0) + (1.6 – 1.5)] (– 0.23558) + [(1.6 – 1.5)(1.6 – 2.0)
EXERCISE 6.4
2. Using Newton’s Divided Difference Formula, find ′(10) from the following data
3 5 11 27 34
o -13 23 99 17315 35606
Time(sec) 0 5 10 15 20
Velocity(m/sec) 0 3 14 69 228
0 1 2 3 4
o = () 1 2.718 7.381 20.086 54.598
1 1.5 2 2.5 3
o = () 0 0.5 1.1 3.2 5.3
2 4 9 10
5. The heat flux R is the quantity of heat following through a unit are of a material per unit time.
It can be computed with Fourier's law
R = − äS
ä
Where T = temperature, X= distance and J has unit S/¯2 .The following temperatures are
measured from the surface (6 = 0) into a stone wall:
6(¬ ¯) 2 4 9 10
INTRODUICTION
Analytical methods will not always give the solution of a definite integral when the function is
difficult to integrate. For instance the following integrals are impossible to integrate analytically
b
I = ∫ f (x )dx …………….. (7.1)
a
The definite integral of a function of a single variable, (), between two limits and can be
viewed as the area under the curve defined by the function.
Sometimes, the function to be integrated is given at discrete data points, and the area under the
curve is needed to be approximated. Here, we will discuss the integration of the region S shaded on
the figure below.
Let o = () be a function, where o takes the values oT , oY , o; , …, ob for x=T , Y , ; ,….., b
I = U? ()è
Let the interval of integration (, ) be divided into n equal subintervals of width
ℎ=(
M?
b
)
Where l = implies dl = dx
8M8¹ Y
{ {
Therefore
This is called Newton-Cotes general quadrature formula. From the general formula eq. (7.3),
we can derive or deduce different integration formula by substituting =1, 2, 3, …, etc.
The trapezoidal approximation of a definite integral is to say that on each subinterval you
approximate by a first-degree polynomial. Putting = 1 in eq. (7.3) and taking the curve
o = () through (T ,oT ) and (Y, oY )so that we have only one interval (T , Y). i.e = 1 with
the two values oT and oY . Consequently we can get only the first differences ∆o . So that
differences of order higher than one vanishs and we get
= ; [ 2oT +(oY − oT )]
{
= [ oT +oY ]
{
;
U8 ()è= [ oY +o; ]
8¹ _;{ {
¹¥ ;
Let o = () be a continuous function with continuous derivatives in the interval [T , b ].
Expanding o in a Taylor’s series around = T , we get
8e 8e
− T ; ÇÇ
V oè = V oT + ( − T )oT Ç + oT + ⋯ è
8¹ 8¹ 2!
Ç ÇÇ ÇÇÇ
oT + oT + oT +……..
{: {= {@
= hoT +
;! E ;A
……..(7.4)
ÇÇ
Likewise, (oT + oY ) = (oT + o(T + ℎ))= (oT +y+ho Ç + oT +……
{ { { {:
; ; ; ;
Ç ÇÇ ÇÇÇ
o + oT + Y; oT +……..
{: {= {@
= hoT + ;! T A
……….(7.5)
Hence, error 7Y in (T , Y ) is obtained from eq. (7.4) and eq. (7.5) as
and so on.
ℎ>
|(, )| ≤ |o ÇÇ + oY ÇÇ + ⋯ . . +ob ÇÇ |
12 T
ℎ>
≤ |o ÇÇ | + |oY ÇÇ | + ⋯ + |ob ÇÇ |
12 T
ℎ> ( − )ℎ;
≤ W = W;
12 ; 12
Now, since ℎ = ( − )/ , the total error in the evaluation of the integral of eq. (7.1) by the
trapezoidal rule is of the order of h2. Therefore, as we increase the number of intervals, the error
decreases.
SOLUTION: let us take h=1/4 .then the values of () = Y_8 : at =0, A,;,A,1 are given below
Y YY>
0 1 3 1
4 4 2
1
By trapezoidal rule
=0.825
¤ MY YT =
A
=0.785398
W; = max L L = 2.
;(>8 @ _;8 : MY)
0 ≤x ≤1 (Y_8 : )@
Y
1 1
ℎ= ; |error| ≤ YE
(2) = = 0.010417
4 12 96
SOLUTION: let us take h=1/4 .then the values of () = √1 − ; at =0, , , ,1 are given
YY>
A;A
below
0 1 3 1
4 4 2
1
By trapezoidal rule
T.;C
;
= [1+2(0.0605+ 0.04133+ 0.2165) + 0] = 0.144895625
Simple Command window on Mat lab to find the Area of f(x) = x2 on [0, 2] with h=0.5
>> y=x.^2;
>> trapz(x,y)
ans =
2.7500
function I = trapzoid(f,a,b,n)
h=(b-a)/n;
x=linspace(a,b,n+1);
fx=feval(f,x);
I=(h/2)*(fx(1)+2*sum(fx(2:1:n))+fx(n+1));
Sample output
>>f = @(x) 1./(1+x.^2);
>> I = trapzoid(f,-5,5,2^4)
I=
2.7468
>> g = @(x) x.^2;
>> I = trapzoid(g,0,2,2^4)
I=
2.6667
Putting = 2 in eq. (7.3) and taking the curve through (T , oT),(Y , oY ) and (; , o; ). Consequently
we can get only up to second order differences. So that differences of order higher than two
vanish, we get
= [ oT +4oY +o; ]
{
>
Similarly,
…………………………………………………..
Remark: To Use Simpson’s 1/3 rule, the given interval of integration must be divided into an
In addition, we have
Hence, from eq. (7.8) and eq. (7.9), the error in the subinterval (T , ; ) is given by
7; = o; ªJ
M{B
UT
7> = o ªJ
M{B
UT A
…………….(7.11)
and so on.
=∑ 7b = + o; ªJ + oA ªJ + ⋯ + o;bM; ªJ ]
M{B
[oT ªJ
UT
ℎC
≤ (o ªJ ( + (oY ªJ ( + ⋯ + (ob ªJ (
90 T
≤ WA = WA
b{B (M?){@
UT YVT
……………(7.12)
()
Mb{B ªJ
YVT
e=
Y YY>
SOLUTION: let us take h=1/4 , then the values of f(x) = Y_8 : at x=0,A,;,A, 1 are given below
By Simpsons'1/3 rule
=0.785
|error| ≤ WA
(M?){@
YVT
WA = max L L = 12.
Y;A8 @ MYCE8: _Y;
0 ≤x ≤1 (Y_8 : )B
Y
1 1
ℎ= ; |error| ≤ ;CE
(12) = = 0.000427
4 180 2340
Remark: From the above example we can conclude that Simpson's 1/3 rule is rapidly
/ 3//2 2/ 5//2 3/
() −1 1 −1
/ 2/ 3/
0 0
oY o>
oT o; oA
Since there are four sub intervals, which is a multiple of two, we can apply Simpson's 1/3 rule
> YZ[ (8)
è => [oT +ob +4{oY +o> + oC +… . + obMY }+ 2{o; +oA + oE +… . + obM; }]
{
I=U 8
MY Y MY
= E [ +4( 0+0)+2(;)+ > ]
≅0.05556
MY Y MY Y
= E[
+ +
>
]=
YV
EXAMPLE 7.5: The velocity of a particle which starts from rest is given by the following table
t(sec) 0 2 4 6 8 10 12 14 16 18 20
v(ft/sec 0 16 29 40 46 51 32 18 8 3 0
S = UT I è¤
;T
\ = > [(0) + 4{(2) + (6) + (10) + (14) + (18)} + 2{(4) + (8) + (12) + (16)} + (20)
{
;
>
= [0 + 4{16 + 40 + 51 + 18 + 3} + 2{29 + 46 + 32 + 8} + 0] = 494.667 feet.
EXAMPLE 7.6: Compute the error in the evaluation of UA R« è by Simpson's 1/3 rule.
C.;
ªJ (x) =
ME
8@
Hence the maximum value of ªJ ()=(A)@. Let the points be taken at width of 0.2.i.e ℎ = 0.2,
ME
then
ªJ ()
M(M?){@
YVT
E=
(( )4 )
M(C.;MA)(T.;)@ −6
YVT 5.2
= = 0.00000025
function In = simpson(f,a,b,n)
% n must be a positive even integer
h=(b-a)/n; %step size
x=linspace(a,b,n+1); %(n+1) equally spaced points
fx=feval(f,x); %function evaluations
%Composite Simpson
In=(h/3)*(fx(1)+4*sum(fx(2:2:n))+2*sum(fx(3:2:n-1))+fx(n+1));
Sample output
>> f = @(x) 1./(1+x.^2);
>> In = simpson(f,-5,5,2^8)
Putting = 3 in eq. (7.3) and taking the curve through (T , oT ),(Y , oY ) ,(; , o;) and (> , o>).
Consequently we can get only up to third order differences. So that differences of order higher
than three vanish. Thus we get
Similarly,
……………………………………………………………
Remark:
(i) To Use Simpson’s 3/8 rule, the given interval of integration must be divided into sub-
intervals whose number n is a multiple of 3.
given below
Here the number of sub interval with h=1/6 is six, which is even and a multiple of three, so that
we can use all the above methods.
By Trapezoidal rule
=
Y
Y;
[1.5+2{3.95543}]=0.784233
=
Y
YV
[1.5+4{2.3631}+2{1.59}]= 0.785389
= 0.7853875
Remark: From the three rules of numerical integration, Simpson's 1/3 rule is more accurate and
it is preferable for numerical integration than the others.
1 4/3 5/3 2
Since the number of subinterval is a multiple of three, we can apply Simpson's 3/8 rule.
Now
Y
YE
= [0.125 + 3 {0.11765 + 0.11111 + 0.10000 + 0.09524} + 2(0.10526) + 0.09091]
= 0.10615.
Therefore Simpson's 3/8 rule give solution is the same with the analytical solution up to five
decimal places
Here we assume that within any four consecutive subintervals of width ℎ each, the interpolating
polynomial approximating () will be of degree 4. Substituting = 4 in the General quadrature
formula given by eq. (7.3) and neglecting all differences above, we get
Likewise
and soon.
REMARK:
1. It should be noted here that the number of subintervals should be taken as a multiple of 4.
ℎX o Jª (Ẋ)
MV
UAC
2. The leading term in the error of formula can be shown as
SOLUTION: Five functional evaluations mean using four sub intervals (i.e. n= 4).
oT oY o; o> oA
Y.;MT YT
A >
h= = , and applying Boole’s rule, we have
= 2.31954
SOLUTION: Since = 4 , ℎ = 3, then the values of () = at =0,3 ,6,9, 12 are given
Y
Y_8 :
below
0 3 6 9 12
() 1 0.1 0.02703 0.01220 0.00690
=1.46174
SOLUTION: The given data is nine that is we have eight sub intervals which is a multiple of four.
Thus we can use the Boole's rule.
or I = U è¤
äJ
ä
a=
I = UT I è¤
VT
+14(37.75) +7(50.67) ]
= 30.87 m/sec
Here we assume that within any six consecutive subintervals of width ℎ each, the interpolating
polynomial approximating () will be of degree 6. Substituting = 6 in the General quadrature
formula given by eq. (7.3) and neglecting all differences above, we have
Eq. (7.18) is known as Weddle’s rule. Weddle’s rule was found to be more accurate than most
of the other rules.
REMARK
o Jª (Ẋ).
M{_
YAT
i. The error estimate is given by
SOLUTION: Dividing the interval [0, 1.5] into 6 equal parts of each of width ℎ =
[Link]
E
=0.25
>(T.;C)
YT
= [0+5(0.0549) +0.1927+6(0.3777)+0.5820 + 5(0.7843)+ 0.9694]
given below
Here the number of sub interval with h=1/6 is six, which is a multiple of six, so that we can use
Widdle's rule.
EXERCISE –7.1
We modify the Trapezoidal rule to find a better approximation to the value of an integral. We
know that the truncation error in the trapezoidal rule is nearly proportional to ℎ2 an interval of
size ℎ. The error in the Trapezoidal rule
ℎ; o ÇÇ ()
(M?)
Y;
E=
o ÇÇ ()
(M?)
Y;
If we put c=
By the Trapezoidal rule with two different intervals say ℎY ,ℎ; . Let Y and Y denote the
approximate values with the corresponding errors E1and E2 respectively.
Then
I= Y + ¦ℎY ; and I = ; + ¦ℎ; ;
Therefore, we get
Y + ¦ℎY ; = ; + ¦ℎ; ;
ae Ma:
)ℎY ;
{: : M{e :
Or c=(
Hence
ae Ma:
I= Y+({ : : )ℎY ;
: M{e
ae {: : Ma: {e :
{: : M{e :
Or I= …………….(7.20)
This will be a better approximation to than Y or Y . The above method is called Richardson's
method.
We get
: e
ae Ma:{: a Ma Aa: Mae
@ @ e :
: =
>
M{:
I= = =
@ @
reduced by a factor . Let åY , å; , å> … denote the results. Let the formula eq. (7.21) be applied
Y
A
åY å; å> åA ……
øY ø; ø> ……
ùY ù; ……
The above computation is continued with two successive values which are very close to each
other. This refinement of Richardson’s method is known as Romberg integration. The values of
the integral, in Romberg integration can be tabulated as follows.
(ℎ)
(ℎ, ℎ/2).
(ℎ/2) (ℎ, ℎ/2, ℎ/4)
(ℎ/2, ℎ/4) (ℎ, ℎ/2, ℎ/4, ℎ/8)
(ℎ/4) (ℎ/2, ℎ/4, ℎ/8)
(ℎ/4, ℎ/8)
(ℎ/8)
Where
(ℎ/4, ℎ/8) = 4 −
Y { {
> V A
…………………………………………….
Y.; ä8
EXAMPLE 7.14: Using Romberg’s method compute I =UT Y_8
correct to 4 decimal places
() = Y_8
Y
SOLUTION: here
(ℎ) = (0.6) =
T.E
;
(1+2(0.6256) +0.4545)
=0.8113
With ℎ =
T.E
;
= 0.3, we get
I(ℎ/2) = (0.3) =
T.>
;
(1+2(0.7692+0.6250+0.5263) +0.4545)
= 0.7943
With ℎ=
T.E
A
=0.15, we get
(ℎ/4) = (0.15) =
T.>
A
(1+2(0.8695+0.7692+0.6896+0.6250+0.5714+0.5263+0.48780) +0.4545)
=0.7899
Now
Similarly
Therefore
I= UT
Y.; ä8
Y_8
=0.7883
SOLUTION:
Here () = §¬
We can take ℎ = //2, //4, //8,
i.e ℎ = //2, ℎ/2 = //4, ℎ/4 = //8
Using Trapezoidal rule with h= //2 , we get
/b /
0 2
() 0 1 0
With ℎ =
A
, we get
/b /b 3/b /
0 4 2 4
() 0 0.707 1 0.707 0
With ℎ =
V
, we get
/ / 3/ / 5/ 3/ 7/
(ℎ/4) = (/ /8 = (0) + 2{ + c d + + c d + c d + c d} + (/)
16 8 8 2 8 4 8
=1.9747
Now
4(1.896) − 1.571
Y
>
=
= 2.004
Similarly
I(; , A ) = I( A , V ) => 4 V − A
{ { Y
= 4(1.9747) − 1.896
Y
>
=2.0009
And
= 4(2.0009) − 2.004
Y
>
= 1.999
Therefore, = UT §¬è= [Link] exact value is 2.
Y ä8
EXAMPLE 7.16: Apply Romberg’s integration method to findUT
Y_8 :
, correct to four decimal
SOLUTION: Let () = . Applying the trapezoidal rule, for ℎ = 0.5, we obtain
Y
Y_8 :
0 0.5 1
() 1 0.8 0.5
(ℎ/2) = UT
Y ä8 T.;C
Y_8 : ;
Hence, = [1+2(0.9412+0.8+0.64) +0.5] =0.7848
Hence,
= UT
Y ä8
Y_8 :
= 0.7855 correct to four decimal places.
ℎ= ∶ = () + () = [ 1+
Y { Y Y/;
e
; ; A [ß^ ( )
]= 0.51072
:
)+(
Y
ℎ=
Y @ :
A V [ß^@
e
[ß^:
e
: I = [ 1+2( ]= 0.507988
e : = e
+ + )+(
Y
ℎ=
Y ç ç ç :
V YE [ß^ç
e
[ß^ç
:
[ß^ç
=
[ß^:
e
:I= [1+2( ] = 0.507298
+ + + +
Y
ℎ=
Y eD eD eD eD eD
YE >; e
[ß^eD
:
[ß^eD
=
[ß^eD
@
[ß^eD
B
[ß^eD
:I= [1+2(
D _ e
+ eD
D + eD
_ )+( :
e
[ß^eD [ß^eD [ß^:
] = 0.507126
ℎ/8
0.507068
0.507298 0.507069
ℎ/16
0.507069
0.507126
Hence,
EXERCISE 7.2
As mentioned earlier, the limits of integration for Gauss-Legendre integration rules are [– 1, 1].
eq. (7.23) can have (2 + 2) unknowns. Hence the Gaussian method can be made exact for
polynomials of degree less than or equal to 2 + 1 and its maximum order is 2 + 1.
Therefore, we transform the limits , to [– 1, 1], using a linear transformation.
Let the transformation be = ͤ + g.
When = , we have ¤ = – 1: = – Í + g.
When = , we have ¤ = 1: = Í + g.
Solving, we get
Í= ,g =
M? _?
; ;
Where c is called the error constant. Then, the error term is given by
h
= (P_Y)! (P_Y) (), a<<b ………….(7.28)
If b ( P_Y )also becomes zero, then the error term is obtained for () = P_; .
Gauss one point rule (Gauss-Legendre one point rule)
The one point rule is given by
() = 1: UMY è = 2 = eT
Y
Error of approximation
The error term is obtained when () = ; . We obtain
c= UMY x ; dx − 0 = 2/3
Y
Remark :
Since the error term contains ÇÇ (), Gauss one point rule integrates exactly polynomials of
degree less than or equal to 1. Therefore, the results obtained from this rule are comparable
with the results obtained from the trapezium rule. However, we require two function
evaluations in the trapezium rule where as we need only one function evaluation in the Gauss
one point rule. If better accuracy is required, then the original interval [a, b] can be subdivided
T ; = , 0r T = ± = −Y
Y Y
> >
Error of approximation
The error term is obtained when () = 4. We obtain
c= UMY x A dx − ¼U + U½ = C − U =AC
Y Y Y ; ; V
where eT ≠0, eY ≠0, e; ≠ 0 and T ≠ Y ≠ ; . The method has six unknowns eT , eY , e; ,T Y , ; .
Making the formula exact for () = 1, , ; , > , A , C , we get
Error of approximation
The error term is obtained when () = 6. We obtain
Remark:
Since the error term contains (6) () , Gauss three point rule integrates exactly polynomials of
degree less than or equal to 5. Further, the error coefficient is very small (1/15750≈0.00006349).
Therefore, the results obtained from this rule are very accurate. We have not derived any
Newton-Cotes rule, which can be compared with the Gauss three point rule. If better accuracy is
required, then the original interval [a, b] can be subdivided and the limits of each subinterval can
be transformed to [– 1, 1]. Gauss three point rule can then be applied to each of the integrals.
Table 7.1: Nodes and Weights for the Gauss-Legendre Integration Methods
n nodes Weightse
1 ± 0.5773502692 1.0000000000
2 0.0000000000 0.8888888889
± 0.7745966692 0.5555555556
3 ± 0.3399810436 0.6521451549
± 0.8611363116 0.3478548451
4 0.0000000000 0.5688888889
± 0.5384693101 0.4786286705
± 0.9061798459 0.2369268851
5 ± 0.2386191861 0.4679139346
± 0.6612093865 0.3607615730
± 0.9324695142 0.1713244924
where (¤) =
V(_ >)
YE_(_>)@
= 2 (0) = 2[YE_VY ]
;A
= UX= 0.494845
AV
Y
= U [5(0.439299) + 8(0.247423) + 5(0.137889)
= 0.540592
The exact values is = 0.540420.
The magnitudes of the errors
In the one point, |0.540420 − 0.494845|=0.045575
= UY = UMY (¤)è¤
; ä Y
_>
where (¤) =
Y
_>
Y
= U [[5(0.449357) + 8(0.333333) + 5(0.264929)]
= 0.693122
The exact solution is = (2) = 0.693147.
The absolute error in the three point Gauss rule is 0.000025.
Y
=; [0.56558356 – 0.15856672]
= 0.20350842
Using the three point Gauss rule, we obtain
Y
= YV [– 1.26018516 + 1.41966658 + 3.48936887]
= 0.20271391
EXERCISE 7.3
Y ä8
1. Using three point Gauss formula, evaluateUT
Y_8 :
.
Y ä8
2. Apply Gauss two point formula to evaluateUMY
Y_8 :
.
Can be reduced to
UMY Ê = ∑b)T e Fl
Y ¸(8)ä8
−±]
……………(7.47)
The abscissa (nodes) ±m and weights nm are unknowns. The nodes ±m ′s are the zeros of the
Chebyshev polynomial
b_Y () = ¦R§ (( + 1) ¦R§ MY ) ……………(7.48)
The first few Chebyshev polynomials are given by
T () = 1, Y ( ) = , ; () = 2 ; − 1
> () = 4 > − 3, A ( ) = 8 A − 8 ; − 1
The Chebyshev polynomials are orthogonal on [– 1, 1] with respect to the weight function
S() =
ÊM±]
. The methods of the form eq. (7.47) are called Gauss-Chebyshev integration
The weights nm′s in eq. (7.47) are equal and are given by
e = ^_Y
o
, where = 0,1,2, … , …………..(7.50)
I= UMY = eT yT
Y ¸(8)ä8
Ê M±]
…….……..(7.51)
UMY + UT
MY T äq Y äq
; √E √E
=
T Y
= -À√1 − ; Á)]MY − À√1 − ; Á]T
= 0 = /T ⇒ T = 0
Therefore
= 2 UT
Y 8 : ä8
ÊM±] ;
=
= UMY = eT yT + eY yY
Y ¸(8)ä8
Ê M±]
..….……..(7.52)
I = UMY = eT yT + eY yY + e; y;
Y ¸(8)ä8
Ê M±]
..………..(7.53)
I= UMY = (− ( )]
Y ¸(8)ä8 √> √>
Ê M±] > ; ;
)+(0)+
I= UT è
Y p^ (8_Y
Ê8(YM8)
e =b_Y , = 0,1, … . , ,
= 1: = ; (− ) + ( ) ] = 1.184022,
Y Y
√; √;
EXERCISE 7.4
Using the trapezoidal rule again to evaluate the integrals in eq. (7.55), we obtain
If discrete value are given we can use the Trapezoidal rule by dividing the interval [a, b] into
4
equal subintervals each of length ℎ and the interval [c, d] into ¯ equal subintervals each of
length . We have
¬ = 0 + ¬ℎ, 0 = and =
os = o0 + s, 0 =¦ and o¯ =è
SOLUTION: Using Trapezoidal rule, with ℎ = = 1/2, the nodal point are (1, 1), (3/2, 1), (2, 1),
(1, 3/2), (3/2, 3/2), (3/2, 2), (1, 2), (3/2, 2), (2, 2).
The values of the integrand at the nodal points are obtained as the following.
y\ 1 1.5 2
1 0.5 0.4 0.33333
1.5 0.4 0.33333 0.285714
2 0.33333 0.285714 0.25
We have
= 0.343303
With h = 0.25, we have the nodal points as shown below
following data.
SOLUTION: We have the step lengths along -axis and o −axis as ℎ = 0.5 and = 1.0
respectively. Also, the number of intervals along x-axis and y-axis are ¯ = 4 and = 2. We
have the following grid.
a. UY UY
Y.C ; ä8är
8_r
.
b. UT UT 7 8_r èèo.
Y Y
c. UY UY
C C ä8är
Ê8 : _r :
,with two and four subintervals
Divide the interval , into 2 equal parts each of length ℎ = (– )/(2). Divide the
interval [c, d] into 2¯ equal parts each of length = (è– ¦)/(2¯). We have odd number of
points on each mesh line and the total number of points, (2 + 1) (2¯ + 1) is also odd.
For ¯ = = 4, the weights in the Simpson’s 1/3 rule are given in Figure.
õ)Y{ fÀT , o;õMY Á + 4 ∑ª)Y (;ªMY , ;õMY )+2 ∑ª)Y (;ª , õMY )+ (;b , ;õMY )}
+ 4∑° b bMY
õ)Y{ fÀT , o;õ Á + 4 ∑ª)Y (;ªMY , ;õ )+2 ∑ª)Y (;ª , ;õ ) + (;b , ;õ )}
+2∑° b bMY
EXAMPLE 7.24: Evaluate the integral I=UY UY , using Simpson's 1/3 rule with ℎ = = 0.5.
; ; ä8är
8_r
SOLUTION: Using Trapezoidal rule, with ℎ = = 1/2, the nodal point are(1, 1), (3/2, 1), (2, 1),
(1, 3/2), (3/2, 3/2), (3/2, 2), (1, 2), (3/2, 2), (2, 2).
The values of the integrand at the nodal points are obtained as the following.
y\ 1 1.5 2
1 0.5 0.4 0.33333
1.5 0.4 0.33333 0.285714
2 0.33333 0.285714 0.25
We have
= 0.33986
y\ 1 1.5 2
1 ½=0.5 2/5=0.4 1/3=0.333
1.25 4/9=0.444 4/11=0.3637 4/13=0.3077
1.5 2/5=0.4 1/3=0.3333 2/7=0.2857
õ)Y{ fÀT , o;õMY Á + 4 ∑ª)Y (;ªMY , ;õMY )+2 ∑ª)Y (;ª , õMY ) + f(;b , ;õMY )}
+ 4∑° b bMY
õ)Y{ fÀT , o;õ Á + 4 ∑ª)Y (;ªMY , ;õ )+2 ∑ª)Y (;ª , ;õ ) + f(;b , ;õ )}
+2∑° b bMY
Hence
0.125
I= [ f (1,1) + f (2,1) + f (1,1.5) + f (2,1.5) + 4{ f (1.5,1) + f (1,1.25) + f (1.5,1.5) + f (2,1.25)}
9
+ 16 f (1.5,1.25)]
0.5 + 0.333 + 0.440.2587 + 4{0.4 + 0.444 + 0.3637 + 0.3077} + 160.3637
T.Y;C
U
=
=0.184432
The exact solution is = 3.5 (3.5) + 4 − 27 − 2.52.5 = 0.184401
EXERCISE 7.6
Evaluate the integral using Simpson's rule
a. UT UT
Y Y ä8är
Y_8_r
.
c. UT UY ,taking ℎ = = 0.5
Y ; ;8rärä8
(Y_8 : )(Y_r : )
In this section we obtain triple integration formula by shifting the single integration formula.
We consider the evaluation of the triple integral
Using the trapezoidal rule again to evaluate the integrals in eq. (7.59), we obtain
Using the trapezoidal rule again to evaluate the integrals in eq. (7.60), we obtain
EXAMPLE 7.27: Evaluate the integral I=UMY UMY UMY (, o, p)èèoèp, using Trapezoidal rule if
Y Y Y
(, o, p) = ; o ; p ;
SOLUTION: Here T = oT = pT = −1, Y = oY = pY = 1, with = ¯ = Í = 1, then
functional values are (−1, −1, −1) = 1, (−1, −1,1) = 1, (−1,1, −1) = 1, (−1,1,1) =
1, (1,1, −1) = 1 (−1, −1,1) = 1 (−1,1, −1) = 1.
Therefore, using Trapezoidal rule
= 0.296296
Where = ( − )/2
Using the Simpson's 1/3 rule again to evaluate the integrals in eq. (7.62), we obtain
EXAMPLE 7.28: Evaluate the integral I=UMY UMY UMY (, o, p)èèoèp, using Simpson's 1/3 rule
Y Y Y
if (, o, p) = ; o ; p ;
SOLUTION: Here T = −1, Y = 0, ; = 1 with = ¯ = Í = 2, then functional values are
(−1, −1, −1) = 1, (−1, −1,1) = 1, (−1,1, −1) = 1, (−1,1,1) = 1, (1,1, −1) = 1,
(1, 1,0) = 0, (1, −1,0) = 0, (0, −1, −1) = 0, (0, −1, 1) = 0 (0, −1,1) =
0, (0,1, −1) = 0, (0,1,1) = 0, (0,1, −1) = 0, (1,0,1) = 0, (−1, 0, −1) = 0, (−1, 0,1) =
0, (0, −1,0) = 0, (0,1,0) = 0, (1,0,0) = 0, (−1,0,0) = 0, (0,0, −1) = 0, (0,0,1) =
0, (0,0,0) = 0
Therefore, using Simpson's 1/3 rule
EXAMPLE 7.29: Evaluate I=UMY UMY UMY (, o, p)èèoèp, using Simpson's rule if (, o, p) =
Y Y Y
=0.1852
+ o + p.
SOLUTION: Here T = −1, Y = 0, ; = 1 with = ¯ = Í = 2, then functional values are
(−1, −1, −1) = −3, (−1, −1,1) = −1, (−1,1, −1) = −1, (−1,1,1) = 1, (1,1, −1
= 1, (1, 1,0) = 2, (1, −1,0) = 0, (0, −1, −1) = −2, (0, −1, 1)
= 0, (0, −1,1) = 0, (0,1, −1) = 0, (0,1,1) = 2, (0,1, −1) = 0, (1,0,1)
= 2, (−1, 0, −1) = −2, (−1, 0,1) = 0, (0, −1,0) = −1, (0,1,0)
= 1, (1,0,0) = 1, (−1,0 = −1, (0,0, −1) = −1, (0,0,1) = 1, (0,0,0) = 0
Therefore, using Simpson's 1/3 rule
= 0.111
b) I = UT §¬dx
c) I = UMY è
Y Y
Y_8 :
2. Determine the area bounded by the curve () = 7 ;8 and the -axis between = 0 and
= 1 using the trapezoidal rule with an interval size of (a) ℎ= 0.5, (b) ℎ= 0.1. Determine the
relative error in each case given that the exact value of the integral I= 2.09726.
3. Evaluate the definite integral using Simpson's 1/3 rule and compare the two results with
analytical solution.
b
a. I= UT ;
√¦R§ è, taking = 6
c. I = UM> A è
>
c. I = UT è
C Y
A8_C
a. I = UT 7 M8 è
Y
b. I = UT 2 ; ¦R§ ; è
√
c. I = UT ( C + 3 > − 2) è
;
9. Use two point and three point Gauss formula to evaluate I= UT è. Compare with the
; Y
>_A8
exact solution.
10. Use two point and three point Gauss formula to evaluate I = UT
; ä8
8 : _;8_YT
11. Evaluate
a) UY UY èèo d) UT UT ; o ; èèo
C C Y ; >
8 : _r :
8_r
f) UT UT èèo
Y ; Y
UT UT( Y_8 : )(Y_r :)doè,
Y ; ;8r
c) with (Y_8r)
ℎ = = 0.25
12. Evaluate UT UT o(2 − ; )(2 − o) èèo, over the region shown.
; >
except possibly at w itself. Then a number x is the limit of v(±) as ± approaches w if for every
y > 0 there is a number . > 0 such that
0 < | I | < c ⟹ | I | <
If such a number can be found, we say that the limit of at exists.
In other words lim8⟶? () = lim8→?¥ () = lim8→? () =
(i.e the left and the right hand limits exist and are each equal to L). Otherwise, the limit does not
exist.
Continuity: A function is continuous at a point in its domain if
() = lim¥ () = lim ().
8→? 8→?
Intermediate value Theorem: Suppose is continuous on a closed interval, . Let be any
number between () and (), so that () ≤ ≤ () or () ≤ ≤ (). Then there exists
a number ¦ in , such that (¦) = .
Rolle's Theorem: If () is continues in the closed interval , , if ′() exists at each value of in
the open interval (, ) and () = (), then there is at least one value ¦ of in (a, b) such that
′(¦) = 0.
Mean value theorem: If f(x) is continuous in the closed interval [a, b] and f'(x) exists in the open
interval (a, b), then there is at least one value c of x in (a, b), such that
() − ()
= Ç (¦)
−
Taylor's series for a function of one variable: If f(x) is continuous and possesses continuous
derivatives of n order in an interval that includes x=a, then in that interval
( − ); ÇÇ ( − )b (b)
() = () + ( − ) Ç () + () + ⋯ + () + b ()
2! ()!
Where is the reminder term, can be expressed in the form
78= 78 U 7 8 è = 7 8 + ¦
ä
ä8
8 = 8 U 8 è = b? + ¦
ä ?
ä8
4. square matrix : If a matrix has equal number of rows and columns we call the matrix as
square matrix.(i.e. An ¯ × matrix å = Àªõ Á is called a square matrix if and only if
¯ = . We represent it as Main diagonal
6. A diagonal matrix whose diagonal elements are equal is called scalar matrix.
8. Let A be × matrix.
If all elements of A below the main diagonal are zero. Then A is called upper triangular
matrix. i.e if ªõ = 0 for ¬ > s.
If all elements of A above the main diagonal are zero. Then A is called lower triangular
matrix if and only if ªõ = 0 for ¬ < s.
Matrix Operations
Addition and scalar multiplication of matrices
DEFINITION: Let å = Àªõ Á and ø = Àªõ Á be two ¯ × matrices . Then
å + ø = (ªõ + ªõ ) for all ¬ and s.
å = ( ªõ ) for all ¬ and s and ∈ ;.
YY Y;
Letå = ;; . Thenè7¤(å) =
;Y
The minor of an element is the determinant obtained from given determinant by deleting the
row and column that contains the element and the minor of an element ªõ is denoted by Wªõ
YY Y; Y>
;> ;Y ;>
;; ;> ÷.Then WYY = L ;;
Let å = ö ;Y >; >> L, WY; = L>Y >> L, …………
>Y >; >>
The cofactor of an element ªõ which is the product of the minor of ªõ and (−1)ª_õ and denoted
by ùªõ
i.e ùªõ = (−1)ª_õ Wªõ
YY Y; Y>
;>
Let å = ö ;Y ;; ;> ÷. Then ùYY = (−1)Y_Y WYY = L ;;
>; >> L,
>Y >; >>
Prepared by Haileyesus T., UOG Page 270
;Y ;>
ùY; = (−1)Y_; WY; = − L >> L, ……….
>Y
REMARK: The cofactor and minor of an element ªõ differ only in sign i.e ùªõ = ±Wªõ .Hence foe
each element ªõ the sign of the cofactor is given as
+ − +
− ....
...
+ − + ....
+ −
.
= YY (;; >> − ;> >; ) − Y; (;Y >> − ;> >Y ) + Y> (;Y >; − ;; >Y )
=YY ;; >> +Y; ;> >Y +Y> ;Y >; − (YY ;> >; +Y; ;Y >> +Y> ;; >Y )
DEFINITION: (Adjoin of a Matrix) Let A be an matrix and B be a cofactor matrix of A. Then the
adjoint matrix of A written as Adj (A) is the transpose of B. i.e. åès (å) = B
Properties of Adjoint of a matrix
1. If A is a square matrix of order n, then
å (åès (å) = è7¤ (A) b = (adjA)A
2. If A and B are two squre matrices of the same order then
åès (åø) = ( åès ø)(åès å)
Solutions of chapter 2
Exercise 2.1
1. 2.6875 2) 0.370 3) 2.88
Exercise 2.2
1. -1.92 2) 2.71 3) 1.571
Exercise 2.3
1. 0.5177573637 2) 4.5616 3) 1.114
Exercise 2.4
1) −0.682327803 2) 5.477225575 3) –0.420365
Exercise 2.5
1) 2.279 2) 0.5885 3) 1.442 4)1.973
Exercise 2.6
1) 2.15466 2) 1.4879
a) = XY ,o = XY ,p = XY b) = 2, o = −1, p = ; c) = 2, o = 3, p = 6
YYX MVY YAV Y
d) = 1, o = 2, p = −1, S = −2
Exercise 3.2
Exercise 3.4
) = 1, o = – 5, p = 1 b) = 1, o = – 1, p = – 1, S = 1
Exercise 3.5
a) = 1, o = 1, p = 1, b) = 1, o = 2, p = 3, c) = 1, o = 1/2, p = – 1/2
Exercise 3.6
a) = 1, o = −1, p = 1 b) = 1.08, o = 1.95, p = 3.16
Exercise 3.7
a) = 1, o = 1, p = 1 b) = 3, o = 2, p = 1 c) = 1, o = −1 , p = 1
Review exercise
1. a) = −1 y = 2, z = 3 b) = 1.5205, o = – 1.0658, p = – 0.0356
= 1, o = 2, p = – 1 b) = 2, o = 2, p = 3
3. a)
= 2, o = 0, p = – 3 b) = 0, o = 1.5, p = 1
4. a)
= 1, o = 2, p = 3 b) = 1, o = – 1, p = 1
5. a)
= 4, o = – 1, p = 3 b) = 1, o = 2, p = 3
6. a)
7. a)
Solutions of chapter 4
Review exercise
2. a)7 ?8 (7 ? − 1) b) §¬(¦R§1 − 1) − §¬1¦R§ c) 8(8_Y)(8_;)
MA8_;
d)
o ⊖; o ⊖> o
3.
ϴo
1 6
9
2 15 12
21 6
3 36 18
39 6
24
63
4 75
5 138
4. ∆AoT = −49
5.
o δo δ; o δ> o
0 1
3
1 4 10
13 6
2 17 16
29
3 46
1 71
22
3 115 17
90 -14/5
5 295 3
117
6 466
7. a) 77 b) 35 c) 3 and 37
8. correct value 63, error 7 =2
Solutions of chapter 5
Exercise 5.1
1) 18.149188 2) 16.25≈ 16
Exercise 5.2
1) ( 3+13 2+ 56 + 28)/2; 100 2) 193.83336
Exercise 5.3
1) 0.15907168 2) 394.6875 3) 32.9453125 4) (31) = 4913, (32) = 5052
Exercise 5.4
Exercise 5.5
Review exercise
1) 23– 7 2+ 6 + 1 2) (i) 28 (ii) 36 3) 1.25 4) 31
5) 103– 272+ 3 + 35; 71 6) 10 7) 44/3 8) 36.05469
9) 36.625 10) 478.375 11) 394.6875 12) 4913, 5052
Solutions of chapter 7
Exercise 7.1
1) 0.2600 2) 0.9871 by trapezoidal , and 1.001 by simpson's rule
3) 53.6540 4) 0.8814 5) 0.1577
Exercise 7.2
3) 0.1577 4) 6.3658
Exercise 7.3
1) with h= π/6, I=1.954097, Exact : 2.0 2) with h= 2, I=360 3) with h= 1, I=35.65
Exercise 7.4
1) for n=1, 0.244956, n=2, 0.355464, n=3, 0.351617
Exercise 7.5
d. 0.178432 b) 3.07627 c) With h= k= 2, I= 4.1345With h= k= 1, I= 3.9975
Exercise 7.6
a) 0.524074 b) 0.232316 c)0.312330
b)
11) (a) when n= 2, I= 4.134 and n= 4, I= 3.997 (b) 0.31233 (c) 0.48997