18.

98 Engineering Mathematics

derivatives of f (x) at x = x/^ are compared. This gives the required number of equations for the
parameters to be determined. Solution of this system of equations gives the required method. The
leading non-vanishing term gives the error term. We illustrate this method through the following
example. Consider the method

= 7 (Xk + h} + bf {x,,} + cf{Xk- /z)].

h
Writing the Taylor series expansion about x^, we get

1 h2 h,3 /zz,

A 2 6

h,2 h3
+ bf{x,} + c {f{x,)-hf'M + — f"(x,)-—f"\x,) + ..}
2 6

2 {a+b + c') f{xi,} + h{a- c} /'{x^) -h (a -t- c) /"(x^.)

h

h.3
+ _(a-c)r'(xj
6
+ ... .
We have three parameters a, b, c to be determined. Comparing the coefficients of/(x^.),
f'\Xk\ we get
a + b + c = 0, a - c = 1, a + c = 0,

whose solution a = \I2, b = Q, c = --\I2.

The method is

In
The error term is given by
3
- -i h-
^*+1-
h 6 6

The order of the method is 2.

18.6.2 Numerical Integration Methods

The problem of numerical integration is to evaluate the integral

r* f{x}dx
Ja
'a
(18.135)

where/(x) is either given explicitly or defined by a data of n + 1 values (x„ /{Xj)}, z = 0, 1/2, ... . n.
We assume that the tabular points are equispaced with spacing h = (b - a^/n. The points are given

1
 Numerical Methods 18.99

^0 + Ih, i \, 2, ..., n and x,^ - b. It is possible to construct numerical intergration

methods when the limits are finite or one of the limits is infinite or both the limits are infinite. The
integral I is approximated by a linear combination of the values of f(x} at the tabular points as

1= JJa f{x)dx = 4/(xo) + A, /(%,) + ... + A„/(x„).

^<7
(18.136)

The tabular points x^.’s are called nodes or abscissas and coefficients A^’s are called weights of the
integration method. The method (18.136) is also called a quadrature rule or an integration rule.

We define the error of approximation as

ZZ

R,n’ (18.137)
'a
k=^
Order of a method

An integration rule is said to be of order p, if it produces exact results for all polynomials of degree
<p. That is, it produces exact results for/(x) = 1, x, ..., x^. Therefore,
n
,(/x")=_(’x"*-£A.xr = o,
R.n for m = 0, 1, 2, ..., p.
'a
A=0

The error term is obtained for /(x) = x^*'. We define

n
b
c = £ f{x)dx - (18.138)
'a
k=0
where c is called the error constant. Then, the error term is given by

c (p+i)
R„{f,x) = f (^)
a e<6. (18.139)
(p + 1)!

18.6.2.1 Newton-Cotes quadrature rules

Consider the case when the abscissas are prescribed and are equi-spaced. That is, the interval (<3, Z?)
is divided into n equal parts. The mesh spacing is given by h = (b - a)ln, and Xq = a,
Xj = Xq + ih, i = 1, 2, 3, ..., n and x„ b. Then the integration rule (18.136) has n + 1 unknowns
A), A,, . A„. Replace the integrand of the rule (18.136) by the Lagrange interpolating polynomial
interpolating at the n + 1 point (x^,/(), A: — 0, 1, 2, 3, n. We obtain
n (»h
>h n

1 = f f{x} dx= \ Y h^^yfkd^ i fk f h = 'L aa- fk^ (18.140)

J'a
a a
J'a h-n
k=Q k=0
t=0 k=0

where
(x-Xo)(x-Xi)...(x-Xjt_|XJC-Xjn.i)...(x-x„)
(Xa- — Xq) {Xk “ -^l)" •(■^A ~ ■’'•A-l) ^^k ~ + • \^k ~ )

are the Lagrange fundamental polynomials, and

2
18.100 Engineering Mathematics

pb
J
Ja'a
dx. (18.141;

Let X = Xq + sh. Then,

(5/;)[(5-1)A]...[{5-(A-1)M][{5-(A-H)}H--[{^-"}^]
W =
(AA) [(A: - 1)A]... (A) (-A) (-2/z)... [(« - A:)/?]

5 (5 - 1) (s - 2) ... (5 - k + 1) (5 - A: — 1) ... (s — n). (18.142)

A!(«-A)!

The error term can be obtained by using (18.138) and (18.139)

Trapezoidal rule

Consider the case n = 1, that xs h = b - a. The data values are (a, f {a}}, {b, f {b}}, and

[ f{x) dx = (a) + Ixfib}.

Ja
We obtain from (18.141) and (18.142)
4 , h 1
J b
= (5-1) ds = ky■1 = h s ds = —
Jo 2 2
The rule is given by

dx=^ =
[/(«) +/(6)]. (18.143)
Ja'a 2
This rule is called the trapezoidal rule.
<•/>
Geometrical interpretation We note that, f (x) dx defines the area under the curve y =f{x}, above
Ja 'a
the x-axis, between the lines X - a and x = b. The end points on the curve are {a, f (a)).
{b,f{by), [SQQ Fig. 18.14). The right hand side of (18.143) is the area of the trapezoid ABQP.
Therefore, trapezoidal rule approximates the area under the curve by the area of the trapezoid.

4
X
O a b

Fig. 18.14. Trapezoidal rule.

3
 Numerical Methods 18.101

Remark 12

The rule can also be obtained by replacing the integrand / (x) by the Newton’s forward difference
interpolating polynomial up to the A/q term.

Error term We first show that the trapezoidal rule integrates exactly polynomials of degree < 1,
that is, R^(f, x) = 0 for/(x) = 1 and x. Now, substituting/(x) - 1, x in (18.143), we get

/(x)=l: b- (2), which is true.

0

/(x) = x; {b^ - = (b-a) (a + Z)) = -^ {b^ — which is true.

1
Let /(x) = x^. Then, from (18.138), we get

•6
{b-a} + ^2) = 1 (/,3 - a^) _ 1 - ab^ - o^)
c = dx -
2

- ?>a^b + 3ab^ ~ b^) = - ^ {b -

The expression for the error is given by

(b-a)^
«,x) = ^r(e) = (18.144)
i2 12

(b-a?
and l^i(/x)|< ^2, where M2= max |/''(x)|. (18.145)
12 o<.v<Z)
Hence, the method is of order 1.

Composite formula If the length of the interval [a, Z?] is large, then b <3 is large and the error
estimate (18.144) becomes meaningless. Therefore, to have meaningful results, we subdivide [a. Z?]
into a number of sub-intervals of equal length and apply the trapezoidal rule to evaluate each integral.
Let [a, b} be sub-divided into N parts of equal length h.

Then, h = {b - a^/N. Let the points be

^0’ X] = Xq b. ^2 = ^0 .. • > ■’‘^0 + = b.

We write

f f{x}dx=
Ja
'fl Jxn
'JCo
f{x}dx=
Jxa
-<0
f
■^2

x.
•’'^1
/{x^dx^ .. . + r
J'x,
x.,,
W-l
dx.

Now, using the trapezoidal rule to evaluate each integral, we obtain

f{x-}dx = [{/(xo) +/(x,)} + W.) +/(X2)} + ... + {/(^v-i) +/Uv)}]

Ja 2

4
18.102 Engineering Mathematics

= h [/(Xo) + 2 {/(X,) +/(X2) + ... +/(X;v-l)} +/(Xaz)]- (18.146)

This formula is called the composite trapezoidal rule and is of order 1.

The error expression (18.144) becomes

h.3
Rt(f, ^) = - ^
12
[/"<41) +/"(& + - in ^N-

Let = max |/"(x)|. Then

[ro,-v„l

\Rx{f, x)i [i/"(e,)i + i/"(e2)i +... + i/''(e7v)i]

A7i,3 _ (b - a)h^ (b-a)^
< —M2 = A/2 M2 (18.147)
12 12 12N^

since Nh = b - a. This expression is a realistic representation of the error in the trapezoidal rule. If
an error tolerance e is prescribed, it is possible to find N, the number of sub-intervals required to
achieve that accuracy. We get from (18.147)

(b-a)^ (Z> - fl)^ M2

M2 £, or #2 > (18.148)
nN^ ne

Examplel8.40 Evaluate the following integrals using trapezoidal rule with N= 1, 4. Compare with
the exact solution.

4 dx (ii) f dx
(i)
'0 3 + 2x ’0 x^ -l-2x -I-10
JO

Find the bound on the error. Also, find the number of sub-intervals required if the error is to be less
than 5 x lO"*.

Solution

(i) We have a — 0, h — 1 and /(x) — 1/(3 -t- 2y). For N ~ 2, we get h = 1/2. The abscissas are
Xq = 0, x,= 1/2 and X2 = 1.0. Hence,

/z 1 _1_ 1 nfl 1 1
/=- /(0) + 2/ - +/(1) — -l- 2 — -(- = 0.25833.
2 2 4 3 4 5

For N = 4, we get h 1/4. The abscissas are = 0, x, = 1/4, X2 = 1/2, X3 = 3/4 and X4 = 1.0.
Hence,

+ /(1)

5
 Numerical Methods 18.103

1 1 2 1 2 1
-- + 2^- + - + -^ + - = 0.25615.
8 3 7 4 8 5

The exact solution is 0.5 In (5/3) 0.25541. The actual errors in magnitude are respectively
0.00292 and 0.00074.

1 2 8
We have /(x) = /'(X) = - /"(x) = and
3 + 2x ’ (3 + 2x)2 ’ (3 + 2x) 3

8 8
M2 = max
[0,1] (3 + 2x) 3 21

{b-a}h^ 2h^
Hence, | Error | < M2 = -----
\2 81
For h = 1/2, I Error] < 0.00617, and for h = 1/4, | Error] < 0.00154.

Note that the actual errors are smaller than the maximum magnitude errors.
If e = 5 X io~*, then the number of sub-intervals required is given by

> {b - M^ jO^ 8
= 49.38 or TV > 7.03.
12£ 60 127

Since TV is an integer, we require TV = 8.

(ii) We have a - 0, b = 2 and /(x) = \l{x^ + 2x + 10). For TV = 2, we get A = 1. The abscissas
are Xq = 0, x, = 1 and X2 = 2. Hence,

h 1 1 2 1
1= - 1/(0) + +/(2)] = - + + = ^.15470.
2 2 10 13 18

For TV = 4, we get h = 1/2. The abscissas are Xq = 0, X] = 1/2, = 1, X3 = 3/2 and X4 =2.
Hence,
Ar f /' i ^3^1
/(0) + 2 / 1 +/(1) + /^ +/(2)

= - [0.1 + 2(0.08889 + 0.07692 + 0.06557) + 0.05556] = 0.15458.

4

, . . . \ n -1
The exact solution is —----- tan = 0.15455.
3 4 3

The actual errors in magnitude are respectively 0.00015 and 0.00003. We have

1 - 2(x + l) 6(x^ +2x-2)

, /"(X) =
x^ +2x +10 (x2+2x + 10) 2 (x2+2x + 10)’

6
18.104 Engineering Mathematics

max Ix^ -I- 2x - 2| = 6 and min (x’ + 2x + 10)^ - 1000.

[0,2] [0,2]

36
Hence, M^ < TTrrr-
1000

{b — a'jh^ 2 f 36 >,2 3
Therefore, | Error | < = =------ h ■
12 12 <1000 J 500

For h= lErrorl < 0.006 and for h = 1/2, |Error| < 0.0015.

If £ = 5 X 10'*, then the number of sub-intervals required is given by

{b-a}^ M.^ _ 8 36
n'^ > 10* = 48, or N-‘I.
12e 60 <1000

Simpson’s rule (Simpson’s 1/3'^'* rule)

We have shown that the trapezoidal rule of integration is of order 1, that is it integrates exactly
polynomials of degree < 1. However, in many applications, we require methods that give higher
accuracy. One such method is the Simpson’s rule of integration.

Consider the case n = 2, that is Zi = (Z? - a)/2. The abscissas are Xq = a, = Xq + h = (a + 6)/2,
M = b. Integrand /(x) is approximated by a quadratic polynomial. From (18.140), the rule is
written as

ayb
f
Ja
f{x) dx= X^f^a} + kJ + A. /(6).

We obtain from (18.141) and (18.142)

■2
h £(^-l)(5-2)rf5=p
h
- A,■1 f 5 (5 - 2) ds = —
Jo ?>
,

h h

The rule is given by

Jffc /(x) h
= - [< (Xo) + (Xi) +/(X2)]

(b-a} a + b]
6 ( + /(6) . (18.149)

This rule is called the Simpson’s 1/3’^'^ rule or simply Simpson’s rule.
Geometrical interpretation The abscissas are Xg = a, x. {a + b}/!, Xt = b. The points on the curve
are P{a, j («)), QGa + b}l2,f {(fl + 6)/2’f), and R^b^f {b}}. The area under the curve y = f{x} , above

7
 Numerical Methods 18.105

the A'-axis, and between the lines .v a and .v - b is approximated by the area under the parabola
passing through the points P, Q, R.

Remark 13

The rule can also be obtained by replacing the integrand /(a) by the Newton’s forward difference
interpolating polynomial up to the term.

Error term We first show that the Simpson's l/3rd rule integrates exactly polynomials of degree
< 3, that is, /?2 (./' 0 for /’(a) 1. .V, .V and A-'. Substituting /'(a) = 1, a, x~ and a'^ in (18.149),
we get

(b-a}
/(X) = 1: Z? - a = (6), which is true.
6

I = (b-a) a+b 1
/(x) = x: a+ 4 +b = — which is true.
6 2 2

/w =
(b-a}
b
a^ +4
a+b
2
I +6' 2

= (^22^ [a^ + fl6 + 6^] = -J (b^ - a^}, which is true.

3

3
.3. /j,4 4x _ (^-q) a+b I 3
./’(■V) = X : - (b - a }-^ +4
4 6 2

= + a"b + ab^ + 6^] = (b^ - a^}, which is true.

Let f(x} = xl Then, from (18.138), we get

4
(b-a} a+b I +Z1'.4
C x^dx- o'*+4
’a 6 2

(b-a} + 40^6 + ba^b^ + ^ab^ +b'^} + b^

= I (*’ -«’) - 6 4

= _L r24(Z?^ - a^} - 5(b - + 4a^Z? + 6a^b^ + Aab^ +56'*}]

120

(b-a} \b!^ - ^ab^ + 6a^h^ - 4a^b + a"^] - -

120 120

The expression for the error is given by

8
18.106 Engineering Mathematics

c (<^)
(b-af
(18.150)
7?2C/; %)
4! ■' 2880 90

since h = {b - a)/2 and Xq < X2-

Therefore, Simpson’s rule is of order 3, which is two orders higher than that of the trapezoidal rule.

Composite rule As in trapezoidal rule, if the length of the interval [fl, 6] is large, then (b - a} is
large and the error estimate (18.150) becomes meaningless. We note that the application of Simpson’s
rule requires three points. Therefore, we subdivide [fl, A] into an even number of sub-intervals. Let
the number of sub-intervals be 2N so that h = {b - fl)/(2A), and we have odd number of points. Let
the abscissas be

fl = Xq, ^1 ~ ^0 ■’^2 ^0 2/z,..., X2N = Xq + TNh = b.

We write

'6
p"
ex2 fx^

J fl
f{x)dx =
Jx„
f(x)dx =
Jao
f{x)dx +
4*2
f(x)dx + ... +
'^^2N-2
f(x)dx.

We note that there are N integrals. Now, using the Simpson’s rule to evaluate each integral, we obtain

h
I f(x)dx = - + 4/(x,) +/(Xj)) + {f{x^ + 4f(x,) +f(x,)}
Ja j

■*■••• 4/(X27V-1) +Z(-1^2y)}]

= I [/(Xo) + 4{/(x,) +/(X3) + ... +/(X2„_,)|

2 tf(x2) + f{x^) + ... +Z(X2;v2)} .A(-^2A')]' (18.151)

This formula in called composite Simpsons rule.

The error expression (18.150) becomes

l5
«2(/. «l) + A’«2) + - + A’ (?,v)]

where Xg e. JC2 < ^2 ^4 etc. Let

A/4 = max I/*'*’(x)|.

Then,

Nh
5 5
{b-a} A (b-a)
M,= ■4 A/4 »(18.152)
90 180 2^8QN 4

9
 Numerical Methods 18.107

since Nh - (b - a}l2. If the error tolerance e is prescribed, it is possible to find the number of sub­
intervals required to achieve that accuracy. We get from (18.152)

{b - g)^ M^ < e, or >

(b -
(18.153)
2^^QN^ 2880e

where TV is an integer.
Example 18.41 Evaluate the following integrals using Simpson’s 1 /3rd rule with two and four sub­
intervals. Compare with the exact solution.

(0 J' ' dx r dx
'o 3 + 2x '0 x^ + 2x + 10
Jo

Solution

(i) We have a = 0, b = 1 and/(.v) = 1/(3 + 2a'). For two sub-intervals, we get h = 1/2. The
abscissas are Xq = 0, x, = 1/2 and X2 = 1.0. Hence,

/(0) + 4/W + /(l) 1 1 1

- + 4 - +- = 0.25556.
3 6 3 4 5

For four sub-intervals, we get h = 1/4. The abscissas are -v,, = 0, v, = 1/4, as = 1/2 = 3!A
and X4 = 1. Hence,

h
I= - +
3

1 1 2 2 1 1
± + 4 - + - +2 - +- = 0.25542.
12 3 7 9 4)5

The exact solution is 0.5 In (5/3) 0.25541. The errors in magnitude are respectively
0.00015 and 0.00001.
(ii) We have a = 0, b = 2 and f{x} = l/(x2 + 2x + 10). For two sub-intervals, we get h = 1. The
abscissas are Xq = 0, Xj = 1 and X2 = 2. Hence,

h 1 1 1 1
+ 4/(1) +/(2)] = - - + 4 - +- = 0.15442.
3 3 10 13 18

For four sub-intervals, we get h = 1/2. The abscissas are Xq = 0, x, = 1/2, v, - 1, A3 - 3/2
and X4 = 2. Hence,

/(0) + 4 / i +/ f +2/(l) + /(2)

3 v^yj

10
18.108 Engineering Mathematics

1 1 = 0.15454.
+ 4(0.08889 + 0.06557) + 2(0.07692) + 0.05556
6 10

The exact solution is 0.15455. The errors in magnitude respectively are 0.00013, 0.00001.

Simpson’s 3/8*** rule

We can derive a fonnula by taking n 3, that is, approximating the integrand /(x) by a cubic
polynomial. We have h = {b - a^/i. The abscissas are Xq = a, Xj = Xq + h, X2 - Xq + 2h,
x-i= Xq + ih = b. From (18.140), the rule is given by

f /(x) dx = Ao/(xq) + A,/ (x,) + ^2/ (^2) +

'a
We obtain from (18.141) and (18.142)

h f3. ih ■ h '^h
- f (5 - 1) (5 - 2) (5 - 3)
6 Jo 8 ’

f3
h_^\^s(s-X)(s-3)ds=—9h , h ih
A2 A3 = — 5 (s - 1) (s - 2) ds
6 Jo 8

The rule is given by

•b
\ fix') dx= — + 3/(x,) + 3/(X2) +/(X3)]
(18.154)
•'i/ o •

This rule is called the Simpson s 3/8"’ rule.

The error term is given by

a
X3. (18.155)
80
Composite rule As in the case of Simpson’s 1/3'''* rule, if the length of the interval [a, 6] is large,
then the error expression given in (18.155) becomes meaningless. In this case, we subdivide [a, 6]
into a number of sub-intervals of equal length such that the number of sub-intervals is divisible by
3. Let the number of sub-intervals be iN. Therefore, the number of sub-intervals can be 6, 9, 12 etc.
Then, the number of nodal points are 7, 10, 13 etc. We evaluate each of the integrals by the
Simpson’s 3/8*" rule and simplify.

Let h - {b - a')l{'iN). The abscissas are given by

a = X,0’ x, - Xq + h, X2 = Xq + Ih, X2 = Xq + 'ih, ..., Xyf^ = Xq + = b.
We write.

f {x} dx /(x) dx = /(x) dx + f (x) dx + ... + f (x) dx

'a
'3N-3

11
 Numerical Methods 18.109

3/z
Y [{/(Xo) + 3/(x,) + 3f[xi) +f{x^'}} + {/(X3) + 3/(X4) + 3/(X5) +/(X6)}

3h
= ^\f + 3/(x,) + 3f{xi) + 2f{x^} + 3/(x4) + 3/(x5) + 2f(x^}
Q

+ ... + 2/(x3^^3) + 3f{x^j^_i) + 3/(X3^_i) +/(X3;v)] (18.156)

The error expression becomes
■3
X3, etc. (18.157)
80

3Nh^
Then, 1^3 (A ^)l M4, where M^ = max
80 N1

Remark 14
Simpson’s S/S**^ rule has a disadvantage. From the error expression (18.157), we conclude that the
rule integrates exactly polynomials of degree < 3, which is same as for the Simpson’s 1/3’’*^ rule.
Simpson’s 1/3^** rule requires computation of three function evaluations, while Simpson's 3/8'*’ rule
requires computation of four function evaluations. Further, the error coefficients in the Simpson’s
1/3'"** and 3/8'*’ rules are —1/90 and -3/80 respectively. Therefore, computationally Simpson's 1/3”* rule
is superior to the Simpson’s 3/8‘*’ rule.

Example 18.42 Using Simpson’s 3/8^*’ rule, evaluate the following integrals with 3 and 6
sub-intervals. Compare with the exact solution.

dx ■3 dx
J'
(i) '0 4 + 3x ’
Jo
JJo 1-1-X 2 ■
1

Solution
(i) Three sub-intervals: h = 1/3. The abscissas are 0, 1/3, 2/3, 1. We have the following data.
X 0 1/3 2/3 1
/(X) 0.250000 0.200000 0.166667 0.142857

l>h
A — [<(0) + 3/(1/3) + 3f{2l3}
o
3
= — [0.250000 + 3(0.200000) + 3(0.166667) + 0.142857] = 0.186607.

Six sub-intervals: h = \I6. The abscissas are 0, 1/6, 2/6, 3/6, 4/6,5/6, 1.
We use the previous data values and the following data values.

X 1/6 1/2 5/6

/(x) 0.222222 0.181818 0.153846
Hence,

12
18.110 Engineering Mathematics

[y (0) + + 3/(1/3) + 2f {MT) + 3/(2/3) + 3f {51b) +/(!)]

8

= A [0.250000 + 3 (0.222222) + 3 (0.200000) + 2 (0.181818)

48
+ 3 (0.166667) + 3 (0.153846) + 0.142867] = 0.186544.
The exact solution is
1
1 = I [log 7 -tog4] = 0.186539.
I = - log(4 + 3x)
3 Jo

The magnitudes of errors in the solution are 0.000068 and 0.000005 respectively.
(ii) Three subintervals: h = 1. The points are 0, 1, 2, 3. We have the following data.

X 0 1 2 3
/(X) 1.0 0.5 0.2 0.1
Hence,

A= V [f(0) +3/(1) +3/(2)+/(3)]

o

3
= - [1.0 + 3(0.5) + 3(0.2) + 0.1] = 1.2.
8
Six subintervals:, h = \I2. The points are 0, 1/2, 1, 3/2, 2, 5/2, 3.
We use the previous data values and the following data values.
X 1/2 3/2 5/2
/(X) 0.8 0.307692 0.137931
Hence,

3/2
4= V + 3/(1) + 2/(3/2) +
(S

3
= — [1.0 + 3 (0.8) + 3 (0.5) T 1 + 3 (0.2) + 3 (0.137931) + 0.1] = 1.242971.
16
The exact solution is

1= Jtan tan ' 3 = 1.249046.

The magnitudes of errors in the solution are 0.049046 and 0.006075 respectively.

Boole’s rule

can derive a formula by taking n = 4. We have h = {b - a}!^. The abscissas are

^0 = a, X, = Xq + h, X2 = Xq + Ih, X3 = Xq + 3h, and X4 = Xq + 4A = b.

From (18.140), the rule is given by

13
 Numerical Methods 18.111

f/>
jy (x) dx = V(^o) + + (x^) +

We obtain from (18.142) and (18.142)

= r(5-l)(5-2)(5-3)(5-4)J5 = —, h P•4
64/2
A,=- 5 (5 - 2) (.V - 3) (.V - 4) r/.v =
24 Jo 45 6 Jo 45
■4 ■4
f V (5 - 1) (5 - 3) (5 - 4) ds = —, h
4 JO 15
^3
6
h(^
Jo
64/2
-\){s~2}is-A}ds = ^,
45

A4 = h f\ (5 - 1) (.s- - 2) (5 - 3) ds =
24 Jo 45

The rule is given by

f*I fix} dx=—

2/z
[Ifix^} + 32fix,} + 12/(x2) + 32/(x3) + 7/(x4)]. (18.158)
’’a 4j

This rule is called the Boole’s rule.

The error term is given by

8 Xo<e
Ra {f,x} = - X4. (18.159)
945
The formula integrates exactly polynomials of degree < 5. (Order = 5).

18.6.2.2 Romberg integration

Romberg integration is a powerful tool which uses the method of extrapolation. Generally, to obtain
accurate results, we compute the integral by a Newton-Cotes formula for a number of values of step
lengths, each time reducing the step length. We stop the computation, when the specified accuracy
is obtained, that is, when the magnitude of the difference of successive values of the integral obtained
by reducing values of the step lengths is less than a given accuracy. Convergence may be obtained
after computing the values of the integral using a number of step lengths. In many cases, convergence
may be slow. It may be noted that while computing the value of the integral with a particular step
length, the values of the integral obtained earlier by using larger step lengths were not used. Romberg
integration method is derived by studying the error of the method that is being used. The method uses
the values of the integrals obtained with various step lengths, to refine the solution such that the new
values are of higher order, that is as if the results are obtained by using a higher order method than
the order of the method used. We illustrate the Romberg method for the trapezoidal rule and the
Simpson’s l/3rd rule.

Romberg method for the trapezoidal rule

Let the integral (18.135) be computed by the composite trapezoidal rule. Let I denote the exact value
of the integral and Ij- denote the value obtained by the composite trapezoidal rule. The error, I 17-,

14
18.112 Engineering Mathematics

in the composite trapezoidal rule in computing the integral is given by

I -Ij- = + C2+ €3/2^+ ..., or /= /y + + C2h^ + ••• (18.160)

where c Co, ... are independent of h.

Let I be evaluated using two step lengths h and qh, 0 7 1. Let these values be denoted by
/j-(A) and I-j-^qh}. The error equations become

I~Ir(h} = c^h^ +C2h^ + ... (18.161)

I - iT^qh} = c^q^h^ + C2q'^h^ + ... (18.162)

Eliminating C| from these two equations (neglecting the higher order terms), we obtain

(1 -q^} I- It {qh} + q\ {K) = C2q^h'^ {q^ - 1).

Solving for I, 'WQ obtain

1= - c^q 21h 4 .
(1-^^)

Note that the error on the right hand side is now of order
Neglecting the C>(A'*) error term, we obtain the new approximation to the value of the integral as

I- =
/r(£^WW (18.163)
(I-9")
This computed result is of order, which is higher than the order of the trapezoial rule, which
is of 6>(A'). For q = 1/2, (computations are done with step lengths h and hl2'), the formula (18.163)
simplifies to

/^(A/2)-(l/4)/^(/;)
4'\/r) «
l-(l/4)

4/^(A/2)-/7-(/;)
(18.164)
4-1
For ease in computations, we normally use the sequence of step lengths as A, A/2, hl2^, hl2^, ...
Formula (18.164) gives the first column of extrapolated values, which are of order Repeating
the extrapolation procedure, we obtain the Romberg method as (when the step lengths are
reduced by the factor 2)

4”*
l^(h)« w = 1, 2, ... (18.165)
4'"-!

where The computed result is of order

The extrapolations using three step lengths h, h/2, hl2^, are given in Table 18.5.

15
 Numerical Methods 18.113

Table 18.5. Romberg method for trapezoidal rule.

Step length Value of I: 0{h^) Value of I: 0{h^} Value of I:

h /(A)

Ml I{hl2) {K} =
15

/'’(A/2) =
M{h/4'}-I{hl2}
3
hl4 I{hlA}

Note that the most accurate value is the values at the end of each column.

Romberg method for the Simpson’s l/3rd rule

We can apply the same procedure as in the trapezoidal rule to obtain the Romberg’s integration
procedure for the Simpson’s 1/3rd rule. The error, I 7^, in the composite Simpson’s l/3rd rule in
computing the integral is given by
I— + C2h^ + CjA® + or I= ,6 + c-^h^ + ...
+ c^h^ (18.166)
Following the extrapolation procedure described above, we obtain the Romberg integration
procedure for the composite Simpson’s l/3rd rule as
4"'+' 77 ”'(A/
, 2)-7*7 ')*(A)
4^+1 m = 1, 2y ... (18.167)
-1

where /^“’(/t) = The computed result is of order

The extrapolations using three step lengths h, hl2, hl2^, are given in Table 18.6. Note that the most
accurate value is the value at the end of each column.

Table 18.6. Romberg method for Simpson’s 1/3rd rule.

Step length Value of I: Value of I: Value of 1:

h iw
167(/;/2)-7(/z)
I^\h} =
15

647^'\/;/2)-7^'\/z)
hl2 I{hl2'} 7(2> (h) =
63

167(/t/4)-7(/t/2)
I^\hl2) =
15
h!^

16
18.114 Engineering Mathematics

Example 18.43 Compute the extrapolated values of the integral given in Example 18.40(i) by
Romberg method and the trapezoidal rule. Compare with the exact solution.

Solution We obtain the following extrapolated values as given in Table 18.7,

Table 18.7. Extrapolated values using the trapezoidal rule.

h values values

0.5 /(A) = 0.25833

4/(/i/2)-/(/i)
= 0.25542
3

0.25 / (A/2) = 0.25615

The magnitudes of the errors in the computed and extrapolated values are 0.00292, 0.00074 and
0.00001 respectively.

18.6.2.3 Gauss quadrature rules

The formulas that we have derived in the previous section (Newton-Cotes rules) are formulas based
on uniform mesh spacing. Gaussian integration rules are fonnulas based on non-uniform mesh
spacing. Consider the general intergration rule as

1= f
Ja
hl
w(x)/(x) dx= Xq/ (xq) + Al/(X|)+...+ A„/(x„). (18.168)

Here, wCx) is called the weight function and .v/s are called abscissas of the formula. When the
abscissas are not prescribed in advance and they are also to be determined, then the formulas using
lesser number of abscissas can produce higher order methods compared to the Newton-Cotes
formulas. Gaussian integration rules can be derived when the limits are finite or one of the limits is
infinite or both the limits are infinite. Gaussian integration rules depend on the limits of integration
and the expression for the weight function w(x). Since the formula (18.168) has 2n + 2 parameters
(n + 1 weights and n + 1 abscissas), the formula can be made exact for polynomials of degree
< 2n + 1.
The derivation of the integration rule can be done using the following theorem.

Theorem 18.1 If the abscissas x/s are chosen as the zeros of an orthogonal polynomial, orthogonal
with respect to the weight function ufx) over [a, 6], then the integration rule (18.168) has order
2n + 1. The weights are given by

Xi= w(x) lj(x} dx,

hl

where /,(%) are Lagrange fundamental polynomials. Further, X,- 0.

1. Gauss-Legendre integration rules

Weight function = w(x) = 1. Limits of integration = [- 1, 1],
Abscissas = Zeros of the corresponding Legendre polynomial.

17
 Numerical Methods 18.115

2. Gauss-Chebyshev integration rules

Weight function = h’(x) = 1/-J1 -x^

. Limits of integration = [- 1, 1],
Abscissas — Zeros of the corresponding Chebyshev polynomial.
3. Gauss-Laguerre integration rules
Weight function = w(x) = e~^. Limits of intergration = [0, c»).
Abscissas = Zeros of the corresponding Laguerre polynomial.
4. Gauss-Hermite integration rules
.2
Weight function = vv(x) = e , Limits of integration (- oo, oo).

Abscissas = Zeros of the corresponding Hermite polynomial.

Note that we can also derive the above methods using the method of undetermined parameters.

Order of a method
A method is said to be of order p if it integrates exactly polynomials of degree < p. That is, it
integrates exactly for f (x) = 1, x, x^, ..., x’’.
For the Gauss-Legendre and Gauss-Chebyshev formulas, define

£ H<X) x”*' dx - [;l„(x„'’*') + A,(xr') + ... + A„(xf')].

c=

c y(n + 1)
Then, error in the formula = (^, -1<^ 1.
(« + l)!
Similarly, we define errors in Gauss-Laguerre and Gauss-Hermite formulas.

1. Gauss-Legendre quadrature rules

To derive the rules, we reduce the interval [a, b] to [-1, 1] by using a linear transformation. Write
the transformation as x = pt + q. "Ne, get, a Z? p + q, whose solution is ;? = (Z? - (7)/2,
q = (b a}l2. The linear transformation is given by x = [(Z> - a)t + {b + a')\l2. Hence,

£,
fl fl
f -{{b-a'jt^^b^a}} dt ~
f /(X) dx =
Ja
'a \
— f
y »— 1 Lz
dt.

Without loss of generality, we may write this integral as J f (x) dx.

Gauss-Legendre one point rule

We use the method of undetermined parameters. The one point rule is given by

J {x) dx^ 2Qf{xQ), Ao * 0.

The method has two parameters Aq, Xq. Making the formula exact for /(x) = 1, x, we get

18
18.116 Engineering Mathematics

■1
/(X) = 1: J* dx — 2 — Aq, f (x) = x: J x^ = 0 = Afl Xq.

Since, Aq 0, we get Xq = 0. Therefore, the Gauss-Legendre one point rule is given by

j’]/(x)dx=2/(0). (18.169)

The error term is given by, error = (1/3) - 1 < 1. Since the error term contains /"(<^),
Gauss-Lagendre one point rule integrates exactly polynomials of degree less than or equal to 1.
Therefore, the results obtained from this rule can be compared with the results obtained from the
trapezoidal rule. However, we require two function evaluations in the trapezoidal rule whereas we
need only one function evaluation in the Gauss-Legendre one point rule. The corresponding
Legendre polynomial is P/x) = x, whose zero is x = 0.

Gauss-Legendre two point rule

The two point rule is given by

J J /(x) dx= (x^ + Al/(xi), Ao ?!: 0, Aj 0, and x,, x,.

The method has four parameters Xq, X], Xq, Xj. Making the formula exact for/(x) = 1, x, x^, x^, we get
•1 ■1
/(X) = 1: die = 2 = An0 + A,■1’ f(x)=x\ J X dtc = 0 = Ao Xq + A, Xj.

pl■1 2 ■t •
/(x) = x^: x^tZx = - = AgXg -F Ajx^, /(x) = X.3.: J'-1
j^dx = 0 = Ao X
q + a, xf.
The solution of these equations is Ag = A] = 1, and Xq = ± (l/->/3j = -Xi.

Therefore, the Gauss-Legendre two point rule i^ given by

+ • (18.170)
3
The error term is given by, error = (1/135) 1 < (^ < 1. Since the error term contains/
Gauss-Legendre two point rule integrates exactly polynomials of degree less than or equal to 3.
Therefore, the results obtained from this rule can be compared with the results obtained from the
Simpson’s l/3rd rule. However, we require three function evaluations in the Simpson’s l/3rd rule
whereas we need only two function evaluations in the Gauss-Legendre two point rule. The
corresponding Legendre polynomial is PjCx) = (3x" - l)/2, whose zeros arc x L I//3.

Gauss-Legendre three point rule

The three point rule is given by
'1
J-i 1
^\f (xi) + ^2 / (xj), Ao 0, Aj 0, A2 0, and Xq Xj x^.

19
 Num erical Methods 18.117

The method has six unknowns £, A,, Ao, Xq, x,,1’ .x^. Making the formula exact for / (x) = 1, x,
x^, x‘*, x^, we get the Gauss-Legendre three point rule as

f/w dx = —9 [5/ (- TsTI)+8/(0) + 5/(7375)]. (18.171)

The error term is given by, error = (1/15750) (^), - 1 < 1. Since the error term contains
j Gauss-Legendre three point rule integrates exactly polynomials of degree less than or equal
to 5. Further, the error coefficient is very small (1/15750). Therefore, the results obtained from this
rule are very accurate. The corresponding Legendre polynomial is ^(x) = x {5x^ - 3)/2, whose zeros
are x = 0, ± 73/5.

Example 18.44 Evaluate the following integrals using Gauss-Legendre one point, two point and
three point rules. Compare with the exact solution.
.2 dx 0. dx
r
I, 1 + X .2 .
Jo'0 x^ -l- 2x + 10

Solution (i) The transformation from [1, 2] to [- 1, 1] is given by

X = {[(/? -a) t + (b + fl)]/2} = (f + 3'}I2. We have dx = dt/2.

.2 dx 1 4dt >1 2dt

Hence,
I, 1 + x 2 [
J-i 4-1-(Z-f3) 1
’-I 4 + (t + 3)
2
dt.

2
One point rule: I, 2f{Q} = 2 = 0.307692.
13

£
Two point rule: h-f +f [3 = 0.202650 + 0.119066 = 0.321716.
[3

Three point rule: [5/- 7575) + 8/(0) + 5/(7775)]

= 1 [5(0.223403) + 8(0.153846) + 5(0.109604)] = 0.321756.
9

Exact solution: /= tan-'(2) - {k/A} = 0.321751.

(ii) The transformation from [0, 2] to [- 1, 1] is given by

X = {[(6 - d)t -F (Z? + a)]/2} = t -t- 1. We have dx = dt.

>2 • 1 J/
dx
Hence, £/(z)t/t.
'0 x^ +2x + 10 U (t + l)2 +2(t + l) + 10

One point rule: /] = 2/(0) = (2/13) = 0.153846.

=/(-1/73 ) + /(l/TI) = 0.090712 + 0.063927 = 0.154639.

Two point rule:

20
18.122 Engineering Mathematics

■2
6. Evaluate cos x^/x using the trapezoidal rule with (i) h = 1, (ii) h= \/2. Compare with the exact
I
solution. Find the maximum error in each case.
>1

7. Find the minimum number of intervals required to evaluate In (1 + x) dx, using trapezoidal rule
'o

with an accuracy of 10 ’ .
8. Using the following data
X 1 2 3 4
/(X) 0.3679 0.1353 0.0498 0.0183

and the trapezoidal rule with n = 1 and « = 3, determine an approximate value of /(x)i/x.
JiI

>!cI2
9. Evaluate the integral e cos X dx, using trapezoidal rule with h = nl2 and /i = kI4.
'o

- x^ cos X dx, using the trapezoidal rule with h

10. Obtain the approximate value of 1,
1/2 and 1/4.
,2
f■2 X
11. Evaluate dx using the Simpson’s l/3rd rule with two and four sub-intervals. Compare with
J' 14-X.3

the exact solution.

•2
12. Evaluate J yj\ + 4x~ sin X dx using the Simpson’s l/3rd rule with A = 1/2 and 1/4.

>1
13. Find the minimum number of intervals required to evaluate In (1 + x}dx using the Simpson’s l/3rd
'o
rule with an accuracy of 10“®.
fS•;r/2
14. Determine the maximum error in evaluating e ' cos X dx, using the Simpson’s l/3rd rule with
Jo
h =
15. Using the Simpson’s 3/8"’ rule, evaluate the following integrals with 4 and 7 nodal points. Compare
with the exact solution.

f. dx (ii) r dx (x -t-1) t/x

® J, 'I 5 + 3x Jo 14-X'.2 •
Jn
Jo x^ -t- 2x -l- 2
16. Using the Gauss-Legendre one point, two point and three point rules, evaluate the following integrals.
Compare with the exact solution wherever possible.

(x+l)’ dx
(i) I
Jo 1 4- (X -F 1)'A
dx. ® £ 25-t-x'2 '

... dx
(ni) ----------
Jo 5 + 4x
(iv) pJo'0 x^ -F dx2x4-10

25