Numerical and Statistical Methods - CRC
Numerical and Statistical Methods - CRC
METHODS
Bharathidasan University
Centre for Distance and Online Education
Chairman:
Dr. M. Selvam
Vice-Chancellor
Bharathidasan University
Tiruchirappalli-620 024
Tamil Nadu
Co-Chairman:
Dr. G. Gopinath
Registrar
Bharathidasan University
Tiruchirappalli-620 024
Tamil Nadu
Course Co-Ordinator:
Dr. A. Edward William Benjamin
Director-Centre for Distance and Online Education
Bharathidasan University
Tiruchirappalli-620 024
Tamil Nadu
The Syllabus is Revised from 2021-22 onwards
Review er
Dr.T.Jai Sankar, Assistant Professor & Head, Department of Statistics, Bharathidasan University, Khajamalai Campus,
Tiruchirappalli-620023
Author
Kalika Patrai, Senior Lecturer in Galgotia Institute of Management and Technology, Greater Noida
All rights reserved. No part of this publication which is material protected by this copyright notice
may be reproduced or transmitted or utilized or stored in any form or by any means now known or
hereinafter invented, electronic, digital or mechanical, including photocopying, scanning, recording
or by any information storage or retrieval system, without prior written permission from the Publisher.
Information contained in this book has been published by VIKAS® Publishing House Pvt. Ltd. and has
been obtained by its Authors from sources believed to be reliable and are correct to the best of their
knowledge. However, the Publisher, its Authors shall in no event be liable for any errors, omissions
or damages arising out of use of this information and specifically disclaim any implied warranties or
merchantability or fitness for any particular use.
Numerical methods: Floating point arithmetic; Basic Unit 1: Errors and Floating Point
concept of floating point numbers systems, implications of Arithmetic
(Pages 3–28);
finite precision, illustrations of errors due to round off.
Interpolation Finite difference calculus, polynomial Unit 2: Interpolation
interpolation, approximation uniform, discrete least square, (Pages 29–84);
polynomial, fourier. Unit 3: Numerical Integration and
Numerical integration and Differentiation interpolatory Differential Interpolation
(Pages 85–123)
numerical integration; Numerical differentiation.
Solution of non-linear: Bisection, fixed point iteration, Unit 4: Solution of Algebraic and
Newton’s Rephsons methods. Transcendental Equations
(Pages 125–163);
Solution of Ordinary Differential Equation:Taylor series
method, Runge-Kutta method, Euler method. Unit 5: Numerical Solution to
Ordinary Differential
Equations
(Pages 165–194)
Self-Instructional Material 1
Introduction and engineering problem accuracy of results, approximations and errors introduced
carry great importance and statistics deals primarily with numerical data gathered
from surveys or collected using various statistical methods. Its objective is to
summarize such data, so that the summary can give us an indication of certain
NOTES characteristics of a population or phenomenon that we wish to study. This book
Scientific and Statistical Computing has been prepared for students in a
professional manner providing helpful and relevant material in a lucid, self-
explanatory and simple language to help them understand the concepts easily.
Relevant statistical tables have been given in the appendix for students refer
to it for solving statistical problems with accuracy. The concepts have been analyzed
in a logical format, beginning with an overview that helps readers to easily grasp
the concept. This is followed by explanations and solved examples. Questions
and examples in the ‘Check Your Progress’ and ‘Questions and Exercises’ sections
will further help in understanding and recapitulation of the subject.
2 Self-Instructional Material
Errors and Floating Point
POINT ARITHMETIC
NOTES
Structure
1.0 Introduction
1.1 Unit Objectives
1.2 Approximate Numbers
1.2.1 Precision
1.2.2 Rounding Off
1.2.3 Defining Precision in Hardware
1.2.4 Compiler Evaluation of Mixed Precision
1.3 Errors
1.3.1 Error in the Approximation of a Function
1.3.2 Error in a Series Approximation
1.3.3 Errors in Numerical Computations
1.4 Floating Point Representation of Numbers
1.4.1 Arithmetic Operations with Normalized Floating Point Numbers
1.4.2 Drawbacks of Floating Point Representation
1.5 Summary
1.6 Key Terms
1.7 Answers To ‘Check Your Progress’
1.8 Questions and Exercises
1.9 Further Reading
1.0 INTRODUCTION
Floating point arithmetic describes a system for representing numbers that would
be too large or too small to be represented as integers.
In this unit, you will learn that there are many situations where analytical
methods are unable to produce the desired results. These limitations of analytical
methods in practical applications have led mathematicians to develop certain
numerical methods.
Hence, the aim of numerical techniques is to provide constructive methods
for obtaining answers to such problems in a numerical form. However, in most
applications that make use of numerical techniques, the input data is not always
exact. It comes from some measurement or the other and hence contains some
error. The input data refers to the approximations, precise upto a certain number
of significant digits. These approximations, as a result, introduce error in the final
computed result.
The error in the final result may be due to an error in the initial data or in the
method or both. Therefore, the use of approximate numbers and/or approximate
methods gives us a solution that is not exact but approximate. So, the main objective
Self-Instructional Material 3
Errors and Floating Point is to obtain an approximate solution of the problem which would be closest to the
Arithmetic
exact solution.
This unit will discuss the concept of approximate numbers, significant digits,
floating point numbers and operations performed on them. Also, you will learn
NOTES about different types of errors, their sources and propagation of errors in numerical
calculations.
4 Self-Instructional Material
• Leading and trailing zeros (unless a decimal point is present) Errors and Floating Point
Arithmetic
where they serve merely as placeholders to indicate the scale of the
number.
• Spurious digits introduced, for example, by calculations carried out to
greater accuracy than that of the original data or measurements NOTES
reported to a greater precision than the equipment supports.
Significant digits are used to express a number. For example, the numbers
3456, 34.56, 0.3456 have four significant digits each, whereas the numbers
0.00345, 0.003456 and 0.00034567 have three, four and five significant digits,
respectively, since zeros are used only to fix the position of the decimal point.
Similarly the number 0.34560 has only five significant digits. So,you can
say that ‘0’ is a significant digit except when it is used to fix decimal points or to fill
the places of unknown or discarded digits.
The best way to identify the significant digits in a given number is to write
the number in the scientific notation with the first digit being non-zero. The numbers
0.003456, 345600 and 345.603 can be expressed by using scientific notation as
0.3456 × 10–2, 0.3456 × 106 and 0.345603 × 103, respectively. The first part is
known as Mantissa and the second part, i.e., 10n is known as exponent.
• 3 4 5 6 ×10 –2
Mantissa Exponent
The digits in mantissa part are the only significant digits in the number.
1.2.2 Rounding Off
22
Consider a number say = 3.142857143.... In practice, you need to limit such
7
numbers to a manageable number of digits such as 3.1428 or 3.143. This process
of dropping the unwanted digits is called rounding off.
There are some guidelines or rules for rounding off a number to n-significant
digits. These guidelines are as follows:
1. Discard the digits to the right of the nth digit.
2. If the discarded number is:
(i) Less than half a unit in the nth place, leave the nth digit as it is
(i.e., unchanged). For example, the number 6.4326 when rounded-
off to three significant digits is 6.43 as the discarded number is
less than half a unit in the nth place.
(ii) Greater than half a unit in the nth place, increase the nth digit by
unity (i.e., add 1 to the nth digit).
For example,
6.43267 ~ 6.4327
3.2589 ~ 3.26
Self-Instructional Material 5
Errors and Floating Point (iii) Exactly half a unit in the nth place, increase the nth digit by unity if
Arithmetic
it is odd otherwise leave it unchanged.
For example,
6.4354 ~ 6.43
NOTES
5.8235 ~ 5.823
Thus, the numbers rounded off to n significant digits are said to be correct
or exact to n significant digits.
1.2.3 Defining Precision in Hardware
In the field of computing, floating point denotes a way to represent numbers that
are too small or too large when we represent these as integers. When using
computers for representing precision such numbers are, generally represented
using some fixed number of significant digits and these are scaled by use of an
exponent. Base used for scaling is usually 2, 10 or 16. Exact representation of a
typical number has the form:
Significant digits × baseexponent
Floating point is so called since the radix point can ‘float’, which means that
it may be placed anywhere in relation to significant digits for the number and such
position is shown separately in some kind of internal representation in the hardware.
Thus, representation of floating-point is a computer implementation of scientific
notation. There are different systems for representing floating-point in computers
and most common representation has been defined by standard IEEE 754.
With representation of integral numbers also, maximum number and minimum
number is limited by capability of computer hardware in storing such numbers. But
use of floating-point can support wider range of numerical values than that of
fixed-point representation. If we use a fixed-point representation having six decimal
digits and we assume positioning of decimal point after fourth digit, we may represent
numbers 1234.56, 3765.43 and like that. But while representing a floating-point
with six decimal digits we may take it as 1.23456, 12345.6, 0.0000123456,
123456000000000, and like that. Computer hardware needs more storage for
representing floating point numbers. Speed for floating-point operations is used to
measure performance of computers and it is measured as FLOPS.
By making radix point adjustable, notation for floating-point permits
calculations of different magnitudes with the use of some fixed number of digits
and this can maintain good precision. For example, a decimal floating-point system
using three digits, as decimal multiplication may be written as
0.11 × 0.11 = 0.0121 and this can be written as:
(1.1 × 10 ) × (1.1 × 10–1) = (1.21 × 10–2).
–1
But when a fixed point system having decimal point put on the left side, it
becomes 0.110 × 0.110 = 0.012.
6 Self-Instructional Material
Thus, due to the limit put for representing digits, one digit in the result has Errors and Floating Point
Arithmetic
been lost. This is so as decimal point and digits do not float relative to each other.
Range of floating-point numbers is dependent on number of digits or bits
used to represent significant digits which is known as significand, and an exponent.
NOTES
In a typical computer system, that uses 64 bits for representing a ‘double precision’
number, the binary floating-point number uses 53 bits with one bit for sign and 11
bits for exponent a positive floating-point number has a range of approximately
10–308 to 10308 and 308 is nearly 1023 × log10(2), as range of exponent is form
–1022 to 1023. This format has complete range of –10308 to +10308 as per IEEE
754.
In the field of computer science, a technique known as arbitrary-precision
arithmetic is used in which calculations are made using where available memory
of the computer system limits digits of precision. This has contrasted that of faster
fixed-precision arithmetic available in most hardware of ALU (Arithmetic and
Logic Unit) having a typical range of 6 to 16 decimal digits.
Standardization has been done by IEEE for computer representation of
binary floating-point numbers. Most modern machines follow this standard. IBM
mainframe is an exception to this and has its own format.
Precision defined by IEEE 754 for floating point representation are:
• Half precision: 16-bit.
• Single precision: 32-bit binary32 and decimal32
• Double precision: 64-bit, binary64 and decimal64
• Quadruple precision: 128-bit: binary128 and decimal128
Number of normalized floating point numbers for a system is given by
F(B, P, L, U), where B denotes the number base, P, the precision of numbers, L,
the smallest exponent and U, the largest exponent that can be represented and is
given as:
2 * (B – 1) * B^(P – 1) * (U – L + 1).
UFL (Underflow Level) = B^L is the least positive normalized floating-
point number and has a 1 as leading digit and 0 for remaining digits in the mantissa
with least possible value for exponent.
Also, OFL (Overflow Level) = B^(U + 1) * (1 – B^(–P)) which is the
maximum floating point number having B - 1 as value for every digit of mantissa
and maximum possible value for exponent. In addition to these, actual values that
are represented lie between –UFL and UFL.
A number is represented by specifying some way for storage of a number
encoded as a string of digits. Logically, a floating-point number has:
• A digit with sign for some given length to a given base also known as radix.
This is called significand, or mantissa or coefficient. Length of significand
gives the precision by which numbers are represented.
Self-Instructional Material 7
Errors and Floating Point • A signed integer exponent called scale or characteristic that modifies the
Arithmetic
number’s magnitude.
The number is represented as: significant digits × baseexponent.
NOTES ‘Machine precision’ shows a quantity characterizing accuracy of a floating
point system. It is also called machine epsilon and generally denoted as
Emach whose value depends on rounding off method and when rounding to zero,
Emach = B^(1– P) but when rounding is done to nearest digits, Emach = (1/2)*B^(1
– P)
Such a concept is important as it bounds relative error in showing a non-
zero real number x within normalized range of a floating point system, such as:
|(f l(x) – x)/x| < = Emach
Storage Layout
There are three basic parts of IEEE floating point numbers; Sign, Exponent, and
Mantissa. Mantissa is fraction with implicit leading digit. If base of the exponent is
2, it is implicit and this does not require storage.
Table below shows layout for single precision using 32 bits and double
precision having 64-bit floating-point values. Numbers of bits in each field have
also been shown. Square brackets have been used to denote bit ranges:
The Exponent
This field has to denote positive as well as negative exponents. Hence, a bias is
put in actual exponent for storing exponent. In case of IEEE single-precision floats,
it is 127. If an exponent is zero it tells that 127 is stored for exponent field. If
stored value is 200, it shows exponent as 200-127 = 73. The exponents of -127,
in which all bits are 0s and +128, in which all bits are 1s and as such these are kept
reserved and stand as special numbers.
In case of double precision, exponent field contains 11 bits having a bias of
1023.
8 Self-Instructional Material
The Mantissa Errors and Floating Point
Arithmetic
Mantissa is known by another name significand, and it shows precision bits in the
number. It contains an implicit leading bit and fraction bits.
1.2.4 Compiler Evaluation of Mixed Precision NOTES
Hardware puts a limit on accuracy depending on how much bits are used for
floating point representations. Different programming languages define different
ways to represent floating points in computer system. A program written in
programming language is compiled to create an object code that is a binary code
accepted by the computer machine.
Precision, in science and engineering, is defined as degree of agreement
amongst a group of individual measurements or results. Precision in computation
and accuracy that is found in the final result has a complicated relationship. Non
linear relations, in most of the cases lead to a situation in which improvement in
accuracy is different for different part of an algorithm. Different methods of numerical
analysis gives knowledge on such relations that can be used by reducing
computational precision in areas that are less sensitive and increasing in areas that
are very sensitive and use logarithm accordingly. Such an approach is known as
mixed precision. Thus, in this approach, different levels of precision are set for
different portions of an algorithm and compiler evaluates all.
Accuracy is given by the degree to which there is an agreement for a
measured value with that of actual value. If actual value is known, there is no need
to measure it as such and you may use it for determining accuracy of the measuring
instrument or tool. If actual value is not known, accuracy can be inferred for the
precision of the measurement. For example, if calculation is made for knowing
expenditure in terms of Rupee and paisa, two decimal digits are enough and it is
wise to calculate up to three decimal digits and then rounding to two digits in the
final result.
Hardware in a computer system gives a limit to which number can be
represented but its full potential is utilized by the software programs written in
some programming language. Support for number representation is provided by
hardware. While computing, precision is used differently. One approach is the
number of bits for representing mantissa of the floating point number. A single-
precision floating point number as per IEEE has 24 bits of precision that is equivalent
to 7 decimal digits. A double-precision has 53 bits of precision equivalent to 15
decimal digits.
Precision is also used to tell as to how many bits amongst these are correct.
In numerical analysis computation is formulated to improve precision and accuracy.
The order in which computations are made has bearing on accuracy. Accordingly,
in a programming language somewhere short integer is used whereas at some
Self-Instructional Material 9
Errors and Floating Point places long integer is used. Similar is the case of decimal numbers. If one is writing
Arithmetic
a program to calculate factorial, he has to see the limit to which a number can be
evaluated without overflow. Since compiler is capable of changing the order of
computations, it can affect the result.
NOTES
Floating point operations are slower. Multiplication operation takes more
time than addition. Square root or division takes ten times longer than addition.
Thus, to improve speed it is wise to reduce multiplication operation by replacing
with addition wherever possible and replacing division by multiplication. Current
microprocessors, Intel (SSE), AMD and PowerPC have built-in instructions for
approximating a single-precision inverse square root. In these cases, instructions
give a low precision approximation that is nearly equal to half the number of bits.
In iterations, used in Newton-Raphson method of numerical computations, results
can be improved to full precision.
Cray 1 implemented a method of reciprocal approximate in place of divide
instruction that required compiler to generate iterations in Newton-Raphson method to
give result with full precision. There are round-off errors in any computer floating point
arithmetic and careful steps has to be followed to avoid results that may be misleading.
There are many formats for floating point arithmetic from base 2 to 16
without or with rounding, with variation in its dynamic range. IEEE in eighties
made an effort for standardizing floating point formats that is now used by all
mainstream processors of today. This provided compatibility such that answers
computed with one type of processor; Intel, IBM or AMD will give the same
result on other brands of microprocessors such as Motorola, SPARC or MIPS.
But in practice, changing the processor changes computation and results produced
by different processors may be nearby but not exactly the same. A compiler
command ‘!PRECISION’sets the value of decimal-places to determine minimum
number when calculating and giving results.
1.3 ERRORS
10 Self-Instructional Material
or, E = | Et – Ea | Errors and Floating Point
Arithmetic
Where, Et is true value or exact value and Ea is approximate value.
The following are the different types of errors that occur during numerical
computations:
NOTES
• Inherent errors
Errors that exist in a problem before it is solved are called inherent errors. These
errors occur due to incorrect measurements or observations that may be due to
the limitations of the measuring instrument such as mathematical tables, calculators
or the digital computer. These errors can be minimized by taking better data or by
using high precision computing aids.
• Truncation errors
The errors which occur when some digits from the number are discarded are
known as truncation errors. Truncation errors occur in two situations:
(i) When the numbers are represented in a normalized floating point form.
(ii) During the conversion of the number from one system to another.
While representing a number in normalized floating point form the error
occurs because only a few digits of the number are accommodated by the Mantissa
part. For example, the number 0.003456789 takes the form of 0.3456 × 10–2 in
its own hypothetical computer during the normalized floating point form.
The error occurs when a number, say 13.1 is represented in a binary form
as 1101.001100110011 ... . It has a repeating fraction and due to this repetition
the conversion is terminated after some digits and hence truncation error is
introduced.
• Round-off errors
These errors occur during the process of rounding off of a number. These are
unavoidable errors due to the limitations of computing. However, these errors can
be reduced by:
(i) Changing the calculation procedure so as to avoid subtraction of nearly
equal numbers.
(ii) Retaining at least one more significant digit at each step than that given in the
data and rounding off at the last step.
• Absolute error
Let Xt be the true value and Xa be the approximate value, then the positive difference
between Xt and Xa, i.e., | Xt – Xa | is known as absolute error, denoted by Ea.
So, E a = | Xt – Xa |
Self-Instructional Material 11
Errors and Floating Point
Arithmetic • Relative error
The relative error denoted by Er is defined as:
Absolute error | X t − X a |
NOTES Er = =
True value Xt
• Percentage error
The percentage error denoted by Ep is defined as:
E p = Relative error × 100
| Xt − Xa |
Ep = × 100
Xt
The relative and percentage errors are independent of the units used while
the absolute error is expressed in terms of these units.
1
Notes: 1. If a number is correct to n decimal places, then the error in it is × 10− n .
2
2. If a number is correct to n-significant digits with the first significant digit
being k, then the relative error Er < 1/(k × 10n – 1).
Example 1.1 If the number 852.47 is correct to five significant digits, then what
will be the relative error?
Solution: Here the first significant digit is k = 8 and total number of significant
digits is n = 5.
1
Then, absolute error = × 10 − m
2
where, m stands for the number that is correct to m decimal places.
1
Here, Ea = × 10−2 = 0.005
2
Ea 0.005 5
Thus, Er ≤ = =
852.47 852.47 852470
1
=
2 × 85247
One can observe that
1 1 1
Er = < =
2 × 85247 2 × 80000 2 × 8 × 104
1 1
< , i.e., k × 10n − 1
8 ×104
12 Self-Instructional Material
Example 1.2 Round off the following numbers correct to four significant Errors and Floating Point
Arithmetic
digits:
0.00032217, 35.46735 and 18.265101 and compute Ea, Er and Ep in each case.
Solution: NOTES
(i) The number rounded off to four significant digits = 0.0003222
So, E a = | 0.00032217 – 0.0003222 | = 0.00000003 = 0.3 × 10–7
0.3 ×10−7
Er = = 0.000093119 = 0.93119 × 10–4
0.00032217
and E p = Er × 100 = 0.0093119 = 0.93119 × 10–2
(ii) The number rounded off to four significant digits = 35.47000
Then, E a = | 35.46735 – 35.47000 | = 0.00265
0.00265
Er = = 0.74717 × 10–4
35.46735
and E p = Er × 100 = 0.74717 × 10–2
(iii) The number is rounded off to four significant digits = 18.270000
Then, E a = | 18.265101 – 18.270000 | = 0.004899
0.004899
Er = = 0.00026822 = 0.26822 × 10–3
18.265101
and E p = Er × 100 = 0.26822 × 10–3 × 100 = 0.26822 × 10–1
1.3.1 Error in the Approximation of a Function
Let y = f(x1, x2), be a function of two variables x1 and x2.
Then, if δx1 and δx2 are the errors in x1 and x2, respectively, the error in δy
in y is given as:
y + δy = f (x1 + δx1, x2 + δx2) ...(1.1)
Expanding the right hand side of Equation (1.1) by using Taylor’s series,
you have
∂f ∂f
y + δy = f (x1 + x2) + . δx1 + . δx2 + ...
∂x1 ∂x2
If the errors δx1 and δx2 are very small, you can neglect the higher order
terms of δx1 and δx2 then,
∂f ∂f
y + δy = f (x1, x2) + . δx1 + . δx2 ...(1.2)
∂x1 ∂x2
Self-Instructional Material 13
Errors and Floating Point
Arithmetic ∂f ∂f
then, δy = . δx1 + . δx2 approximately..
∂x1 ∂x2
NOTES In general, the error δy in the function y = f (x1, x2, ..., xn) corresponding to
the errors δxi (i = 1, 2, ..., n) in xi is given by,
∂f ∂f ∂f
δy ∼ . δx1 + . δx2 + ... + . δxn
∂x1 ∂x2 ∂xn
δy ∂f δ x1 ∂f δ x2 ∂f δ xn
Er = = . + . + ... + .
y ∂x1 y ∂x2 y ∂xn y
( x − a)n n
where, Rn(x) = f (θ), a < θ < x
n!
If the series is convergent, Rn(x) → 0 as n → ∞ and hence if f(x) is
approximated by the first n terms of this series, then the remainder term Rn(x) of
the series gives the maximum error. If you have already received the accuracy that
is required in series approximations you can find out the number of terms, i.e., n to
obtain the desired accuracy.
Example 1.3 Find the number of terms of the exponential series such that their
sum gives the value of ex correct to six decimal places at x = 1.
x 2 x3 x n −1
Solution: ex = 1 + x + + + ... + + Rn ( x)
2! 3! (n − 1)!
xn θ
where, Rn(x) = e ,0<θ< x
n!
xn x
Maximum absolute error at θ = x is e
n!
14 Self-Instructional Material
Errors and Floating Point
xn Arithmetic
and, Maximum relative error =
n!
1
Hence, (er)max. at x = 1 is =
n! NOTES
For accuracy upto six decimal places at x = 1, you have
1 1
< × 10−6 , i.e., n! > 2 × 106
n! 2
This gives n = 10.
Hence, you need ten terms of the series to ensure that its sum is correct to
six decimal places.
Example 1.4 The function f(x) = tan–1x can be expanded as,
x3 x5 x2n − 1
tan–1x = x − + − ... + (−1)n − 1.
3 5 2n − 1
Find n such that the series determines tan–1(1) correct upto eight
significant digits.
x 2n + 1
Solution: If you retain n terms, then (n + 1)th term= (–1)n .
2n + 1
(−1)n
For x = 1, (n + 1)th term =
2n + 1
For the determination of tan–1 (1) correct up to eight significant digits accuracy,
(−1) n 1 −8
< × 10
2n + 1 2
⇒ 2n + 1 > 2 × 108
This gives n = 108.
Example 1.5 Given a value of x = 3.51 with an error δx = 0.001, estimate the
resulting error in the function y = x3.
Solution: You have y = x3
δy = 3x2 . δx
= 3(3.51)2 . (0.001) = 0.03696
Since, f (3.51) = 43.2436
Let us assume that f (3.51) = 43.2436 ± 0.03696
In other words, you can say that the value of the function y = x3 at x = 3.51
lies between (43.2436 + 0.03696) = 43.28056 and (43.2434 – 0.03696)
= 43.20664 if x has an error of 0.001.
Self-Instructional Material 15
Errors and Floating Point
Arithmetic 8xy 2
Example 1.6 Let R = 3 and ∆x = ∆y = ∆z = 0.001. Find the maximum
z
relative error in R at x = y = z = 1.
NOTES Solution: Given that ∆x = ∆y = ∆z = 0.001 then the relative error in R is given as,
∂R ∂R ∂R
δR = · δx + · δy + · δz
∂x ∂y ∂z
8 y2 16 xy 24 xy 2
= 3 · δx + 3 · δy − · δz
z z z4
Since the errors δx, δy and δz (i.e. ∆x, δy, ∆z) may be positive or negative,
you take the absolute values of the terms on the right hand side which gives
8 y2 16 xy 24 xy 2
(δR)max = δx + δy + δz
z3 z3 z4
= 8(0.001) + 16(0.001) + 24(0.001)
= 0.048
Hence, the maximum relative error is given as,
δR 0.048
= = 0.006
R 8
1.3.3 Errors in Numerical Computations
16 Self-Instructional Material
Notes: Let us follow the given procedure for adding up – numbers of different Errors and Floating Point
Arithmetic
absolute accuracies:
1. Isolate the number with the greatest absolute error.
2. Round off all other numbers, retaining in them one digit more than in the NOTES
isolated number.
3. Add up.
4. Round off the sum by discarding one digit.
∆X 1 ∂X 1 ∂X 1 ∂X
Now, = ∆x1 + ∆x2 + ... + ∆xn
X X ∂x1 X ∂x2 X ∂xn
1 ∂X x . x ... xn 1
Now, = 2 3 =
X ∂x1 x1. x2 . x3...xn x1
1 ∂X x .x ... xn 1
= 1 3 =
X ∂x2 x1 . x2 . x3...xn x2
M M
1 ∂X 1
=
X ∂xn xn
Self-Instructional Material 17
Errors and Floating Point
∆X ∆x ∆x ∆x
= 1 + 2 + ... + n
Arithmetic
∴
X x1 x2 xn
∴ The relative and absolute errors are given by,
NOTES
∆X ∆ x1 ∆ x2 ∆ xn
Maximum relative error = ≤ + + ... +
X x1 x2 x3
∆X
Maximum absolute error = X
X
∆X
= . (x1x2x3 ... xn)
X
x1
Let, X =
x2
∆X 1 ∂X 1 ∂X
∴ = ∆x1 + . ∆x2
X X ∂x1 X ∂x2
∆X ∆x1 ∆x2
∴ ≤ + which is relative error..
X x1 x2
∆X
Absolute error = | ∆X | ≤ .X
X
Error in evaluating xk
X = xk, where k is an integer or fraction
dX
∆X = ∆x = kx k − 1. ∆x
dx
∆X ∆x
= k.
X x
∆x
The relative error in evaluating xk = k .
x
18 Self-Instructional Material
Errors and Floating Point
Arithmetic
CHECK YOUR PROGRESS
4. How are errors introduced in numerical calculations?
5. What is the result when two numbers are added? NOTES
6. Name the two formats in which real numbers can be written.
Self-Instructional Material 19
Errors and Floating Point For positive numbers mantissa ranges from 0.1 to 1.0 and for negative
Arithmetic
numbers it ranges from –1.0 to –0.1. In other words, mantissa of a normalized
floating point number must satisfy the following condition:
0.1 ≤ | Mantissa | ≤ 1.0
NOTES
For example, the number 0.3245 × 105 is represented in this form as: 0.3245
E5 where E5 is used to represent 105. Thus, here the mantissa is 0.3245 and
exponent is 5. This number is stored in hypothetical computer as shown in Figure
1.2.
Sign of → + + ← Sign of exponent
Mantissa • 3 2 4 5 0 5
↑ Mantissa | Exponent
Implied decimal point
Figure 1.2 A Memory Locating Storing Number
20 Self-Instructional Material
(ii) Here exponents are not equal. The operand with the larger exponent Errors and Floating Point
Arithmetic
is kept as it is.
0.5567 E 7
+ 0.0045 E 7 | 0.4546 E 5 = 0.0045 E 7 NOTES
0.5612 E 7
(iii) The addition will be as follows:
0.5567 E 7
+ 0.0000 E 7 | Q 0.4546 E 3 = 0.0000 E 7
0.5567 E 7
(iv) 0.6434 E 3
0.4845 E 3
1.1279 E 3 (not normalized)
Here exponents are equal but when the mantissas are added, the sum
is 1.1279 E 3. As the mantissa has five digits and is > 1, it is shifted
right one place before it is stored.
Hence, Sum= 0.1127 E 4
(v) 0.6434 E 99
0.4845 E 99
1.1279 E 99
Here, again the sum of the mantissas exceeds 1. The mantissa is shifted
right and the exponent is increased by 1, resulting in a value of 100 for
the exponent. The exponent part cannot store more than two digits.
This condition is called an overflow condition and the arithmetic unit
will intimate an error condition.
Example 1.8 Subtract the following floating point numbers:
(i) 0.9432 E – 4 from 0.5452 E – 3 (ii) 0.5424 E 3 from 0.5452 E 3
(iii) 0.5424 E – 99 from 0.5452 E – 99.
Solution: (i) 0.5452 E – 3
– 0.0943 E – 3
0.4509 E – 3
(ii) 0.5452 E 3
– 0.5424 E 3
0.0028 E 3
In a normalized floating point, the mantissa is ≥ 0.1
Hence, the result is 0.28 E 1
Self-Instructional Material 21
Errors and Floating Point (iii) 0.5452 E – 99
Arithmetic
– 0.5424 E – 99
0.0028 E – 99
NOTES For normalization, the mantissa is shifted left and the exponent is reduced
by 2. The exponent would thus become – 100 with the first left shift, which
cannot be accommodated in the exponent part of the number.
This condition is called an underflow condition and the arithmetic unit will
signal an error condition.
Notes: 1. If the result of an arithmetic operation gives a number smaller than
0.1000 E – 99 then it is called an underflow condition. Similarly, any
result greater than 0.9999 E 99 leads to an overflow condition.
2. Mantissa of the sum (before normalization) can be a maximum of
+1.9999,where at the most you need to shift the decimal point to the
left by one position in order to normalize it. As a result of this the
exponent of the sum is incremented by 1.
Example 1.9 In normalized floating point mode, carry out the following
mathematical operations:
(i) (0.4546 E 3) + (0.5454 E 8) (ii) (0.9432 E – 4) – (0.6353 E – 5).
Solution: (i) 0.5454 E 8
+ 0.0000 E 8 | Q 0.4546 E 3 = 0.0000 E 8
0.5454 E 8
(ii) 0.9432 E – 4
– 0.0635 E – 4 | Q 0.6353 E – 5 = 0.0635 E – 4
0.8797 E – 4
Multiplication
Two numbers are multiplied in the normalized floating point mode by multiplying
the mantissas and adding the exponents. After the multiplication of the mantissas,
the resulting mantissa is normalized as in an addition or subtraction operation and
the exponent is appropriately adjusted.
Example 1.10 Multiply the following floating point numbers:
(i) 0.5543 E 12 and 0.4111 E – 15
(ii) 0.1111 E 10 and .1234 E 15
(iii) 0.1111 E 51 and 0.4444 E 50
(iv) 0.1234 E – 49 and .1111 E – 54.
Solution: (i) 0.5543 E 12 × 0.4111 E – 15 = 0.2278 E – 3
(ii) 0.1111 E 10 × 0.1234 E 15 = 0.1370 E 24
22 Self-Instructional Material
(iii) 0.1111 E 51 × 0.4444 E 50 = 0.4937 E 100 Errors and Floating Point
Arithmetic
The result overflows.
(iv) 0.1234 E – 49 × 0.1111 E – 54 = 0.1370 E – 104
The result underflows. NOTES
Example 1.11 Apply the procedure for the following multiplications:
(i) (0.5334 × 109) × (0.1132 × 10–25)
(ii) (0.1111 × 1074) × (0.2000 × 1080)
Indicate if the result is overflow or underflow.
Solution: (i) 0.5334 E 9 × 0.1132 E – 25 = 0.6038 E – 17
(ii) 0.1111 E 74 × 0.2000 E 80 = 0.2222 E 153
Hence, the above result overflows.
Example 1.12 Multiply 0.1234 E – 75 by 0.1111 E – 37.
Solution: 0.1234 E – 75 × 0.1111 E – 37
= (0.1234 × 0.1111) E (–75 – 37)
= 0.01370974 E – 112 (before normalization)
= 0.1370 E – 113 (after normalization)
This is the case of underflow as the exponent is 113.
Division
In division, the mantissa of the numerator is divided by that of the denominator. The
denominator exponent is subtracted from the numerator exponent. The quotient mantissa
is normalized to make the most significant digit non-zero and the exponent is
appropriately adjusted. The mantissa of the result is chopped down to four digits.
Example 1.13 Perform the division operation on following operations:
(i) 0.9998 E 1 ÷ 0.1000 E – 99
(ii) 0.9998 E – 5 ÷ 0.1000 E 98
(iii) 0.1000 E 5 ÷ 0.9999 E 3.
Solution:
(i) 0.9998 E 1 ÷ 0.1000 E – 99 = 0.9998 E 101
Hence the result overflows.
(ii) 0.9998 E – 5 ÷ 0.1000 E 98 = 0.9998 E – 104
Hence the result underflows.
(iii) 0.1000 E 5 ÷ 0.9999 E 3 = 0.1000 E 2.
Example 1.14 Evaluate, applying normalized floating point arithmetic, for the
following:
1 – cos x at x = 0.1396 Radian
Self-Instructional Material 23
Errors and Floating Point Assume, cos (0.1396) = 0.9903
Arithmetic
24 Self-Instructional Material
If you substitute x = 0.7320 E 0 in equation x2 + 2x – 2, you obtain: Errors and Floating Point
Arithmetic
(0.7320 E 0)2 + 2(0.7320 E 0) – 2
= (0.7320 E 0 × 0.7320 E 0)
+ (0.2000 E 1 × 0.7320 E 0 – 0.2000 E 1 NOTES
= (0.5358 E 0) + (0.1464 E 1) – (0.2000 E 1)
= (0.0535 E 1 + 0.1464E 1) – (0.2000 E 1)
= (0.1999 E 1 – 0.2000E 1)
= –0.0001 E 1 = 0.1000 E – 2
Ideally, the expression should yield the value zero. Thus, in an algorithm it is
not advisable to write a branching or looping instruction that compares the value
of an expression with zero.
1.5 SUMMARY
Self-Instructional Material 25
Errors and Floating Point • Error is defined as the difference between exact and approximate values.
Arithmetic
• Errors that exist in a problem before it is solved are called inherent errors.
These errors occur due to incorrect measurements or observations that
may be due to the limitations of the measuring instrument such as
NOTES
mathematical tables, calculators or the digital computer. These errors can
be minimized by taking better data or by using high precision computing
aids.
• The errors which occur when some digits from the number are discarded
are known as truncation errors.
• Round-off errors occur during the process of rounding-off a number.
Rounding-off errors are unavoidable errors due to the limitations of
computing.
• If two numbers represented in a normalized floating point notation are to be
added, the exponents of the two numbers must be made equal and the
mantissa should be shifted appropriately.
• Significant digits: These are digits that carry meaning, contributing to their
precision.
• Rounding off: It is the process of dropping the unwanted digits.
• Truncation errors: The are errors that occur when some digits from the
number are discarded.
• Normalization: It is the shifting of a decimal point to the left of the most
significant digit.
26 Self-Instructional Material
6. The two formats in which real numbers can be written are the fixed-point Errors and Floating Point
Arithmetic
(without exponent) format and the floating point (with exponent) format.
7. If two numbers represented in normalized floating point notation are to be
added, the exponents of the two numbers must be made equal and the
NOTES
mantissa shifted appropriately.
8. In division, the mantissa of the numerator is divided by that of the
denominator.
9. In algorithms it is not advisable to write a branching or looping instruction
that compares the value of an expression with zero.
Short-Answer Questions
1. What are the guidelines or rules for rounding off a number to n-significant
digits?
2. Write a short note on truncation errors.
3. Round off the following numbers to four significant digits: 0.00032217,
35.46735 and 18.265101 and compute Ea, Er and Ep in each case.
4. Which standard defines form at for represented floating point format.
5. On what depends the range of a floating point number.
Long-Answer Questions
1. Discuss the different types of errors that occur during numerical
computations.
2. Describe errors in numerical computations.
3. Add the following floating point numbers:
(i) 0.4546 E 5 and 0.5433 E 5
(ii) 0.4546 E 5 and 0.5567 E 7
(iii) 0.4546 E 3 and 0.5567 E 7
(iv) 0.6434 E 3 and 0.4845 E 3
(v) 0.6434 E 99 and 0.4845 E 99.
4. Given a value of x = 3.51 with an error δx = 0.001, estimate the resulting
error in the function y = x3.
5. Find the relative error in the function
m m m
y = ax1 2 x2 2 ... xn n .
6. Explain some drawbacks of floating point representation with the help of an
example.
Self-Instructional Material 27
Errors and Floating Point
Arithmetic 1.9 FURTHER READING
Mott, J.L. Discrete Mathematics for Computer Scientists, 2nd edition. New
NOTES Delhi: Prentice-Hall of India Pvt. Ltd., 2007.
Bonini, C.P., W.H. Hausman H. and Har Bierman. Quantitative Analysis for
Business Decisions. Illinois: Richard D. Irwin, 1986.
Charnes, A., W.W. Cooper, and A. Henderson. An Introduction to Linear
Programming. New York: John Wiley & Sons, 1953.
Gupta, S.C. and V.K. Kapoor. Fundamentals of Mathematical Statistics, 9th
revised edition. New Delhi: Sultan Chand & Sons, 1997.
28 Self-Instructional Material
Interpolation
UNIT 2 INTERPOLATION
Structure NOTES
2.0 Introduction
2.1 Unit Objectives
2.2 Polynomial Interpolation
2.2.1 Finite Differences
2.2.2 Differences of a Polynomial
2.2.3 Some Useful Symbols
2.3 Missing Term Technique
2.3.1 Effect of an Error on a Difference Table
2.4 Newton’s Formulae for Interpolation
2.4.1 Newton–Gregory Forward Interpolation Formula
2.4.2 Newton–Gregory Backward Interpolation Formula
2.5 Central Difference Interpolation Formula
2.5.1 Gauss’s Forward Difference Formula
2.5.2 Gauss’s Backward Difference Formula
2.5.3 Stirling’s Formula
2.5.4 Bessel’s Interpolation Formula
2.5.5 Lagrange’s Interpolation Formula
2.6 Newton’s Divided Difference Interpolation Formula
2.7 Discrete Least Squares Approximation
2.7.1 Fourier Series
2.8 Summary
2.9 Key Terms
2.10 Answers to ‘Check Your Progress’
2.11 Questions and exercises
2.12 Further reading
2.0 INTRODUCTION
The process of estimating or approximating the value of independent variable y for
a given value of x in the specified range say xi ≤ x ≤ xn is known as ‘interpolation’.
If the approximating function φ(x) is a polynomial, the process is known as
‘polynomial interpolation’.
Sometimes, you want to estimate the value of y for any given value of x
which is outside the range xi ≤ x ≤ xn, then the process of approximating y = f(x)
by another function φ(x) is known as extrapolation.
Interpolation is a very important topic in numerical analysis as it provides
the base for numerical differentiation and integration.
The study of interpolation is based on the calculus of finite differences. Finite
differences deal with the variations in a function corresponding to the changes in
the independent variables. These changes occur when independent variable x
takes finite jumps which may be equal or unequal. In this unit, you will study the
Self-Instructional Material 29
Interpolation variations in the function when the independent variable changes by equal intervals.
Since finite differences are the base of interpolation, our discussion will begin with
the definition of finite differences only.
In this unit, you will learn about polynomial interpolation, missing term
NOTES technique, Newton’s formulae for interpolation, Newton’s divided difference in
interpolation formula and discrete least squares approximation.
30 Self-Instructional Material
Interpolation provides the base for numerical differentiation and integration. Interpolation
∆ 2 yn −1 = ∆yn − ∆yn −1
and so on.
In general, the rth forward differences are given as
∆ r yk = ∆ r −1 yk +1 − ∆ r −1 yk .
These differences are systematically set out in Table 2.1.
Self-Instructional Material 31
Interpolation Table 2.1 Forward Difference Table
The higher order differences can be expressed in terms of y, which are here
known as entries, whereas, x is called the argument in the difference table.
The first value y0 in Table 2.1 is called the leading term and the differences
∆y0 , ∆ 2 y0 , ...,etc. are known as the leading differences.
32 Self-Instructional Material
Interpolation
∆ 4 y0 = y4 − 4C1 y3 + 4C2 y2 − 4C3 y1 + 4C4 y0
Thus, in general, you can say that,
∇ 2 y2 = y3 − 2 y2 − y1
and so on.
Self-Instructional Material 33
Interpolation Table 2.2 shows the backward difference table:
Table 2.2 Backward Difference Table
Self-Instructional Material 35
Interpolation The second order divided differences are obtained as:
∆| y2 − ∆| y1 ∆| y3 − ∆| y2
∆| 2 y 1 = ; ∆| 2y2 =
x3 − x1 x4 − x2
NOTES In general,
i ∆| i −1yk +1 − ∆| i −1yk
∆| yk =
xi + k − xk
Table 2.4 represents the divided difference table.
Again, ∆ 2 f ( x) = ∆f ( x + h) − ∆ f ( x)
= b1 + b2 ( x + h) + b3 ( x + h)2
+ ... + bn −1 ( x + h)n − 2 + nh an ( x + h) n −1
36 Self-Instructional Material
Interpolation
– [b1 + b2 x + b3 x 2 + ... + bn −1x n − 2 + nan hx n −1 ]
Thus, E n y r = yr + n
Now, ∆y0 = y1 − y0 ⇒ Ey0 − y0 = ( E − 1) y0
⇒ ∆ =E–1 or E ≡1+ ∆
= E 3 y0 − 3E 2 y0 + 3Ey0 − 1
= y3 − 3 y2 + 3 y1 − 1
(2) Mean operator (µ): Defined by equation,
1
( yr + ½ + yr – ½ )
µyr =
2
Relation between operators: (i) E = 1 + ∆
Proof: ∆yx = y x + h − y x
= Ey x − yn = ( E − 1) y x
Self-Instructional Material 37
Interpolation
⇒ ∆ = E – 1 or E = 1 + ∆
= y x − E −1 y x = (1 − E −1 ) y x ⇒ ∇ ≡ 1 − E −1 .
⇒ δ ≡ E1/2 − E – 1/2
1 1/2
(iv) µ = ( E + E1/2 )
2
1
µyx = ( y x − h/2 + y x − h/2 )
2
1 1/2
⇒ µyx = ( E + E −1/2 ) y x
2
1 1/2
⇒ µ≡ ( E + E −1/2 )
2
(v) ∆ = E∇ = ∇E = δE1/2
Proof: E (∇y x ) = E ( y x − y x − h ) = y x + h − y x = ∆y x
⇒ E∇ ≡ ∆
Similarly,
∇( Ey x ) = ∇( y x + h ) = y x + h − y x = ∆
⇒ ∇E ≡ ∆
δ( E1/2 ( y x )) = δ( y x + h/2 )
= y h h −y h h = yx + h − yx = ∆
x+ + x+ −
2 2 2 2
δE1/2 ≡ ∆
(vi) E ≡ ehD
Ef(x) = f(x + h)
38 Self-Instructional Material
h2
Interpolation
= f ( x) + hf ′( x) + f ′′( x) + ...
2
h2 2
= f ( x) + hD f ( x) + D f ( x) + ... NOTES
2
h2 D 2
= 1 + hD +
2
+ ... f ( x )
= e hD f ( x)
⇒ E ≡ e hD or 1 + ∆ = e hD or ∆ = e hD − 1
Example 2.1: Using the method of separation of symbols, show that:
n(n − 1)
∆ nu x − n = u x − nu x −1 + u x − 2 + ... + ( −1) n u x − n
2
n( n − 1)
Solution: R.H.S. u x − nu x −1 + u x − 2 + ... + ( −1) n u x − n
2
−1 n(n − 1) −2
= 1 − nE + E + ... + (−1) n E − n u
2 x
n n
−1 n 1 E −1
(
= 1− E ) . ux = 1 − ux =
E E
ux
∆n
= n
u x = ∆ n . E − nu x = ∆ nu x − n .
E
Hence Proved. ( Q E – 1 = ∆)
Example 2.2: Show that:
x2 2 2
e x u0 + x ∆u0 + ∆ u0 + ... = u0 + u1 x + u2 x + ...
2
2
x x2∆ 2
e
Solution: L.H.S. 1 + x ∆ + + ... . u0 = e x . e x∆ . u = e x + x∆ . u
2
0 0
= e x (1 + ∆ ) . u0 = e xE u0
x2 E 2
= 1 + xE +
2
+ ... u0
Self-Instructional Material 39
Interpolation
x2 2
= u0 + xEu0 + E u0 + ...
2
x2
NOTES = u0 + x u1 + u2 + ...
2
x+h−x −1 h
= tan −1 = tan
2
1 + x ( x + h) 1 + hx + x
(ii) ∆2 cos 2x = ∆[cos 2(x + h) – cos 2x]
= [cos 2(x + 2h) – cos 2(x + h)] – [cos 2(x + h) – cos 2x]
= – 2 sin (2x + 3h) sin h + 2 sin (2x + h) sin h
= – 2 sin h [2 cos (2x + 2h) sin h] = – 4 sin2 h cos 2(x + h).
Example 2.4: Evaluate:
5 x + 12
∆2 2
x + 5 x + 6
; the interval of differenciation being unity..
5 x + 12
Solution: ∆ 2
( x + 2) ( x + 3)
2 3 2 3
= ∆2 + = ∆ ∆ + ∆
x + 2 x + 3 x + 2 x + 3
1 1 1 1
= ∆ 2 − + 3 −
x+3 x + 2 x + 4 x + 3
1 1
= −2∆ − 3∆
( x + 2) ( x + 3) ( x + 3) ( x + 4)
1 1
= −2 −
( x + 3) ( x + 4) ( x + 2) ( x + 3)
1 1
−3 −
( x + 4) ( x + 5) ( x + 3) ( x + 4)
40 Self-Instructional Material
Interpolation
4 6
= +
( x + 2) ( x + 3) ( x + 4) ( x + 3) ( x + 4) ( x + 5)
2(5 x + 16)
= . NOTES
( x + 2) ( x + 3) ( x + 4) ( x + 5)
Example 2.5: Prove that hD = –log (1 – ∇) = sin h–1 (µδ).
Solution: hD = log E = – log (E–1) = – log (1 – ∇) [ Q E–1 = 1 – ∇]
1 1/2
Also, µ= ( E + E −1/2 )
2
δ = E1/2 − E −1/2
1 −1 1 hD − hD
∴ µδ = ( E − E ) = (e − e ) = sin h ( hD )
2 2
or hD = sin h −1 (µ δ) .
Example 2.6: (i) Find the function whose first difference is ex.
(ii) If ∆ 3u x = 0, prove that,
u 1 1
x+
1
= (u x + u x + 1 ) − ( ∆ 2u x + ∆ 2u x + 1 ).
2 2 16
∆e x
⇒ e =x
eh − 1
ex x ex
Hence, ∆ h = e or f ( x) = .
e −1 eh − 1
u 1
(ii) x+ = E1/2u x = (1 + ∆ )1/2 u x
2
1 1 2
= 1 + ∆ − ∆ u x [ Q ∆3ux = 0]...(1)
2 8
Now, ∆ 3u x = 0
⇒ ∆3u x + 1 − ∆ 2u x = 0
⇒ ∆ 2u x + 1 = ∆ 2u x and ∆u x = u x + 1 − u x
Self-Instructional Material 41
Interpolation ∴ From Equation (1),
1 ∆ 2u x ∆ u x + 1
2
1
u = u x + (u x + 1 − u x ) − +
1 2 8 2 2
NOTES
x+
2
1 1
= (u x + u x + 1 ) − ∆ 2u x + ∆ 2 u x + 1 .
( )
2 16
Example 2.7: (i) Find f(6) if f(0) = –3, f(1) = 6, f(2) = 8, f(3) = 12; third
difference being constant.
(ii) Find ∆10 (1 − ax)(1 − bx 2 )(1 − cx3 )(1 − dx 4 ).
E1/2 + E – 1/2
= f ( x).( E1/2 − E −1/2 ) g ( x)
2
E1/2 + E – 1/2 NOTES
+ g ( x) ( E1/2 − E −1/2 ) f ( x)
2
1 1 1 1 1
= f x + + f x – g x + − g x −
2 2 2 2 2
1 1 1 1
+ g x + + g x − f x + − f x −
2 2 2 2
1 1 1 1 1 1 1
= f x + g x + − f x + g x − + f x– gx+
2 2 2 2 2 2 2
1 1 1 1
f x– gx– + f x+ gx+
2 2 2 2
1 1 1 1 1 1
+ f x + g x − − f x − g x + − f x − g x −
2 2 2 2 2 2
1 1 1 1
= f x+ gx+ − f x− gx −
2 2 2 2
= E1/2 f ( x) g ( x) − E −1/2 f ( x) g ( x)
= ( E1/2 − E −1/2 ) f ( x) g ( x) = δf ( x) g ( x).
µ g ( x) δf ( x) − µ f ( x) δg ( x)
(ii) R.H.S. =
1 1
gx − gx +
2 2
Numerator of R.H.S.
E1/2 + E −1/2
= g ( x) ( E1/2 − E −1/2 ) f ( x)
2
E1/2 + E −1/2
– f ( x) ( E1/2 − E −1/2 ) g ( x)
2
1 1 1 1 1
= g x + + g x − f x+ − f x −
2 2 2 2 2
1 1 1 1
− f x+ + f x − g x + – g x –
2 2 2 2
Self-Instructional Material 43
Interpolation
1 1 1 1 1 1 1
= f x+ gx+ + f x+ gx− − f x− gx+
2 2 2 2 2 2 2
NOTES 1 1 1 1 1 1 1
– f x – g x – − f x+ gx+ − f x+ gx–
2 2 2 2 2 2 2
1 1 1 1
+ f x− gx+ − f x − g x −
2 2 2 2
1 1 1 1
= f x+ gx– − f x− gx+
2 2 2 2
1 1 1 1
f x+ gx− − f x− gx+
R.H.S. =
2 2 2 2
∴
1 1
gx− gx+
2 2
1 1
f x+ f x−
=
2
−
2
1 1
gx+ gx−
2 2
f ( x) f ( x)
= E1/2 − E1/2 g ( x)
g ( x)
1/2 – 1/2 f ( x ) f ( x)
= (E − E ) =δ .
g ( x) g ( x)
Suppose n values out of (n + 1) values of y = f(x) are given, the values of x being
equidistant.
Let the unknown value be N.
44 Self-Instructional Material
Since, only n values of y are known, you can assume y = f(x) to be a Interpolation
polynomial of degree (n – 1) in x.
By equating to zero the nth difference, you can get the value of N.
Example 2.9: Express f(x) = x4 – 12x3 + 24x2 – 30x + 9 and its successive NOTES
difference in factorial notation. Hence, show that ∆5f(x) = 0.
∆ 4 f ( x) = 24
and ∆5 f ( x) = 0.
Example 2.10: Using the method of separation of symbols, show that
1 1 1
u0 − u1 + u2 − u3 + ... = u0 − ∆u0 + ∆ 2u0 − ... .
2 4 8
1 1
2 3
1 1
Solution: R.H.S. = 1 − ∆ + ∆ − ∆ + ... u0
2 2 2 2
−1
1 1 1 1
= . u0 = 1 + ∆ u0
2 1 2 2
1 + ∆
2
Self-Instructional Material 45
Interpolation
= (2 + ∆ )−1 u0 = (1 + E )−1 u0
= (1 – E + E2 – E3 + ...) u0
= u0 − u1 + u2 − u3 + ... = L.H.S.
NOTES
Example 2.11: Use the method of separation of symbols to prove the following
identities:
x x
2 E + ∆2
−1 x E2 − E +1
= (1 + ∆ E ) u x = u x = u x
E E
1
= x
[1 + E ( E − 1)]x u x
E
= E − x (1 + ∆E ) x u x = (1 + ∆E ) x u0
= (1 + x C1 ∆E + x C2 ∆ 2 E 2 + ...) u0
x x 2
= u + C1 ∆u1 + C2 ∆ u2 + ... = R.H.S.
−1 x +1
(ii)
x
R.H.S. = un + C1 ∆E un + C2 ∆ 2 E −2u2
x+2
+ C3 ∆ 3 E −3un + ...
−1 x +1
x
= (1 + C1 ∆E + C2 ∆ 2 E −2 + ....) un
= (1 − ∆E −1 )− x un
−x −x
∆ E−∆
= 1 − un = un
E E
−x
1
= un = E x un = un + x = L.H.S.
E
= (1 + E + E 2 + ... + E n ) u0
46 Self-Instructional Material
E n +1 − 1 (1 + ∆ ) n +1 − 1
Interpolation
= E − 1 0 u = u0
∆
1
= (1 + n +1 n +1
C2 ∆ 2 + n +1
C3 ∆ 3 + ... + ∆ n + 1 ) – 1 u0 NOTES
∆
C1∆ +
n +1 n +1 n +1
= C1 u0 + C 2 ∆u 0 + C3 ∆ 2u0 + ... ∆ nu0 = R.H.S.
2.3.1 Effect of an Error on a Difference Table
Suppose there is an error ε in the entry of yS of a table. As higher differences
are formed this error spreads out and is considerably magnified, as shown in
Table 2.5.
Table 2.5 Effect of Error on a Difference Table
x y ∆y ∆ 2y ∆ 3y ∆ 4y
x0 y0
∆y 0
x1 y1 ∆2 y 0
∆y 1 ∆3 y 0
2
x2 y2 ∆ y1 ∆4 y 0
3
∆y 2 ∆ y1
x3 y3 ∆2 y 2 ∆4y1 + ε
3
∆y 3 ∆ y2 – ε
x4 y4 ∆2y3 + ε ∆4y2 – 4ε
∆y4 + ε ∆3y3 – ε
x5 y5 + ε ∆2 y4 –2ε ∆4y3 + 6ε
3
∆y5 – ε ∆ y4 + 3ε
2
x6 y6 ∆ y5 + ε ∆4y4 – 4ε
3
∆y 6 ∆ y5 – ε
x7 y7 ∆2 y 6 ∆4y5 + ε
3
∆y 7 ∆ y6
x8 y8 ∆2 y 7
∆y 8
x9 y9
2 2 2 2
(1 – ε)2 = 1 − C1 ε + C1 ε = 1 − 2ε + ε
4 4 2 4 3 4 4
Similarly, (1 – ε)4 = 1 − C1 ε + C2 ε – C3 ε + C4 ε
2 3 4
= 1 − 4ε + 6ε − 4ε + ε .
Self-Instructional Material 47
Interpolation (iii) The algebraic sum of all the errors in any column is zero.
(iv) The maximum error in each column occurs opposite to the entry
containing the error, i.e., here opposite to y5.
NOTES Example. 2.12: Find the error and correct the wrong figure in the following
functional values.
2, 5, 10, 18, 26, 37, 50
Solution: Following table shows the functional values which are as follows:
Table 2.6 Table for the Functional Values 2, 5, 10, 18, 26, 36, 50
x y ∆y ∆2 y ∆3 y ∆4 y ∆5 y ∆6 y
1 2
3
2 5 2
5 1
3 10 3 –4
8 –3 10
4 18 0 6 –20
8 3 –10
5 26 3 –4
11 –1
6 37 2
13
7 50
Sum of all the third differences is zero.
Let ε be the error. Now the two adjacent values are equal in magnitude,
i.e., – 3 and 3. The coefficients in this column = Binomial coefficient in the expression
of (1 – ε)3 = 1 – 3ε + 3ε2 – ε3
⇒ –3ε = –3 ⇒ ε =1
⇒ Error lies in the value of y4, i.e., 18.
So corrected value is 18 – 1 = 17
Example. 2.13: If f(x) = y is a polynomial of degree 3 and following gives the
values of x and y then locate and correct the error.
Values of x and y
x: 0 1 2 3 4 5 6
y: 4 10 30 75 160 294 490
48 Self-Instructional Material
Solution: Following table represents the correction of errors: Interpolation
Table 2.12 Correction of Errors in the Table
x y ∆y ∆2 y ∆3 y
0 4 NOTES
6
1 10 14
20 11+ (1)
2 30 25
45 15 – 3(1)
3 75 40 Sum = 11 + 15 + 9 + 13
+ε 85 9 + 3(1) = 20 + 15 + 13
4 160 49 = 48/4 = 12
134 13 – (1)
5 294 62
196
6 490
Since, the polynomial is of degree 3, ∆3y should be a constant and the error
corresponds to the coefficient.
So, the entry corresponding to x = 3 is incorrect, it should be 75 – 1 = 74.
Self-Instructional Material 49
Interpolation
For x = a + h, f(a + h) = A0 + A1h
= f (a ) + 2h
∆f ( a ) 2
+ 2h A2
h
⇒ 2h 2 A2 = f (a + 2h) − 2 f (a + h) + f (a) = ∆ 2 f (a)
∆ 2 f (a)
⇒ A2 =
2! h2
∆3 f (a )
Similarly, A3 = and so on.
3! h3
∆ n f (a)
Thus, An = .
n! hn
From Equation (2.2),
∆ f (a) ∆ 2 f (a)
f(x) = f (a) + ( x − a) h + ( x − a) ( x − a − h) + ...
2! h 2
∆ n f (a)
+( x − a) ... ( x − a − n − 1 h)
n! hn
x−a
Put, x = a + hk ⇒ k= , you will have
h
∆ f (a) (hk ) (hk − h) 2
f(a + hk) = f (a) + hk + ∆ f (a) + ...
h 2! h 2
( hk ) ( hk − h) ( hk − 2h) ... ( hk − n − 1 h)
+ n
∆ n f (a)
n! h
k (k − 1) 2
⇒ f (a + hk ) = f (a) + k ∆ f (a) + ∆ f (a) + ...
2!
k (k − 1) (k − 2) ... (k − n + 1) n
+ ∆ f (a)
n!
this is the required formula.
50 Self-Instructional Material
This formula is particularly useful for interpolating the values of f(x) near the Interpolation
beginning of the set of values given. h is called the interval of difference, while ∆ is
the forward difference operator.
Example.2.14: From the following table, estimate the number of students who
NOTES
weight between 45 and 50.
x yx ∆y ∆2 y ∆3 y ∆4 y
45 20
45
55 65 –10
35 –13
65 100 –23 34
12 21
75 112 –2
10
85 122
By taking x0 = 45, you shall find the number of students with weight less
than 50.
50 − 45
So, k = = 0.5
10
Using Newton’s Forward Interpolation formula, you get
k ( k − 1) 2
y50 = y45 + k ∆y45 + ∆ y45 + ...
2!
0.5 (0.5 − 1)
∴ y50 = 20 + 0.5 × 45 + × ( −10)
2!
Self-Instructional Material 51
Interpolation But the number of students with weight less then 45 is 20.
Hence, the number of students with weight between 45 kg and 50 kg are
42 – 20 = 22.
NOTES Example 2.15: Using Newton’s Forward interpolation formula obtain a polynomial
in x which takes the following values.
x: 4 6 8 10
y: 1 3 8 16
x−4
Here, x0 = 4, K=
2
x−4 x−4
− 1
x−4 2 2 × (3) + 0
So, y(x) = 1 + × 2 +
2 2!
( x − 4) ( x − 6)
= 1 + (x – 4) + ×3
2
( x − 6).3
= 1 + (x – 4) 1 +
2
x−4 1
=1+ [3 x − 16] = 1 + (3 x 2 − 28 x + 64)
2 2
3 2
y(x) = x − 14 x + 33
2
52 Self-Instructional Material
Interpolation
3
So, y(5) = (5) 2 − 14(5) + 33 = 0.5
2
2.4.2 Newton–Gregory Backward Interpolation Formula
NOTES
Let y = f(x) be a function of x which assumes the value f(a), f(a + h), f(a + 2h),
..., f(a + nh) for (n + 1) equidistant values a, a + h, a + 2h, ..., a + nh of the
independent variable x.
Let f(x) be a polynomial of the nth degree.
+ An ( x − a − nh) ( x − a − n − 1 h) + ... ( x − a − h)
∇f ( a + nh)
⇒ A1 = (2.7)
h
Put, x = a + (n – 2)h, then,
= ∇ 2 f (a + nh)
∇ 2 f (a + nh)
A2 = (2.8)
2! h 2
Proceeding, you get,
∇ n f (a + nh)
An = (2.9)
n ! hn
Substituting the values in Equation (2.6), you get,
∇f ( a + nh)
f(x) = f ( a + nh) + ( x − a − nh) + ...
h
∇ n f (a + nh)
+ ( x − a − nh) ( x − a − n − 1 h) ... (2.10)
n! hn Self-Instructional Material 53
Interpolation Put x = a + nh + kh, then,
x – a – nh = kh
and x – a – (n – 1)h = (k + 1)h
NOTES M
x – a – h = ( k + n − 1) h
∴ Equation (2.10) becomes,
k ( k + 1) 2
f(x) = f ( a + nh) + k ∇f ( a + nh) + ∇ f ( a + nh)
2!
∇ n f (a + nh)
+ ... + kh .(k + 1)h ... (k + n − 1) (h)
n ! hn
Which is the required formula.
Or k (k + 1) 2
f ( xn + kh) = f ( xn ) + k ∇f ( xn ) + ∇ f ( xn )
2!
k (k + 1)(k + 2) ... (k + n − 1) n
+ ... + ∇ f ( xn )
n!
This formula is useful when the value of f(x) is required near the end of the
table.
Where xn = (x0 + nh) and a = x0, so f(a + nh) = f(xn).
Example 2.16: Using Newton’s backward interpolation formula, obtain the value
of tan 22°, given that:
θ° 0 4 8 12 16 20 24
tan θ : 0 0.0699 0.1405 0.2126 0.2867 0.3640 0.4452
x 104f(x) = y ∆y ∆2 y ∆3 y ∆4 y ∆5 y
0 0
699
7
4 699 8
706 –3
8 1405 15 10
721 5
12 2126 20 7
741 12 –12
16 2867 32 –5
773 7
20 3640 39
812
24 4452
54 Self-Instructional Material
Here, x n = 24, x = 22 Interpolation
x − xn
and k = = – 0.5.
h
Thus, by Newton’s backward interpolation formula, you have NOTES
k ( k + 1) 2 k ( k + 1) ( k + 12) 3
y(x) = yn + k ∇yn + ∇ yn + ∇ yn + ...
2! 3!
y(22) = 4452 + (–0.5) × 812
( −0.5).(0.5) ( −0.5).(0.5).(1.5)
+ × 39 + × (7)
2! 3!
You shall now study the central difference formulae most suited for Interpolation
near the middle of a tabulated set.
2.5.1 Gauss’s Forward Difference Formula
Newton–Gregory forward difference formula is,
k ( k − 1) 2 k ( k − 1) ( k − 2) 2
f(a + hk) = f(a) + k∆f(a) + ∆ f (a) + ∆ f (a)
2! 3!
k ( k − 1) ( k − 2) ( k − 3) 4
+ ∆ f ( a ) + ... (2.11)
4!
Given a = 0, h = 1, you get
k ( k − 1) 2 k ( k − 1) ( k − 2) 3
f(k) = f(0) + k∆f(0) + ∆ f (0) + ∆ f (0)
2! 3!
k ( k − 1) ( k − 2) ( k − 3) 4
+ ∆ f (0) + ... (2.12)
4!
Now,
k ( k − 1) (k − 2) 3
+ {∆ f ( −1) + ∆ 4 f ( −1)}
3!
k ( k − 1) (k − 2) ( k − 3) 4
+ {∆ f ( −1) + ∆5 f ( −1)} + ...
4!
k ( k − 1) 2 k ( k − 1) k − 2 3
= f (0) + k ∆f (0) + ∆ f ( −1) + 1 + ∆ f ( −1)
2! 2 3
k ( k –1) ( k − 2) k − 3 4 k ( k − 1) ( k − 2) ( k − 3) 5
+ 1 + ∆ f ( −1) + ∆ f (−1) + ...
6 4 4!
k (k − 1) 2 ( k + 1) k ( k − 1) 3
= f (0) + k ∆ (0) + ∆ f ( −1) + ∆ f (−1)
2! 3!
(k + 1) k (k − 1) ( k − 2) 4 k (k − 1) (k − 2) (k − 3) 5
+ ∆ f (−1) + ∆ f (−1) + ...
4! 4!
(2.13)
56 Self-Instructional Material
Interpolation
k (k − 1) 2 (k + 1) k (k − 1) 3
f (u ) = f (0) + k ∆f (0) + ∆ f (−1) + ∆ f (−1)
2! 3!
k (k + 1) (k − 1) (k − 2) 4
+ ∆ f (−2) + ...
4! NOTES
x: 30 35 40 45 50
f ( x) : 3678.2 2995.1 2400.1 1876.2 1416.3
( k + 1) k ( k − 1) ( k − 2) 4
+ ∆ y−2 + ...
4!
The central difference table is shown below:
Central Difference Table
x k f(x) ∆y ∆2 y ∆3 y ∆4 y
30 –2 3678.2
–683.1
35 –1 2995.1 88.1
–595 –17
40 0 2400.1 71.1 9.9
–523.9 –7.1
45 1 1876.2 64
–459.9
50 2 1416.3
Thus,
41 − 40 1
Value of k = = = 0.2
45 − 40 5
0.2 (0.2 − 1)
Now, f(41) = 2400.1 + (0.2) × (–523.9) + × (71.1)
2!
Self-Instructional Material 57
Interpolation
(0.2 + 1) (0.2) (0.2 − 1) (0.2 + 1) (0.2) (0.2 − 1) (0.2 − 2)
+ × ( −7.1) + × (9.9)
3! 4!
= 2400.1 – 104.78 – 5.688 + 0.2272 + 0.14256
NOTES = 2290.001760.
2.5.2 Gauss’s Backward Difference Formula
Newton–Gregory forward difference formula is,
k ( k − 1) 2
f(a + hk) = f ( a ) + k ∆f ( a ) + ∆ f (a)
2!
k ( k − 1) ( k − 2) 3
+ ∆ f ( a ) + ... (2.14)
3!
Putting a = 0, h = 1, you get
k ( k − 1) 2 k ( k − 1) ( k − 2) 3
f(k) = f (0) + k ∆f (0) + ∆ f (0) + ∆ f (0)
2! 3!
k ( k − 1) ( k − 2) ( k − 3) 4
+ ∆ f (0) + ... (2.15)
4!
Now, ∆f (0) = ∆f (−1) + ∆ 2 f (−1)
k −1 2
= f (0) + k ∆f ( −1) + k 1 + ∆ f ( −1)
2
k ( k − 1) k −2 3
+ 1 + ∆ f ( −1)
2 3
58 Self-Instructional Material
Interpolation
k ( k − 1) ( k − 2) k − 3 4
+ 1 + ∆ f ( −1)
6 4
k ( k − 1) ( k − 2) ( k − 3) 5
+ ∆ f ( −1) + ... NOTES
4!
(k + 1) k 2
= f (0) + k ∆f ( −1) + ∆ f ( −1)
2!
k ( k + 1) ( k − 1) 3
+ ∆ f ( −1)
3!
( k + 1) k ( k − 1) ( k − 2) 4
+ ∆ f ( −1) + ... (2.17)
4!
( k + 1) k ( k − 1)
+ {∆3 f (−2) + ∆ 4 f (−2)}
3!
( k + 1) k ( k − 1) ( k − 2) 4
+ {∆ f ( −2) + ∆ 5 f ( −2)} + ...
4!
Thus, Gauss’ backward formula is,
k ( k + 1) 2 ( k + 1) k ( k − 1) 3
y = y0 + k ∆y−1 + ∆ y −1 + ∆ y −2
2! 3!
( k + 2) ( k + 1) k ( k − 1) 4
+ ∆ y−2
4!
( k + 2) ( k + 1) k ( k − 1) ( k − 2) 5
+ ∆ y−2 + ...
5!
Example 2.18: Apply Gauss’s backward formula to compute sin 45° from the
following table:
θ0 : 20 30 40 50 60 70 80
sin θ : 0.34202 0.502 0.6479 0.76604 0.86603 0.93969 0.98481
Self-Instructional Material 59
Interpolation Solution: The difference table is shown below:
Difference Table to Compute sin 45°
Here, x 0 = 40, x = 45
45 − 40
and k = = 0.5.
10
Thus, by Gauss’ backward formula, you will have
k ( k + 1) 2 ( k + 1) k ( k − 1) 3
y(x) = y0 + k ∆y−1 + ∆ y −1 + ∆ y −2
2! 3!
( k + 2) ( k + 1) k ( k − 1) 4 ( k + 2) ( k + 1) k ( k − 1) ( k − 2) 5
+ ∆ y −2 + ∆ y−3 + ...
4! 5!
you have,
1.5 × 0.5
y(45) = 0.64279 + 0.5 × 0.14079 + × (–0.01754)
2!
1.5 × 0.5 ( −0.5) (2.5) (1.5) (0.5) ( −0.5)
+ × (0.00165) + × ( −0.0074)
3! 4!
= 0.64279 + 0.070395 + (–0.0065775) – (0.000103125)
– 0.00028789
= 0.70679.
2.5.3 Stirling’s Formula
Gauss’ forward and backward formula are used to derive Stirling formula.
Gauss’ forward formula is,
k ( k − 1) 2 ( k + 1) k ( k − 1) 3
f(k) = f (0) + k ∆f (0) + ∆ f ( −1) + ∆ f ( −1)
2! 3!
( k + 1) k ( k − 1) ( k − 2) 4
+ ∆ f ( −2) + ... (2.18)
60 Self-Instructional Material
4!
Gauss’ backward formula is, Interpolation
( k + 1) 2 ( k + 1) k ( k − 1) 3
f(k) = f (0) + k ∆ ( f − 1) + ∆ f ( −1) + ∆ f ( −2)
2! 3!
( k + 2) k ( k + 1) ( k − 1) 4 NOTES
+ ∆ f ( −2) + ... (2.19)
4!
Take the mean of Equations (2.18) and (2.19),
∆f (0) + ∆f ( −1) k2 2
f ( x ) = f (0) + + ∆ f ( −1)
2 2!
(k + 1) k (k − 1) ∆ 3 f ( −1) + ∆ 3 f ( −2)
+
3! 2
k 2 ( k 2 − 1) 4
+ ∆ f ( −2) + ... (2.20)
4!
1 1 1
This is called Stirling’s formula. It is useful when | k | < or − < k < . It
2 2 2
1 1
gives the best estimate when – <k< .
4 4
Example 2.19: Use Stirling’s formula to evaluate f(1.22) given,
x − 1.2
Taking x0 = 1.2 and h = 0.2 p =
0.2
Self-Instructional Material 61
Interpolation Using Stirling’s formula,
2
0.031 + 0.041 (0.2)
y(1.22) = 0.932 + (0.2) × + × (−0.01)
2 2!
NOTES
(0.2) [(0.2)2 − 12 ] 0.001 − 0.001 (0.2) 2 [(0.2)2 − 12 ]
+ × + × (0.002)
3! 2 4!
= 0.932 + 0.0072 – 0.0002 – 0.0000032
= 0.9389 (approx.)
2.5.4 Bessel’s Interpolation Formula
Gauss’ forward formula is,
k ( k − 1) 2 ( k + 1) k ( k − 1) 3
f(k) = f (0) + k ∆ (0) + ∆ f ( −1) + ∆ f ( −1)
2! 3!
( k + 1) k ( k − 1) ( k − 2) 4
+ ∆ f ( −2) ... (2.21)
4!
Gauss’s backward formula is,
k (k + 1) 2 ( k + 1) k ( k − 1) 3
f(k) = f (0) + k ∆f (−1) + ∆ f (−1) + ∆ f (−2)
2! 3!
( k + 2) k ( k + 1) ( k − 1) 4
+ ∆ f ( −2) + ... (2.22)
4!
In Equation (2.22), shift the origin to 1 by replacing k by k – 1 and adding
1 to each argument 0, –1, –2, ..., you get
k ( k − 1) 2
f(k) = f (1) + ( k − 1) ∆f (0) + ∆ f (0)
2!
k ( k − 1) ( k − 2) 3
+ ∆ f ( −1)
3!
( k + 1) k ( k − 1) ( k − 2) 4
+ ∆ f ( −1) + ... (2.23)
4!
By taking mean of Equations (2.21) and (2.23), you get
f (0) + f (1) k + ( k − 1)
f(k) = + ∆f (0)
2 2
k ( k − 1) ∆ 2 f ( −1) + ∆ 2 f (0)
+
2! 2
62 Self-Instructional Material
Interpolation
k ( k − 1) ∆ 3 f ( −1)
+ ( k + 1 + k − 2)
3! 2
f (0) + f (1) 1
f (k ) = + k − ∆f (0)
2 2
k (k − 1) ∆ 2 f (−1) + ∆ 2 f (0)
+
2! 2
1
(k − 1) k − u
+ 2
∆3 f (−1) (2.24)
3!
(k + 1) k (k − 1) (k − 2) ∆ 4 f (−2) + ∆ 4 f (−1)
+ + ...
4! 2
Example 2.20: Using Bessel’s formula obtain y26. Given that y20 = 2854, y24 =
3162, y28 = 3544 and y32 = 3992.
x − 24
Solution: With x0 = 24 and k = the central difference table is shown
4
below:
Central Difference Table for Obtaining y26
x k y ∆y ∆2 y ∆3 y
20 –1 2854
308
24 0 3162 74
382 –8
28 1 3544 66
448
32 2 3992
Using Bessel’s formula,
0.5 (0.5 − 1) 74 + 66
y26 = 3162 + 0.5 × 382 + × + 0 × (–8)
2! 2
= 3162 + 191 + (–8.75)
= 3344.25.
Self-Instructional Material 63
Interpolation 2.5.5 Lagrange’s Interpolation Formula
Let f ( x0 ), f ( x1 ) ,..., f ( xn ) be (n + 1) entries of a function y = f(x), where
f(x) is assumed to be a polynomial corresponding to the arguments x0, x1, x2,
NOTES ..., xn .
The polynomial f(x) may be written as,
f(x) = A0 ( x − x1 ) ( x − x2 ) ... ( x − xn )
+ A1 ( x − x0 ) ( x − x2 ) ... ( x − xn )
f(x1) = A1 ( x1 − x0 ) ( x1 − x2 ) ... ( x1 − xn )
f ( x1 )
∴ A1 = (2.27)
( x1 − x0 ) ( x1 − x2 ) ... ( x1 − xn )
M M M
f ( xn )
Similarly, An = (2.28)
( xn − x0 ) ( xn − x1 ) ... ( xn − xn −1 )
Substituting the values of A0, A1, ..., An in Equation (2.25), you will get,
( x − x1 ) ( x − x2 ) ... ( x − xn )
f ( x) = f ( x0 )
( x0 − x1 ) ( x0 − x2 ) ... ( x0 − xn )
( x − x0 ) ( x − x2 ) ... ( x − xn )
+ f ( x1 ) ..(2.29)
( x1 − x0 ) ( x1 − x2 ) ... ( x1 − xn )
( x − x0 ) ( x − x1 ) ... ( x − xn −1 )
+ ... + f ( xn )
( xn − x0 ) ( xn − x1 ) ... ( xn − xn −1 )
64 Self-Instructional Material
Interpolation
f ( x) f ( x0 ) 1
= .
( x − x0 ) ( x − x1 ) ... ( x − xn ) ( x0 − x1 ) ( x0 − x2 ) ... ( x0 − xn ) ( x − x0 )
f ( x1 ) 1
+ . + ... NOTES
( x1 − x0 ) ( x1 − x2 ) ... ( x1 − xn ) ( x − x1 )
f ( xn ) 1
+ . . ... (2.30)
( xn − x0 ) ( xn − x1 ) ... ( xn − xn −1 ) ( x − xn )
n
Where, φ(x) = ∏ ( x − xr )
r =0
d
and φ′( xr ) = [φ( x)]
dx x = xr
Proof:
You have the Lagrange’s formula,
n ( x − x0 ) ( x − x1) ... ( x − xr −1) ( x − xr +1) ... ( x − xn )
Pn(x) = ∑ (x − x0 ) ( xr − x1) ... ( xr − xr −1) ( xr − xr +1) ... ( xr − xn )
f ( xr )
r =0 r
n φ( x) f ( xr )
= ∑
r = 0 x − xr ( xr − x0 ) ( xr − x1 ) ... ( xr − xr −1) ( xr − xr +1) ... ( xr − xn )
...(2.31)
Now,
n
φ( x ) = ∏ ( x − xr )
r =0
= ( x − x0 ) ( x − x1 ) ... ( x − xr − 1 ) ( x − xr ) ( x − xr + 1 ) ... ( x − xn )
∴ φ′( x ) = ( x − x1 ) ( x − x2 ) ... ( x − xr ) ... ( x − xn )
+ ( x − x0 ) ( x − x2 ) ... ( x − xr ) ... ( x − xn ) + ...
+ ( x − x0 ) ( x − x1 ) ... ( x − xr −1 ) ... ( x − xr +1 ) ... ( x − xn ) + ...
+ ( x − x0 ) ( x − x1 ) ... ( x − xr ) ... ( x − xn +1 )
Self-Instructional Material 65
Interpolation
⇒ φ′( xr ) = [ φ′( x)]x = xr
Hence proved.
Example 2.21: Find the unique polynomial P(x) of degree 2 such that,
P(1) = 1, P(3) = 27, P(4) = 64
Use the Lagrange’s method of Interpolation.
Solution: Here, x0= 1, x 1 = 3, x2 = 4
f(x0) = 1, f(x1) = 27, f(x2) = 64
Lagrange’s Interpolation formula is,
( x − x1 ) ( x − x2 ) ( x − x0 ) ( x − x2 )
P(x) = f ( x0 ) + k f ( x1 )
( x0 − x1 ) ( x0 − x2 ) ( x1 − x0 ) ( x1 − x2 )
( x − x0 ) ( x − x1 )
+ f ( x2 )
( x2 − x0 ) ( x2 − x1 )
( x − 3) ( x − 4) ( x − 1) ( x − 4) ( x − 1) ( x − 3)
= (1) + (27) + (64)
(1 − 3) (1 − 4) (3 − 1) (3 − 4) (4 − 1) (4 − 3)
1 2 27 64
= ( x − 7 x + 12) − ( x 2 − 5 x + 4) + ( x 2 − 4 x + 3)
6 2 3
= 8x2 – 19x + 12
Hence, the required unique polynomial is,
P(x) = 8x2 – 19x + 12
66 Self-Instructional Material
Interpolation
2.6 NEWTON’S DIVIDED DIFFERENCE
INTERPOLATION FORMULA
Let y0, y1, ..., yn be the values of y = f(x) corresponding to the arguments x0, x1, NOTES
..., xn, then from the definition of divided differences, you have,
y − y0
[x, x0] =
x − x0
so that, y = y0 + ( x – x0 ) [ x, x0 ] ...(2.33)
[ x, x0 ] − [ x0 , x1 ]
Again, [x, x0, x1] =
x − x1
Which gives, [x, x0] = [x0, x1] + (x – x1) [x, x0, x1] (2.34)
From Equations (2.33) and (2.34)
y = y0 + (x – x0) [x0, x1] + (x – x0) (x – x1) [x, x01, x1] ...(2.35)
[ x, x0 , x1 ] − [ x0 , x1 , x2 ]
Also [x, x0, x1, x2] =
x − x2
Which gives, [x, x0, x1] = [x0, x1, x2] + (x – x2) [x, x0, x1, x2] (2.36)
From Equations (2.35) and (2.36),
y = y0 + ( x − x0 ) [ x0 , x1 ] + ( x − x0 ) ( x − x1 ) [ x0 , x1 , x2 ]
+ ( x − x0 ) ( x − x1 ) ( x − x2 ) [ x, x0 , x1 , x2 ]
Proceeding in this manner, you get,
y = f(x) = y0 + ( x − x0 ) [ x0 , x1 ] + ( x − x0 ) ( x − x1 ) [ x0 , x1 , x2 ]
+ ( x − x0 ) ( x − x1 ) ( x − x2 ) [ x0 , x1 , x2 , x3 ]
+ ... + ( x − x0 ) ( x − x1 ) ( x − x2 )
... ( x − xn −1 ) [ x0 , x1 , x2 , x3 , ..., xn ]
+ ( x − x0 ) ( x − x1 ) ( x − x2 )
... ( x − xn ) [ x, x0 , x1 , x2 , ..., xn ]
This is called Newton’s General Interpolation formula with divided differences.
The last term being the remainder term after (n + 1) terms.
Self-Instructional Material 67
Interpolation Example 2.22: Referring to the following table, find the value of f(x) at point
x = 4:
x: 1.5 3 6
NOTES f(x): –0.25 2 20.
Solution: The divided difference table is shown below:
Divided Difference Table for Finding the Value of f(x)
1.5 − 0.25
1.5
3 2 1
6
6 20
( x + 1) x ( x − 1) 3
+ ∆ f ( −2) + ...
3!
Solution: Taking the arguments, 0, – 1, 1, –2, ... the divided Newton’s difference
formula is,
2
f(x) = f (0) + x ∆| f (0) + x ( x + 1) ∆| f (0)
−1 −1, 1
2
= f (0) + x ∆| f ( −1) + x( x + 1) ∆| f ( −1)
0 0, 1
+ ( x + 1) x ( x – 1) ∆| 3 f ( −2) + ...
−1, 0, 1
68 Self-Instructional Material
Interpolation
f (0) − f (−1)
Now, ∆| f (−1) = = ∆f (−1)
0 0 − (−1)
1
∆| 2 f ( −1) = [ ∆| f (0) − ∆| f (−1)] NOTES
0, 1 1 − (−1) 1 01
1 1
= [ ∆f (0) − ∆f ( −1)] = ∆ 2 f ( −1)
2 2
1
∆| 3 f (−2) = [ ∆| 2 f (−1) − ∆| 2 f (−2)]
−1, 0, 1 1 − (−2) 0, 1 −1, 0
1 ∆ 2 f (−1) ∆ 2 f (−2)
= −
3 2 2
∆ 3 f (−2) ∆ 3 f (−2)
= = and so on.
3.2 3
Substituting these values, you get,
( x + 1) x 2
f(x) = f (0) + x∆f ( −1) + ∆ f ( −1)
2!
( x + 1) x ( x − 1) 3
+ ∆ f ( −2) + ...
3!
xi yi xi yi
1 1.3 6 8.8
2 3.5 7 10.1
3 4.2 8 12.5
4 5.0 9 13.0
5 7.0 10 15.6
Self-Instructional Material 69
Interpolation y
16
14
NOTES
12
10
6
4
x
2 4 6 8 10
Interpolation requires a function that assumes the value of yi and xi. Each i
= 1, 2,...,10. From this graph, it appears that the actual relationship between x
and y is linear. However, it is likely that no line precisely fits the data, because of
errors in the data.
The least squares approach to this problem involves determining the best
approximating line when the error involved is the sum of the squares of the differences
between the y-values on the approximating line and the given y-values. Hence,
constants a0 and a1 must be found that minimize the total least squares error.
10 10
E = ∑ | yi − (a1xi + a0 ) | = ∑ ( yi − (a1xi + a0 ))2
2
i =1 i =1
where, ai xi + a0 denotes the ith value on the approximating line and yi is the ith
given y-value. E is the error.
The least squares method is the most convenient procedure for determining
best linear approximations, and there are also important theoretical considerations
that favour this method. The least squares approach puts substantially more weight
on a point that is out of line with the rest of the data but will not allow that point to
dominate the approximation.
The general problem of fitting the best least squares line to a collection of
data {( xi , yi )}im= 1 involves minimizing the total error
m
E = ∑ ( yi − (a1xi + a0 ))2
i =1
with respect to the parameters a0 and a1. For minimum, find the partial derivatives
of a1 and a0 and put them equal to zero.
m m
∑ ( yi − (a1xi + a0 )) = 2 ∑ ( yi − a1 xi + a0 ) (−1)
∂ 2
0=
∂a0 i =1 i =1
70 Self-Instructional Material
Interpolation
∂ m m
and 0= ∑ i 1i 0
∂a1 i =1
( y − ( a x + a )) 2
= 2 ∑ ( yi − a1xi − a0 ) (− xi )
i =1
These equations simplify to the normal equations
NOTES
m m m m m
a0 ∑ xi + a1 ∑ xi2 = ∑ xi yi and a0 . m + a1 ∑ xi = ∑ yi .
i =1 i =1 i =1 i =1 i =1
10(572.4) − 55(81)
and a1 = = 1.538.
10(385) − (55)2
Self-Instructional Material 71
Interpolation
So, P(x) = 1.538x – 0.360.
The problem of approximating a set of data, {(xi, yi) | t = 1, 2, ..., m}, with
an algebraic polynomial
NOTES
Pn(x) = an x n + an − 1 x n − 1 + ... + a1 x + a0
of degree n < m – 1 using least squares is handled in a similar manner. It requires
choosing the constants a0 , a1 , ..., an to minimize the total least squares error..
m
E = ∑ ( yi − Pn ( xi ))2 .
i =1
m m m m m
a0 ∑ x1i + a1 ∑ xi2 + a2 ∑ xi3 + ... + an ∑ xin + 1 = ∑ yi xi1,
i =1 i =1 i =1 i =1 i =1
M
m m m m m
a0 ∑ xin + a1 ∑ xin +1 + a2 ∑ xin + 2 + ... + an ∑ xi2n = ∑ yi xin .
i =1 i =1 i =1 i =1 i =1
The normal equations will have a unique solution, provided that the xi are
distinct.
Continuous Least Squares Approximation
Suppose f ∈ C| a, b | and you want a polynomial of degree at most n,
n
n
Pn(x) = an x + an −1x
n −1
+ ... + a1x + a0 = ∑ ak x k .
k =0
72 Self-Instructional Material
Interpolation
∂E
(a0 , a1 , ..., an ) = 0 for each j = 0, 1, ..., n.
∂a j
for each j = 0, 1, ..., n. Setting these to zero and rearranging, you obtain the
(n + 1) linear normal equations
n
∑ ak ∫a x j + k dx = ∫a x j f ( x) dx,
b b
for each j = 0, 1, ..., n,
k =0
which must be solved for the n + 1 unknown a0, a1, ..., an. The normal equations
have a unique solution provided that f ∈ C [a, b].
Example 2.24: Find the least squares approximating polynomial of degree 2 for
the function f(x) = sin πx on the interval [0, 1].
1 1 1 2 1
a0 ∫ 1dx + a1 ∫ xdx + a2 ∫0 x dx = ∫0 sin πx dx,
0 0
1 1 1 3 1
a0 ∫ xdx + a1 ∫ x 2 dx + a2 ∫ x dx = ∫0 x sin πx dx,
0 0 0
1 1 1 4 1 2
a0 ∫ x2dx + a1 ∫ x3dx + a2 ∫0 x dx = ∫0 x sin πx dx,
0 0
1 1 2 1 1 1 1 1 1 1 π2 − 4
a0 + a1 + a2 = , a0 + a1 + a2 = , a0 + a1 + a2 = .
2 3 π 2 3 4 π 3 4 5 π3
These three equations in three unknown can be solved to obtain
Self-Instructional Material 73
Interpolation Consequently, the least squares polynomial approximation of degree 2 for
2
f(x) = sin πx on [0, 1] is P2 ( x) = − 4.12251x + 4.12251x − 0.050465. (See
Figure 2.1)
NOTES y
1.0
0.8
0.6
P2(x)
0.4
0.2
x
0.2 0.4 0.6 0.8 1.0
b b j + k +1 − a j + k +1
∫a x j + k dx =
j + k +1
.
k ≠ j and j ≠ 0.
π π
∫−π φn + k ( x) φ j ( x) dx = ∫−π sin kx cos jx dx
NOTES
The trigonometric identity
1 1
sin kx cos jx = sin ( k + j ) x + sin ( k − j ) x
2 2
can now be used to give
1 π
∫−π φn + k ( x) φ j ( x) dx = 2 ∫−π (sin(k + j ) x + sin(k − j ) x) dx
π
π
1 − cos(k + j ) x − cos(k − j ) x
= − = 0,
2 k+ j k− j −π
since cos(k + j)π = cos(k + j) (–π) and cos(k – j)π = cos (k – j) (–π). The result
also holds when k = j, for in this case you have sin(k – j)x = sin 0 = 0.
Showing orthogonality for the other possibilities from {φ0 , φ1 , ..., φ2n − 1}
is similar and uses the appropriate trigonometric identities from the collection
1
sin jx cos kx = (sin( j + k ) x + sin( j − k ) x )
2
1
sin jx sin kx = (cos( j − k ) x + cos( j + k ) x )
2
1
cos jx cos kx = (cos( j + k ) x + cos( j − k ) x )
2
to convert the products into sums.
Given f ∈ C[–π, π], the continuous least squares approximation by functions
in Tn is defined by
2n − 1 2n − 1
1 1
Sn(x) = 2 0 ∑ k k
a + a φ ( x ) =
2
a0 + ∑ ak .cos kx
k =1 k =1
1 π
Where, ak =
π ∫−π f ( x) φk ( x) dx for each k = 0, 1, ..., 2n – 1.
π 2 n
(−1)k − 1
Sn(x) = +
2 π
∑ k2
cos kx.
k =1
y = |x|
4 4
y = S3(x) = cos x – cos 3x
2 9
y = S1(x) = S 2(x) = 4 cos x
2
2
y = S0(x) =
2
2 2
Figure 2.2 First few trigonometric polynomials for f(x) = | x |
Since | cos kx | ≤ 1 for every k and x, the series converges and S(x) exists
for all real numbers x.
2.7.1 Fourier Series
There is a discrete analog to Fourier series that is useful for the least squares
approximation and interpolation of large amounts of data when the data are given
76 Self-Instructional Material
at equally spaced points. Suppose that a collection of 2m paired data points Interpolation
{( x j , y j )}2j m= −0 1 is given, with the first elements in the pairs equally partitioning a
closed interval. For convenience, you may assume that the interval is [–π, π] and
that, NOTES
j
x j = −π + π for each j = 0, 1, ..., 2m – 1.
m
If this were not the case, a linear transformation could be used to change
the data into this form.
0 = xm
The orthogonality follows from the fact that if r and m are positive integers
with r < 2m, you have
2m − 1 2m − 1
∑ cos rx j = 0 and ∑ sin rx j = 0.
j=0 j=0
by setting to zero the partial derivatives of E with respect to the ak’s and the bk’s.
NOTES This implies that
1 2m − 1
ak = ∑ y j cos kx j ,
m j=0
for each k = 0, 1, ..., n
1 2m − 1
and bk = ∑ y j sin kx j ,
m j=0
for each k = 1, 2, ..., n – 1.
Let f(x) = x4 – 3x3 + 2x2 – tan x(x – 2). To find the discrete least squares
approximation S3 for the data {( x j , y j )}9j = 0 , where xj = j/5 and yj = f(xj), first
requires a transformation from [0, 2] to [–π, π]. The required linear transformation
is
zj = π(xj – 1),
and the transformed data are of the form
9
z j
z j , f 1 + .
π
j=0
Consequently, the least squares trigonometric plynomial is
2
+ a3 cos 3 z+ ∑ (ak cos kz + bk sin kz ),
a0
S3(z) =
2 k =1
1 9 zj
Where, ak =
5
∑ f 1 + cos kz j , for k = 0, 1, 2, 3
j=0 π
1 9 zj
and bk =
5
∑ f 1 + sin kz j , for k = 1, 2.
j=0 π
78 Self-Instructional Material
Table 2.10 lists values of f(x) and S3(x). Interpolation
Example 2.25: Find the value of log 323.5 from the following data:
x: 321 322.8 324.2 325
log10 x: 2.55051 2.50893 2.51081 2.51188
Self-Instructional Material 79
Interpolation Thus,
log 323.5 =2.50893 + (323.5 – 322.8) × (0.00134286) + (323.5 – 322.8)
× (323.5 – 324.2) × (–0.00000244) + (323.5 – 322.8) ×
(323.5 – 324.2) . (323.5 – 325) × (–0.00000022)
NOTES
= 2.50987
2.8 SUMMARY
In this unit, you have learned that:
• The process of estimating or approximating the value of independent variable
y for a given value of x in the specified range say xi ≤ x ≤ xn is known as
Interpolation. If the approximating function φ(x) is a polynomial, the process
is known as Polynomial Interpolation.
• If you want to estimate the value of y for any given value of x which is
outside the range xi ≤ x ≤ xn, then the process of approximating y = f(x)
by another function φ(x) is known as Extrapolation.
• Interpolation is very important in numerical analysis as it provides the base
for numerical differentiation and integration.
• The study of interpolation is based on the calculus of finite differences.
• Finite differences deal with the variations in a function corresponding to the
changes in the independent variables.
• y = f(x) is a function where x is the independent and y is the dependent
variable.
• Finite differences are the differences between the values of the function or
differences between the past differences. There are four types of differences
which are forward difference, backward difference, central difference and
divided difference.
• Central difference is denoted by δ and is defined as:
– δy1/2 = y1 – y0, δy3/2 = y2 – y1
δyn – 1/2 = yn – yn – 1 and so on.
80 Self-Instructional Material
• Higher order difference are obtained as: Interpolation
1. If you want to calculate the value of independent variable x for a given value
of dependent variable y, the process is known as inverse interpolation.
Self-Instructional Material 81
Interpolation 2. The differences y1 – y0, y2 – y1, y3 – y2, ..., yn – yn – 1 when denoted as
∇y1, ∇y2,..., ∇yn, respectively, are called the first backward differences.
3. Shift operator (E) is defined by the equation
NOTES Eyr = yr + 1 or Ef(x) = f(x + h)
which shows that the effect of E is to shift the functional value yr to the next
higher value yr + 1.
4. Newton’s formula is used for constructing the interpolation polynomial. It
makes use of divided differences.
5. Gauss’ forward difference formula is applicable when a lies between 0 and
1.
2
6. Gauss’s forward and backward formulas are used to derive Stirling’s
formula.
7. The merits of Lagrange’s formula are as follows:
• The formula is simple and easy to remember.
• There is no need to construct the divided difference table and you
can directly interpolate the unknown value with the help of given
observations.
8. The least square approach involves determining the best approximating line
when the error involved is the sum of the squares of the differences between
the y-values on the approximating line and the given y-values.
9. Trigonometric functions are used to approximate functions that have periodic
behaviour, functions with the property that for some constant
T, f(x+T) = f(x) for all x.
10. Fourier series are used to describe the solution of various ordinary and
partial-differential equations that occur in physical situations.
82 Self-Instructional Material
Long-Answer Questions Interpolation
= u x − ∆ 2u x −1 + ∆ 4u x − 2 − ∆ 6u x − 3 + ∆8u x −4 − ....
∞
1 ∞ 1 ∆ ∆2
(ii) ∑ ∑ x 4 2 4 u0 .
u2 x =
2 x=0
u + 1 − + − ...
x=0
8. If f(E) is a polynomial in E such that,
f(E) = a0 E n + a1E n −1 + a2 E n − 2 + ... + an
Prove that f(E) ex = ex f(e), taking the interval of differencing unity.
We now proceed to study the use of finite difference calculus for the purpose
of interpolation. This we shall do in three cases as follows:
(i) The value of the argument in the given data varies by an equal
interval. The technique is called an interpolation with equal
intervals.
Self-Instructional Material 83
Interpolation (ii) The values of argument are not at equal intervals. This is known
as interpolation with unequal intervals.
(iii) The technique of central differences.
NOTES 9. The following table gives the distance in nautical miles of the visible horizon
for the given heights in feet above the earth’s surface.
x: 100 150 200 250 300 350 400
y: 10.63 13.03 15.04 16.81 18.42 19.9 21.27
Use Newton’s forward formula to find y when x = 218 ft.
10. Use Newton’s divided difference formula to find f(7) if, f(3) = 24, f(5) =
120, f(8) = 504, f(9) = 720, and f(12) = 1716.
11. Find the least squares polynomials of degrees 1, 2, and 3 for the data in the
following table. Compute the error E in each case.
xi 1.0 1.1 1.3 1.5 1.9 2.1
yi 1.84 1.96 2.21 2.45 2.94 3.18
Mott, J.L. Discrete Mathematics for Computer Scientists, 2nd Edition. New
Delhi: Prentice-Hall of India Pvt. Ltd., 2007.
Bonini, Charles P., Warren H. Hausman, and Harold Bierman. Quantitative
Analysis for Business Decisions. Illinois: Richard D. Irwin, 1986.
Charnes, A., W.W. Cooper, and A. Henderson. An Introduction to Linear
Programming. New York: John Wiley & Sons, 1953.
Hogg, Robert. V., Allen T. Craig and Joseph W. McKean. Introduction to
Mathematical Statistics. New Delhi: Pearson Education, 2005.
84 Self-Instructional Material
Numerical Integration
AND DIFFERENTIAL
NOTES
INTERPOLATION
Structure
3.0 Introduction
3.1 Unit Objectives
3.2 Numerical Differentiation
3.3 Differentiating a Tabulated Function
3.4 Differentiating a Graphical Function
3.5 Numerical Integration Interpolation Formulae
3.5.1 Newton–Cotes Quadrature Formulae
3.5.2 Trapezoidal Rule (n = 1)
3.5.3 Simpson’s (1/3) Rule (n = 2)
3.5.4 Boole’s Rule (n = 4)
3.5.5 Weddle’s Rule (n = 6)
3.5.6 Gauss Quadrature Formulae
3.5.7 Errors in Quadrature Formulae
3.6 Approximation Theory
3.7 Summary
3.8 Key Terms
3.9 Answers to ‘Check Your Progress’
3.10 Questions and Exercises
3.11 Further Reading
3.0 INTRODUCTION
Sir Isaac Newton had proposed Interpolation formulae for forward and backward
interpolation. These are used for numerical differentiation. Such tools are widely
used in the field of engineering, statistics and other branches of mathematics.
Computer science also uses these concepts to find nearly accurate solution for
differentiation.
Forward interpolations
Sir Isaac Newton had proposed a formula for forward interpolation that bears his
name. It is expressed as a finite difference identity from which an interpolated
value, in between tabulated points using first value y0 with powers of forward
difference is used. Forward difference is shown by using ∆ which is known as
forward difference, operator. Forward difference is defined as the value obtained
by subtracting the present value from the next value. If initial value is y0 and next
value is y1 then ∆y0 = y1 – y0. In a similar way ∆2 is used. The ∆2y0 = ∆y1 – ∆y0.
Proceeding this way, you may write for first forward difference, second forward
difference and like wise of nth forward difference as follows:
∆y0 = y1 – y0 ∆2y0 = ∆y1 – ∆y0 .... ∆ny0 = ∆n – 1y1 – ∆n – 1y0
Taking this difference you may denote the next term (s) and thus,
y 1 = y0 + ∆y0 ⇒ y1 = (1 + ∆)y0
86 Self-Instructional Material
Here, 1 + ∆ shows a forward shift and a separate operator E, known as Numerical Integration
and Differential
forward shift operator is used as E = 1 + ∆. Now in the light of this fact, you may Interpolation
write y1 = Ey0 and y2 = Ey1 = E(Ey0) = E2y0. Proceeding this way, you may write
yn = En – 1 y0.
NOTES
Backward interpolations
Just as there is forward different operators for forward interpolations, there are
backward difference operator for backward difference interpolation. This is also
credited to Newton. In forward you think to next, but in backward you think of
the preceding term, i.e., the one earlier to it. Backward difference are denoted by
backward difference operator ∇ and is given as:
yn – 1 = yn – ∇yn ⇒ ∇yn = yn – yn – 1 and yn – 1 = (1 – ∇)yn
Just as in forward difference it was y0 in backward difference operator it is
yn.
Thus, ∇y1 = y1 – y0 ∇y2 = y2 – y1 ⇒ ∇2y2 = ∇y2 – ∇y1 and proceeding this
way you get, ∇nyn = ∇n – 1 yn – ∇n – 1 yn – 1
Relation between difference operators
E =1+∇ ⇒ ∇=E–1
∇(Ey0) = ∇(y1) = y1 – y0
Thus, (1 – ∇) Ey0 = Ey0 – ∇(Ey0) = y1 – ∇(y1) = y1 – (y1 – y0) = y0
or, (1 – ∇) (1 + ∆)y0= y0 which is true for all the terms of y, i.e., y0, y1, y2, .... yn.
Thus, (1 – ∇) (1 + ∆) = 1 and (1 – ∇)–1 = (1 + ∆) = E
and also, ∆ = (1 – ∇) –1 –1.
Self-Instructional Material 87
Numerical Integration (1) Newton’s forward difference interpolation formula
and Differential
Interpolation k ( k − 1) 2 k ( k − 1) ( k − 2) 3
y = y0 + k∆y0 + ∆ y0 + ∆ y 0 + ... ...(3.1)
2! 3!
NOTES where, k =
x−a
...(3.2)
h
Differentiating Equation (3.1) with respect to k, you get,
dy 2k − 1 2 3k 2 − 6 k + 2 3
= ∆ y0 + ∆ y0 + ∆ y0 + ... ...(3.3)
dk 2 6
Differentiating Equation (3.2) with respect to x, you get,
dk 1
= ...(3.4)
dx h
You know that,
dy dk 1 3k 2 − 6k + 2 3
= ∆y0 +
dy 2k − 1 2
= . ∆ y0 + ∆ y0 + ...
dx dk dx h 2 6
...(3.5)
dy
Equation (3.5) provides the value of at any x which is not tabulated.
dx
Equation (3.5) becomes simple for tabulated values of x in particular when x = a
and k = 0.
Putting k = 0 in Equation (3.5), you get,
dy
= ∆y0 − ∆ 2 y0 + ∆3 y0 − ∆ 4 y0 + ∆ 5 y0 − ... ...(3.6)
1 1 1 1 1
dx x = a h 2 3 4 5
Differentiating Equation (3.5) with respect to x, you get
d2y d dy d dy dk
2 = =
dx dx dx dk dx dx
1 2 3 6 k 2 − 18k + 11 4 1
= ∆ y 0 + ( k − 1) ∆ y 0 + ∆ y0 + ...
h 12 h
1 2 3 6 k 2 − 18k + 11 4
= 2 ∆ y 0 + ( k − 1) ∆ y 0 + ∆ y0 + ...
h 12
...(3.7)
Putting k = 0 in Equation (3.7), you get
d2y
= 2 ∆ 2 y0 − ∆ 3 y0 + ∆ y0 + ...
1 11 4
2 ...(3.8)
dx x = a h 12
Similarly, you get
d3y 1 ∆ 3 y − 3 ∆ 4 y + ...
3 = 3 0 0 ...(3.9)
dx x = a h 2
and so on.
88 Self-Instructional Material
Aliter: You know that, Numerical Integration
and Differential
E = ehD ⇒ 1 + ∆ = ehD Interpolation
∆ 2 ∆3 ∆ 4
∴ hD = log (1 + ∆) = ∆ –
+ − + ...
2 3 4 NOTES
1 1 2 1 3 1 4
⇒ D = ∆ − ∆ + ∆ − ∆ + ...
h 2 3 4
Similarly,
2
D = 2 ∆ − ∆ 2 + ∆3 − ∆ 4 + ...
2 1 1 1 1
h 2 3 4
= 2 ∆ 2 − ∆3 + ∆ 4 − ∆5 + ...
1 11 5
h 12 6
D3 = 3 ∆ 3 − ∆ 4
1 3
and
h 2
(2) Newton’s backward difference interpolation formula
k ( k + 1) 2 k ( k + 1) ( k + 2) 3
y = yn + k ∇yn + ∇ yn + ∇ yn + ... ...(3.10)
2! 3!
x − xn
where, k = ...(3.11)
h
Differentiating Equation (3.10) with respect to k, you get,
dy 2k + 1 2 3k 2 + 6k + 2 3
= ∇yn + ∇ yn + ∇ yn + ... ...(3.12)
dk 2 6
Differentiating Equation (3.11) with respect to x, you get,
dk 1
= ...(3.13)
dx h
dy dy dk
Now, = .
dx dx dx
1 2k + 1 2 3k 2 + 6k + 2 3
= ∇ yn + ∇ yn + ∇ yn + ...
h 2 6
...(3.14)
dy
Equation (3.14) provides the value of at any x which is not tabulated.
dx
At x = xh, you have k = 0
∴ Putting k = 0 in Equation (3.14), you get,
dy
= ∇yn + ∇ 2 yn + ∇3 yn + ∇ 4 yn + ...
1 1 1 1
...(3.15)
dx x = xn h 2 3 4
Differentiating Equation (3.14) with respect to x, you get,
d2y d dy dk
=
dx 2
dk dx dx
Self-Instructional Material 89
Numerical Integration
and Differential 1 2 3 6k 2 + 18k + 11 4
Interpolation = ∇ y + ( k + 1)∇ y + ∇ y + ...
h2
n n n
12
...(3.16)
NOTES Putting k = 0 in Equation (3.16), you get,
d2y
= 2 ∇ 2 yn + ∇3 yn + ∇ 4 yn + ...
1 11
2 ...(3.17)
dx x = x0 h 12
d3y 1 ∇3 y + 3 ∇ 4 y + ...
Similarly, you get, 3 = 3 ...(3.18)
n n
dx x = x0 h 2
and so on.
Formulae for computing higher derivatives may be obtained by successive
differentiation.
Aliter: You know that,
E –1 = 1 – ∇
e–hD = 1 – ∇
k ( k − 1) ( k − 1) ( k − 2) k −
NOTES 1
+ 2 ∆5 y
−2
5!
k ( k + 1) ( k + 2) ( k − 1) (k − 2) ( k − 3) ∆ 6 y−3 + ∆ 6 y−2
+ + ...
6! 2
...(3.27)
x−a
where, k = ...(3.28)
h
Differentiating Equation (3.27) with respect to k, you get,
3k 2 − 3k + 1
dy 2k − 1 ∆ 2
y −1 + ∆ 2
y0
2 ∆2 y
= ∆y0 + + −1
dk 2! 2 3!
4k 3 − 6k 2 − 2k + 2 ∆ 4 y−2 + ∆ 4 y−1 5k 4 − 10k 3 + 3k − 1 5
+ + ∆ y−2
4! 2 5!
6k 5 − 15k 4 − 20k 3 + 45k 2 + 8k − 12 ∆ 6 y−3 + ∆ 6 y−2
+ + ... ...(3.29)
6! 2
Differentiating Equation (3.28) with respect to x, you get,
dk 1
=
dx h
dy dy dk
Now, = .
dx dk dx
3k 2 − 3k + 1
2 2 2 ∆3 y
= ∆y0 +
1 2k − 1 ∆ y −1 + ∆ y0
+ −1
h 2! 2 3!
1 ∆ y−3 + ∆ y−2
6 6
1 5
+ ∆ y−2 + + ...
24 90 2
...(3.33)
and so on.
Self-Instructional Material 93
Numerical Integration
and Differential 3.3 DIFFERENTIATING A TABULATED FUNCTION
Interpolation
2k − 1 ∆ 2 y + 3k − 6k + 2 ∆3 y + ...
2
dy
= ∆y0 + 0 0 ...(3.35)
dk 2 6
dy
For maxima or minima, = 0. Hence, by equating Equation (3.35) to zero (i.e.,
dk
right hand side of Equation (3.35)) and retaining up to third difference, you have,
3k 2 − 6k + 2 3
∆y0 +
2k − 1 2
∆ y0
+ ∆ y0 = 0
2 6
1 ∆3 y k 2 + (∆ 2 y − ∆ 3 y )k + ∆y − 1 ∆ 2 y + 1 ∆3 y = 0 ...(3.36)
or 0 0 0 0 0 0
2 2 3
Which is a quadratic equation in k. Thus, substitute the values of ∆y0, ∆2y0 and
∆3y0 from the difference table and solve Equation (3.36) for k. The corresponding
values of x are given by x = x0 + kh, at which y is maximum or minimum.
Example 3.1: Find f ′(1.0), f ′′(1.0), f ′(2.0) and f ′′(2.0) from the following table:
x f(x) = y ∆y ∆2 y ∆3 y ∆4 y ∆5 y
1.0 0.0
0.1280 Forw
ard d
1.2 0.1280 0.298 irecti
on
0.4260 0.018
1.4 0.5540 0.316 0.06
0.7420 0.078 –0.1
1.6 1.2960 0.394 –0.04
1.1360 0.038
ion
1.8 2.4320 0.432 r d direct
wa
1.5680 Back
2.0 4.0000
94 Self-Instructional Material
Since you have to find f ′(x) and f ′′(x) at the start of the table, you will Numerical Integration
and Differential
differentiate Newton’s forward difference formula at x = x0. Interpolation
You know that,
dy 1 ∆y − 1 ∆ 2 y + 1 ∆ 4 y − 1 ∆5 y + ... NOTES
= 0 2
dx x = x0 h
0
3
0
4
0
d2y
and, 2 = 1 ∆ 2 y − ∆3 y + 11 ∆ 4 y − 5 ∆5 y
dx x = x0 h2 0 0
12
0
6
0
∆ y0 − ...
137 6
+
180
Thus,
= 10.4583
Now to find f ′(2.0) and f ′′(2.0) you will differentiate Newton’s backward
interpolation formula at x = xn. You get,
dy 1 1 2 1 3
dx x = xn = h ∇yn + 2 ∇ yn + 3 ∇ yn + ...
d2y
and, 2
1
= 2 ∇2 y + ∇3 y + 11 ∇4 y + 5 ∇5 y + 137 ∇6 y + ...
dx x = xn h n n
12
n
6
n
180
n
Thus,
1
1.5680 + × 0.432 + × 0.038 + × (−0.04) + × (−0.1)
1 1 1 1
f ′(2.0) =
0.2 2 3 4 5
= 8.8333
1
0.432 + 0.038 + × (−0.04) + × (−0.1)
11 5
and f ′′(2.0) = 2
(0.2) 12 6
= 8.7500
Self-Instructional Material 95
Numerical Integration Example 3.2: Obtain the value of f ′(0.04) using Bessel’s formula from the
and Differential
Interpolation following .
x : 0.01 0.02 0.03 0.04 0.05 0.06
f(x) : 0.1023 0.1047 0.1071 0.1096 0.1122 0.1148.
NOTES
x − x0
Solution: Using Bessel’s formula at x0 = 0.04 with h = 0.01 and k =
h
x − 0.04
=
0.01
[f ′(x)]x = 0.04 =
1 ∆y − 1 (∆ 2 y + ∆ 2 y ) + 1 ∆ 3 y
h 0 4 −1 0
12
−1
1 1 5
+ (∆ 4 y−2 + ∆ 4 y−1 ) − ∆ y−2
24 120
x k f(x) = y ∆y ∆2 y ∆3 y ∆4 y ∆5 y
0.01 –3 0.1023
0.0024
0.02 –2 0.1047 0.0000
0.0024 0.0001
0.03 –1 0.1071 0.0001 –0.0001
0.0025 0.000 0.0000
0.04 0 0.1096 0.0001 –0.0001
0.0026 –0.0001
0.05 1 0.1122 0.000
0.0026
0.06 2 0.1148
By using the above shown values in the Equation (1), you have,
1 1 1
[f ′(x)]0.04= 0.0026 − × (0.0001 + 0.000) + × (−0.001
0.01 4 12
1
+ × (−0.0001 + 0)
24
1
= [0.0026 – 0.000025 – 0.000008333 – 0.000004167]
0.01
1
= [0.0025625] = 0.25625
0.01
96 Self-Instructional Material
Example 3.3: Given are the following pairs of values of x and y: Numerical Integration
and Differential
Interpolation
x: 1 2 4 8 10
y: 0 1 5 21 27
NOTES
dy
Determine at x = 4 numerically.
dx
Solution: Since the values are not equispaced you will have to use Newton’s
divided difference formula. The difference table is shown below:
Difference Table for Newton’s Divided Difference Formula.
x y ∆y ∆2 y ∆3 y ∆4 y
1 0
1
2 1 1/3
2 0
4 5 1/3
4 –1/16 –1/144
8 21 –1/6
3
10 27
1
+ (x – 1) (x – 2) (x – 4) (x – 8) × −
144
x2 − 3x + 2
= 0 + (x – 1) + +0
3
x 4 − 15 x3 + 70 x 2 − 120 x + 64
–
144
Self-Instructional Material 97
Numerical Integration Now differentiating f(x) with respect to x, you have,
and Differential
2 x − 3 4 x3 − 45 x 2 + 140 x − 120
Interpolation
f ′(x) = 1 + +
3 144
NOTES Now f ′(x) at x = 4,
8 − 3 + 256 − 720 + 560 − 120
[f ′(x)]4 = 1 +
3 144
5 24
= 1+ − = 2.5
3 144
Example 3.4: From the following table find the minimum and maximum value of
the function.
x: 0 1 2 3 4 5
f(x): 0 0.25 0 2.25 16.0 56.25
Solution: The forward difference table is given below (as x are equispaced)
Forward Difference Table with Minimum and Maximum Value
0.25
1 0.25 − 0.5
–0.25 3.0
2 0 2.5 6
2.25 9.0 0
3 2.25 11.5 6
13.75 15.0
4 16.0 26.5
40.25
5 56.25
Newton’s forward difference formula is
k ( k − 1) 2
f(x) = f ( x0 ) + k ∆f ( x0 ) + ∆ f ( x0 ) + ... ...(1)
2!
Differentiating with respect to k, you have,
df ( x ) (2k − 1) 2 3k 2 − 6k + 2 3
= ∆f ( x0 ) + ∆ f ( x0 ) + ∆ f ( x0 ) + ...
dk 2 6
...(2)
98 Self-Instructional Material
Numerical Integration
df ( x ) and Differential
For maxima and minima, = 0. Interpolation
dk
Hence, equating the right-hand side of Equation (2) to zero and retaining
only upto third differences, you obtain NOTES
2k − 1 2 3k 2 − 6k + 2 3
∆f ( x0 ) + ∆ f ( x0 ) + ∆ f ( x0 ) = 0
2 6
1 ∆3 f ( x ) k 2 + ∆ 2 f ( x ) − ∆3 f ( x ) k
or, 0
2
0 0
1 2 1 3
+ ∆f ( x0 ) − ∆ f ( x0 ) + ∆ f ( x0 ) = 0
2 3
Substituting these values from Table, you have,
1 2 1 1
× 3 k + (−0.5 − 3)k + 0.25 − × (− 0.5) + × (3.0) = 0
2 2 3
3 2
or, k − 3.5k + 1.5 = 0
2
Solving for k, you get
3
−(−3.5) ± (3.5)2 − 4 × ×1.5
2
k =
3
2×
2
Self-Instructional Material 99
Numerical Integration
and Differential d2 f d2 f
Interpolation Now you will find . If 2 > 0 at k, then the value of the function
dx 2 dk
NOTES d2 f
at that k is minimum and if 2 < 0 at k, then the value of the function is
dk
maximum at that k. Thus,
d2 f 2 6k − 6 3 12k 2 − 36k + 22
2 = ∆ f ( x0 ) + ∆ f ( x0 ) + × ∆ 4 f ( x0 ) + ...
dk 6 24
at k = 1.7676,
d2 f 6 × 1.7676 − 6
2 = (–0.5) +
× (3) +
dk 6
12 × (1.7676)2 − 36 × (1.7676) + 22
× 6 = 0.7676
24
This is positive. Thus, the value of the function at x = 1.7676 is minimum.
Also at k = 0.5656
d2 f
= (–0.5) + (0.5656 – 1) × (3)
dk 2
12 × (0.5656) 2 − 36 × (0.5656) + 22
+ × 6
24
= –0.4339
d2 f
As the value 2 at k = 0.5656 is negative, hence the value of the function at
dk
x = x0 + kh = 0.5656 is maximum.
3
Graph of f(x)
2
1
x
0 1.5
0.5 1 2 2.5 3
a
3
Graph of f(x)
2
1
x
0 1.5
0.5 1 2 2.5 3
As we know, f ′(a) shows slope of the tangent at point (x, f(x)) on this
graph f(x). For sketching graph of the derivative of f(x), we select few points on
this graph. We take values of x at 0, 0.5, 1.0, 1.5, 2.0, 2.5 and 3.0 and make a
table by finding slopes of tangents at these points. The table is drawn as below:
m = f'(x) 3 0 –4 –3 0 1 0
NOTES
For measuring exact slope, use a ruler and a grid. But in the absence of a
ruler and a grid we may make a rough sketch. Now with these tabulated values a
graph can be made which is shown below.
y
3
Graph of f (x)
2
1
x
0 0.5 1 2 3
1.5 2.5
–1
–2
–3
–4
With the help of this graph of derivative, one can find derivative at any
point. Also, one may note that the graph intersects x-axis at points corresponding
to low and high points on graph of the function f(x). By taking more points, more
accurate graph can be drawn.
The process of evaluating a definite integral from a set of tabulated values (xi, yi)
of the integrand (or function) f(x) is known as quadrature, when applied to a
function of single variable.
Numerical integration of a function f(x) is done by representing f(x) as an
interpolation formula and then integrating it between the prescribed limits. By this
approach, a quadrature formula is obtained for approximate integration of the
function f(x).
3.5.1 Newton–Cotes Quadrature Formulae
b
Let I= ∫ f ( x) dx. ...(3.37)
a
y = f(x)
NOTES
y0 y1 y2 yn
O x0 x0 + h x0 + 2h X
x0 + nh
We divide the internal (a, b) into n-equal parts of width h such that,
a = x0, x1 = x0 + h, x2 = x0 + 2h, ..., xn = x0 + nh = b
Then,
x0 + nh
I= ∫ x0 f ( x) dx.
Put x = x0 + kh ⇒ dx = h.dk
When x = x0 ⇒ k = 0
and when x = xn ⇒ k = n
n n
so, I= ∫ f ( x0 + kh). h . dk = ∫ yk .dk
0 0
Replacing f (x0 + kh) or yk by Newton’s forward difference interpolation
formula, you obtain
n
k (k − 1) 2 k (k − 1) (k − 2) 3
I = h ∫ y0 + k .∆y0 + ∆ y0 + ∆ y0 + ... dk
0
2! 3!
n 1 n2 n 2 1 n3 2 3
I = nh y0 + ∆y0 + − ∆ y0 + − n + n ∆ y0
2 2 3 2 6 4
1 n 4 3n3 11n 2
+ − + − 3n ∆ 4 y0 + ... .
24 5 2 3
This is known as Newton–Cote’s quadrature formula. By substituting
n = 1, 2, 3, ... you will deduce the following important quadrature formulae.
f ( x) dx = h y0 + ∆y0
x0 + h 1
∫x
x1
0
f ( x) dx = ∫x
0 2
h
= ( y0 + y1 ) {As ∆y0 = y1 – y0}
2
Similarly from (x1, y1) to (x2, y2),
f ( x) dx = h y1 + ∆y1 = ( y1 + y2 )
x0 + 2 h 1
∫x + h
h
{As ∆y1 = y2 – y1}
0 2 2
and so on,
x0 + nh
∫x + (n − 1) h f ( x) dx =
h
( yn − 1 + yn )
0 2
Adding these n-integrals, you have,
x0 + nh
∫x0
f ( x) dx = h y0 + 2( y1 + y2 + ... + yn −1 ) + yn
2
This is known as trapezoidal rule.
1
3.5.3 Simpson’s Rule (n = 2)
3
By fitting a parabola through the points (x0, y0), (x1, y1) and (x2, y2), i.e., putting
n = 2 such that the differences of order higher than second vanish, you obtain,
x0 + 2 h 1 h
I= ∫x0 f ( x) dx = 2h y0 + ∆y0 + ∆ 2 y0 = [ y0 + 4 y1 + y2 ]
6 3
As ∆y0 = y1 − y0
2
∆ y0 = y2 − 2 y1 + y0
x0 + 4 h h
Similarly, ∫x0 +2h f ( x) dx =
3
[ y2 + 4 y3 + y4 ]
x0 + nh h
and so on ∫x0 +(n−2)h f ( x) dx =
3
[ y n − 2 + 4 yn − 1 + y n ]
k (k − 1) (k − 2) 3 k (k − 1) (k − 2) (k − 3) 4
+ ∆ y0 + ∆ y0 dk
3! 4!
(By Newton’s forward interpolation formula)
5 3 7
= 4h y0 + 2∆y0 + ∆ 2 y0 + ∆ 3 y0 + ∆ 4 y0
3 2 90
2h
= (7 y0 + 32 y1 + 12 y2 + 32 y3 + 7 y4 )
45
x0 +8h
Similarly, ∫x0 +4h f ( x) dx = 2h (7y4 + 32y5 + 12y6 + 32y7 + 7y8) and
45
so on.
Adding all these integrals from x0 to x0 + nh, where n is a multiple of 4, you
get,
x0 + nh 2h
∫x0 f ( x) dx =
45
[7 y0 + 32 y1 + 12 y2 + 32 y3 + 14 y4 +
32 y5 + 12 y6 + 32 y7 + 7 y8 + ...]
6 k (k − 1) 2 k (k − 1) (k − 2) 3
= h ∫ y0 + k ∆y0 + ∆ y0 + ∆ y0
0 2! 3!
k ( k − 1) ( k − 2) ( k − 3) 4 k ( k − 1) ( k − 2) ( k − 3) ( k − 4) 5
+ ∆ y0 + ∆ y0
4! 5!
k (k − 1) (k − 2) (k − 3) (k − 4) (k − 5) 6
+ + ∆ y0 dk
6!
k2 1 k3 k2 2 1 k4 3 2 3
= h ky0 + ∆y0 + − ∆ y0 + − k + k ∆ y0
2 2 3 2 6 4
1 k 5 3k 4 11k 3
+ − + − 3k 2 ∆ 4 y0
24 5 2 3
1 k6 5 35k
4
50k 3
+ − 2k + − + 12k 2 ∆5 y0
120 6 4 3
6
1 k 7 5k 6 225k 4 274k 3
+ − + 17 k 5 − + − 60k 2 ∆ 6 y0
720 7 2 4 3 0
9 2 3 41 4
= 6h y0 + 3∆y0 + ∆ y0 + 4∆ y0 + ∆ y0
2 20
11 5 41 6
+ ∆ y0 + ∆ y0
20 840
6h 2 3 4
= 20 y0 + 60∆y0 + 90∆ y0 + 80∆ y0 + 41∆ y0
20
41 6
+ 11 ∆5 y0 + ∆ y0
42
3h
= [ y0 + 5 y1 + y2 + 6 y3 + y4 + 5 y5 + y6 ]
10
Similarly,
x0 +12 h 3h
∫x0 +6h f ( x) dx =
10
[ y6 + 5 y7 + y8 + 6 y9 + y10 + 5 y11 + y12 ]
.
.
.
x0 + nh 3h
∫x0 +(n − 6)h f ( x) dx = 10 [ yn − 6 + 5 yn − 5 + yn − 4 + 6 yn − 3 + yn − 2
+ 5 yn −1 + yn ]
Adding the above integrals, you get,
x0 + nh 3h
∫x0 f ( x) dx =
10
[ y0 + 5 y1 + y2 + 6 y3 + y4 + 5 y5 + 2 y6
1
Use (i) Trapezoidal rule (ii) Simpson’s rule (iii) Weddle’s rule to find the
3
distance covered by the car in 12 minutes.
Solution: Let k be the distance travelled by the car in t min.
dk
then, = v, where v is the velocity of the car.
dt
12
so, [ k ]tt == 12
0 = ∫0 v. dt
Here n = 6. so by
12 h
∫0 v. dt =
3
[v0 + v6 + 2(v2 + v4 ) + 4(v1 + v3 + v5 )]
2
= [0 + 0 + 2(30 + 18) + 4(22 + 27 + 7)] = 213.33 km
3
(iii) By Weddle’s rule:
12 3h
∫0 v. dt =
10
[v0 + 5v1 + v2 + 6v3 + v4 + 5v5 + v6 ]
3× 2
= [0 + 5 × 22 + 30 + 6 × 27 + 18 + 5 × 7 + 0]
10
= 213.00 km
Example 3.6: A curve is given by the following table:
x: 0 1 2 3 4 5 6
y: 0 2 2.5 2.3 2 1.7 1.5
The x-coordinate of the area is bounded by the curve. The end coordinates
6
and the x-axis is given by Ax = ∫0 xy dx, where A is the area. Find x using
1
Simpson’s rule.
3
6
Solution: Given A.x = ∫0 xy dx, where A is the area
6
∫0 xy dx
∴ x = 6
∫0 y dx
108 Self-Instructional Material
Numerical Integration
6
as A = ∫0 y dx is the area now, you will compute xy as shown in Table given and Differential
Interpolation
below.
Computing Values for x, y NOTES
x: 0 1 2 3 4 5 6
y: 0 2 2.5 2.3 2 1.7 1.5
X = xy: 0 2 5.0 6.9 8 8.5 9.0
1
Then, using Simpson’s rule (with h = 1)
3
6
∫0 xy dx
h
= [ X 0 + X 6 + 2( X 2 + X 4 ) + 4( X1 + X 3 + X 5 )]
3
1
= [0 + 9 + 2(5 + 8) + 4(2 + 6.9 + 8.5)]
3
= 34.8667
6 h
and, ∫0 ydx =
3
[ y0 + y6 + 2( y2 + y4 ) + 4( y1 + y3 + y5 )]
1
= [0 + 1.5 + 2(2.5 + 2) + 4(2 + 2.3 + 1.7)]
3
= 11.500
Thus, x =
∫0 xy dx = 34.8667 3.0319
6
∫0 y dx 11.50
Example 3.7: A river is 60 ft wide. The depth ‘d’ infect at a distance x ft from one
bank is given by the following table:
x: 0 10 20 30 40 50 60
d: 0 5 9 15 10 6 2
Find approximately the area of the cross-section by, (i) Trapezoidal rule
1
(ii) Simpson’s rule (iii) Weddle’s rule.
3
Solution: The required area of the cross-section of the river is given as,
60
∫0 d .dx
NOTES
10
= [(0 + 2) + 2(5 + 9 + 15 + 10 + 6)] = 460
2
Hence, the required area of the cross-section of the river using Trapezoidal
rule is 460 sq. m.
1
(ii) By Simpson’s rule:
3
60
∫0
h
d .dx = [( d 0 + d 6 ) + 4( d1 + d3 + d5 ) + 2( d 2 + d 4 )]
3
10
= [(0 + 2) + 4(5 + 15 + 6) + 2(9 + 10)]
3
= 480
Hence, the required area of the cross-section of the river using Simpson’s
1
rule is 480 sq.m.
3
(iii) By Weddle’s rule:
60 3h
∫0 d .dx =
10
[ d0 + 5d1 + d 2 + 6d3 + d 4 + 5d5 + d 6 ]
3 × 10
= [0 + 5 × 5 + 9 + 6 × 15 + 10 + 5 × 6 + 2]
10
= 3[25 + 9 + 90 + 10 + 30 + 2] = 498.
Thus, the required area of the cross-section of the river using Weddle’s rule
is 498 sq. m.
3.5.6 Gauss Quadrature Formulae
We have already discussed the Trapezoidal rule and Simpson’s rule of integration
of a function. Now, we will discuss the Gauss Quadrature rule to obtain the formula
for polynomial having certain maximum power.
A general quadrature rule in numerical analysis is an approximation of the definite
integral of a function. It is also stated as a weighted sum of function values at
specified points within the domain of integration.
The Gauss quadrature formula or Gaussion quadrature rule is named after Carl
Friedrich Gauss and is defined as a n-point quadrature rule that yields the exact
result of the polynomials of degree 2n –1 for n points xi and n weight wi. The
110 Self-Instructional Material
domain of integration for the rule is taken as [a, b] or [–1, 1]. Hence, the Gauss Numerical Integration
and Differential
quadrature rule is stated as, Interpolation
1 n
∫ f ( x ) dx ≈ ∑ wi f ( xi )
–1
i =1 NOTES
The evaluation points are just the roots of a polynomial which belong to a class
of orthogonal polynomial.
Legendre polynomials are taken as the associated polynomials for the integration
of Pn(x). The nth polynomial is normalized to give Pn(1) = 1, the ith Gauss node,
xi, is the ith root of Pn; and its weight wi is given by,
2
wi = (1 – x 2 ) ( P ' ( x ))2
i n i
I= ∫ f ( x)dx
a
∫ f ( x)dx = ∫ (a
a
0 + a1 x + a2 x 2 + a3 x 3 ) dx
a
b
x2 x3 a4
= 0
4 a
a x + a1 + a 2 + a3
2 3
b2 – a2 b3 – a 3 b4 – a4
a
= 0 ( b − a ) + a1 2
+ a 3
+ a
2 3 4
(3.39)
b
∴ ∫ f ( x)dx
a
= ≈ c1 f ( x1 ) + c2 f ( x2 )
b−a
c2 = NOTES
2
b − a 1 b + a
x1 = – + 2
2 3
b − a 1 b + a
x2 = + 2
2 3
Substituting these values in the equation (9.9), we get,
b
∫ f ( x)dx
a
= ≈ c1 f ( x1 ) + c2 f ( x2 )
where P(x) is the polynomial representing the function y = f(x) in the given
interval [a, b].
Table 3.2 gives the error in each of the rule of numerical integration.
Table 3.2 Error in each Rule of Numerical Integration
( x − x0 )n n ( x − x0 )n + 1 ( n + 1)
+ ... + f ( x0 ) + ·f (k )
n! (n + 1)!
( x − x0 ) n + 1 f n + 1 (k )
The term is called the ‘remainder term’. Here, k is a
(n + 1)!
function of x and lies between x and x0, i.e., x ≤ k ≤ x0.
( n + 1) ( x − x0 ) n + 1 | ( x − x0 ) n + 1 | NOTES
Truncation Error = f ( k ) . ≤ M
( n + 1)! (n + 1)!
Example 3.8: Give a Taylor series representation of f(x) = sin x and compute sin
x correct to three significant digits.
Solution: The Taylor series representation of sin x is (about x0 = 0)
( x − x0 ) 2
f(x) = f ( x0 ) + f ′( x0 ) ( x − x0 ) + f ′′( x0 ) + ...
2!
x2 x3
= sin (0) + cos (0) ( x) – sin (0). – cos (0) + ...
2! 3!
x3 x 5 x 7
Thus, sin x = x – + − + ...
3! 5! 7!
x3 x5 x7
or, sin x = x − + − + ... ...(1)
6 120 5040
Since you require sin x correct to three significant digits and the Truncation
error after three terms of equation is
1
T.E. ≤ = 0.000198
5040
Thus, you truncate the series after three terms.
x3 x5
So, sin x = x − +
3! 5!
is the required representation to obtain the value of sin x correct to three significant
digits.
Example 3.9: Express the Taylor Series expansion of
x 2 x3 x 4
e− x = 1 − x + − + + ....
2! 3! 4!
in terms of Chebyshev polynomials.
1 1
+ [3T0 ( x ) + 4T2 ( x ) + T4 ( x )] − [10T1 ( x ) + 5T3 ( x ) + T5 ( x )] + ...
195 1920
3.7 SUMMARY
1. Sir Isaac Newton had proposed Interpolation formulae for forward and
backward interpolation. These are used for numerical differentiation.
2. It is expressed as a finite difference identity from which an interpolated
value, in between tabulated points using first value y0 with powers of forward
difference is used. Forward difference is shown by using ∆ which is known
as forward difference, operator. Forward difference is defined as the value
obtained by subtracting the present value from the next value.
3. Backward difference is denoted by backward difference operator ∇ and is
given as:
yn – 1 = yn – ∇yn ⇒ ∇yn = yn – yn – 1 and yn – 1 = (1 – ∇)yn
4. In forward difference it was y0 whereas in backward difference operator it
is yn.
5. Central difference operator is defined as δf(x) = f(x + h/2) – f(x – h/2).
This is known as first central difference operator.
k ( k − 1) 2 k ( k − 1) ( k − 2) 3
6. y = y0 + k∆y0 + ∆ y0 + ∆ y 0 + ...
2! 3!
x−a
where, k=
h
k ( k + 1) 2 k ( k + 1) ( k + 2) 3
7. y = yn + k ∇yn + ∇ yn + ∇ yn + ...
2! 3!
x − xn
where, k =
h
8. Formulae for computing higher derivatives may be obtained by successive
differentiation.
x−a
where, k=
h
10. Newton’s divided difference formula is,
f(x) = y0 + (x – x0). ∆y0 + (x – x0) (x – x1) ∆2y0
+ (x – x0) (x – x1) (x – x2) ∆3y0 + (x – x0) (x – x1)
× (x – x2) (x – x3) ∆4y0 + ...
d2 f
11. If 2 > 0 at k, then the value of the function at that k is minimum and
dk
d2 f
if 2 < 0 at k, then the value of the function is maximum at that k.
dk
12. The process of evaluating a definite integral from a set of tabulated values
(xi, yi) of the integrand (or function) f(x) is known as quadrature, when
applied to a function of single variable.
13. Numerical integration of a function f(x) is done by representing f(x) an
interpolation formula and then integrating it between the prescribed limits.
By this approach, a quadrature formula is obtained for approximate
integration of the function f(x).
14. Newton Cote’s quadrature formula is as follows:
n 1 n2 n 2 1 n3 2 3
I = nh y0 + ∆y0 + − ∆ y0 + − n + n ∆ y0
2 2 3 2 6 4
1 n 4 3n3 11n 2
+ − + − 3n ∆ 4 y0 + ... .
24 5 2 3
15. The trapzoidal is as follows:
x + nh
0 h
f ( x) dx = y2 + 2 ( y1 + y2 + ... + yn – 1 ) + yn
x 0 2
x0 + nh 2h
16. Boole’s rule ∫ x0 f ( x) dx =
45
[7 y0 + 32 y1 + 12 y2 + 32 y3 + 14 y4 +
32 y5 + 12 y6 + 32 y7 + 7 y8 + ...]
17. Weddle’s rule
x0 + nh 3h
∫ x0 f ( x) dx =
10
[ y0 + 5 y1 + y2 + 6 y3 + y4 + 5 y5 + 2 y6
Short-Answer Questions
1. What is trapezoid rule?
1
2. What is Simpson’s rule?
3
3. What is Weddle’s rule?
21 1
4. Evaluate ∫1
x
dx by Simpson’s
3
rule with n = 4 and determine the error
by direct integration.
Long-Answer Questions
1. The population of a certain town is shown in the table below:
6. From the table given below, find for what value of x, y is maximum. Also
find this value of y.
x: 3 4 5 6 7 8 9
y: 0.205 0.240 0.259 0.262 0.250 0.224 0.201
7. Using the following data find the value of x for which y is minimum. Find this
value of y also.
x: −1 1 2 3
y: − 21 15 12 3
Estimate the area bounded by the curve, x-axis and the lines x = 1 and
x = 4.
4
[Hint: Area A = ∫1 y . dx ]
x3 1
17. Estimate the length of the arc of the curve y = from (0, 0) to 1,
3 3
1
using Simpson’s rule taking 8 subintervals.
3
2
dy 1
[Hint: Length of the arc S = ∫0 1 + dx . dx ]
18. Express the following as polynomials in x
1 1
(i) T0 ( x) + 2T1 ( x) + T2 ( x) (ii) 2T0 ( x) − T2 ( x) + T4 ( x)
4 8
(iii) 5T0 ( x) + 2T1 ( x) + 4T2 ( x) − 8T3 ( x)
19. Compute cos x correct to 3 significant digits. Obtain a series with minimum
number of terms using Taylors series and Chebyshev series.
20. Obtain a Taylor series expansion of ex and represent it in terms of Chebyshev
polynomials.
Mott, J.L. Discrete Mathematics for Computer Scientists, 2nd Edition. New
Delhi: Prentice-Hall of India Pvt. Ltd., 2007. NOTES
Chance, William A. Statistical Methods for Decision Making. Illinois: Richard
D Irwin, 1969.
Bonini, Charles P., Warren H. Hausman, and Harold Bierman. Quantitative
Analysis for Business Decisions. Illinois: Richard D. Irwin, 1986.
Charnes, A., W.W. Cooper, and A. Henderson. An Introduction to Linear
Programming. New York: John Wiley & Sons, 1953.
Hogg, Robert. V., Allen T. Craig and Joseph W. McKean. Introduction to
Mathematical Statistics. New Delhi: Pearson Education, 2005.
AND TRANSCENDENTAL
NOTES
EQUATIONS
Structure
4.0 Introduction
4.1 Unit Objectives
4.2 Types of Non-linear Equations
4.3 Intermediate Value Theorem
4.4 Methods of Finding Solutions of Algebraic and Transcendental Equations
4.4.1 Direct Methods
4.4.2 Iterative Method
4.4.3 Initial Guess
4.4.4 Rate of Convergence
4.5 Bisection Method
4.5.1 Minimum Number of Iterations Required in Bisection Method to Achieve
the Desired Accuracy
4.5.2 Convergence of Bisection Method
4.6 Regula–Falsi Method or Method of False Position
4.6.1 Convergence of Regula–Falsi Method
4.7 Secant Method
4.8 Fixed Point Iteration Method
4.8.1 Condition for Convergence of Iteration Method
4.8.2 Convergence of Iteration Method
4.9 Newton–Raphson Method
4.9.1 Geometrical Interpretation
4.9.2 Convergence of Newton–Raphson Method
4.9.3 Newton–Raphson Method for System of Non-Linear Equations
4.9.4 Generalized Newton’s Method for Multiple Roots
4.10 Summary
4.11 Key Terms
4.12 Answers to ‘Check Your Progress’
4.13 Questions and Exercises
4.14 Further Reading
4.0 INTRODUCTION
In this unit, you will learn that the solution of algebraic and transcendental equations
is a topic of interest not only to mathematicians but also to scientists and engineers.
There are analytical methods to solve algebraic equations of degree two or three
or four. However, the need often arises to solve higher degree algebraic equation
or transcendental equations for which direct or analytical methods are not available.
Thus, numerical techniques or approximate methods are required to solve such
Self-Instructional Material 125
Solution of Algebraic and kind of equations and before you proceed to solve such kind of equations you
Transcendental Equations
should find out the difference between algebraic and transcendental equations.
Let f(x) = a0 + a1 x + a2 x 2 + ... + an x n be a polynomial of degree n, where
NOTES a’s are constants such that a0 ≠ 0 and n is a positive integer.
The polynomial f(x) = 0 is called an algebraic equation of degree n. If f(x)
contains some other functions such as log x, ex, trigonometric functions or ax, etc.,
then f(x) = 0 is called a transcendental equation.
You will learn that the solution of these equations is nothing but to obtain a
value of x which satisfies f(x) = 0. This value is called the root of the equation f(x)
= 0. In other words, x = α is called a root of equation f(x) = 0 if f(α) = 0. The
process of finding the roots of an equation is known as the solution of that equation.
This is a problem of great importance in applied mathematics. To find out the roots
numerically you should make use of the fundamental theorem on roots of f(x) = 0
in a ≤ x ≤ b.
Quintic Equations
A quintic equation is formed as a polynomial function of degree five, equated to
zero and has the form:
ax5 + bx4 + cx3 + dx2 + ex + f = 0,
Where a, b, c, d, e, f are constants that may real or complex numbers where
a ≠ 0. These contain polynomial of odd degree and hence, their nature has some
similarity with cubic function and their graph may have additional local maxima
and local minima. There are some quintic equations such as x5 – x – 1 = 0 which
is not exactly solvable and these can be analyzed qualitatively analyzed.
Sextic Equations
A sextic equation contains polynomial of degree six and has the form:
Logarithmic Equations
A logarithmic equation involves one or more terms having logarithmic function. A
logarithmic function of x can be expanded as an expansion in terms of different
powers of x and hence these are non-linear. Exapnsion of logarithmic function is
shown as
x 2 x3 x 4 xn
log (1 + x) = x − + – + – ... + (– 1) n – 1 + Rn
2 3 4 n
Putting x – 1 in place of x, we get expansion for log x. Thus,
log x = (x – 1) – (x – 1)2/2 + (x – 1)3/3 – (x – 1)4/4 + …..
Exponential Equations
If f(x) is continuous in a closed interval [a, b] and f(a), f(b) have opposite signs,
i.e., f(a) . f(b) < 0, then the equation f(x) = 0 has atleast one root between x = a
and x = b.
This method is based on the repeated application of the Intermediate Value Theorem.
Let f(x) be a continuous function between a and b, such that f(a) . f(b) < 0.
For definiteness, let f(b) > 0 and f(a) < 0. Then the first approximation to the root
α is given by ,
a +b
x0 = i.e., average of the ends of the range
2
Now, if f(x0) = 0, x0 is the root of f(x) = 0. Otherwise the root will lie
between x0 and a or between x0 and b depending upon the sign of f(x0).
Continue the same process of bisecting the interval until the root is found to
desired accuracy (see Figure 4.2).
1
xi + 1 = ( xi − 1 + xi )
2
Now, xi + 1 is called the root of the equation f(x) = 0 if f(xi + 1) = 0 or
|(xi + 1) – (xi)| ≤ ε where ‘ε’ is the desired accuracy. Clearly, ε is a very small
positive number.
Example 4.1: Find the root of the equation x3 – 3x – 7 = 0, using Bisection
method, correct upto 3 decimal places.
Solution: Let f(x) = x3 – 3x – 7 for initial guess.
x: 0 1 2 3
f(x): –7 –9 –5 11
Sign: –ve –ve –ve +ve
Since f(2) is –ve and f(3) is +ve, a root lies between 2 and 3 by Intermediate
Value Theorem.
So, a = 2, b = 3, the first approximation to the root is,
a+b 2+3
x0 = = = 2.5
2 2
Now f(2.5) = (2.5)3 – 3(2.5) – 7 = 1.125 which is +ve. So, 2.5 is not the
root of f(x) = 0 and now the root lies between 2 and 2.5 as f(2) . f(2.5) < 0. The
second approximation is given as,
2 + 2.5
x1 = = 2.25
2
1
S.No. a b x= (a + b) f(a) f(b) f(x)
2
3. 2.375 2.5 2.4375 –ve +ve +ve
interval
4. 2.375 2.4375 2.40625 –ve +ve –ve
interval
5. 2.40625 2.4375 2.4219 –ve +ve –ve
interval
6. 2.4219 2.4375 2.4297 –ve +ve +ve
new interval
7. 2.4219 2.4297 2.4258 –ve +ve –ve
a1 + b1
and x2 = ...(4.1)
2
where, x2 is the next approximation to the root α.
Now if x2 is the desired root of f(x) = 0, then f(x2) = 0. Otherwise either
f(x2) . f(a1) < 0 or f(b1) . f(x2) < 0, giving the interval (a1, x2) or (x2, b1) for the
root 2. After renaming the interval as (a2, b2), you have,
1
b2 – a2 = (b1 − a1 ) ...(4.2)
2
By substituting the value of (b1 – a1) from Equation (4.1) to Equation (4.2),
you have,
1
(b2 – a2) = (b0 − a0 )
22
Proceeding in the same way, you find,
an + bn
xn + 1 =
2
Which gives the (n + 1)th approximation to the root and root lie in (an, bn).
1
where, (bn – an) = (b0 − a0 ) ...(4.3)
2n
Here, (bn – an) gives the length of the interval in which root lies after
(n + 1)th approximation. Minimum length of the interval gives the maximum
accuracy in the root, if the desired accuracy is ∈, then,
| bn – an | ≤ ε
b0 − a0
⇒ ≤ε
2n
ei + 1
Comparing this with ≤ k, you have p = 1. Hence, Bisection method
eip
is linearly convergent.
Example 4.2: Using Bisection method, compute the smallest positive root of the
equation xex – 2 = 0, correct upto three decimal places.
Solution: Let f(x) = xex – 2
f(0) = –2, i.e., negative and f(1) = 0.7183, i.e, positive.
Thus, the root exist between 0 and 1.
You know that Bisection method is,
xn −1 + xn
xn + 1 =
2
Here xn – 1 = 0 and xn = 1
x0 + x1 0 + 1
so, x2 = = = 0.5.
2 2
and, f(x2 ) = –1.1756, i.e., negative.
Thus, the root exists between (0.5, 1).
0.5 + 1
Also, x3 = = 0.75
2
Self-Instructional Material 135
Solution of Algebraic and and f(x3 ) = –0.4122, i.e., negative.
Transcendental Equations
So root lies between (0.75, 1).
Next iteration is,
NOTES 0.75 + 1
x4 = = 0.875
2
and f(x4 ) = 0.875.e0.875 – 2 = 0.0990, i.e, positive.
So, the next interval for the root is (0.75, 0.875) and hence,
0.75 + 0.875
x5 = = 0.8125
2
and f(0.8125) = – 0.1690, i.e, negative.
Next interval (0.8125, 0.875).
So, next iteration is,
0.8125 + 0.875
x6 = = 0.8438
2
and f(x6 ) = –0.0381, i.e, negative.
Next interval is (0.8438, 0.875).
0.8438 + 0.875
So, x7 = = 0.8594
2
and f(x7 ) = 0.0297, i.e, positive.
Next interval for root is (0.8438, 0.8594).
0.8438 + 0.8594
and x8 = = 0.8516
2
∴ f(x8 ) = –0.0044, i.e, negative.
Next interval is (0.8516, 0.8594).
0.8516 + 0.8594
and x9 = = 0.8555
2
∴ f(x9 ) = 0.0126, i.e, positive.
So, the next interval is (0.8516, 0.8555).
0.8516 + 0.8555
and x10 = = 0.8536
2
∴ f(x10) = 0.0042, i.e, positive.
Next interval is (0.8516, 0.8536)
0.8516 + 0.8536
and x11 = = 0.8526
2
136 Self-Instructional Material
∴ f(x11) = –0.000024, i.e, negative. Solution of Algebraic and
Transcendental Equations
Giving the next interval for root is (.8526, .8536).
0.8526 + 0.8536
and x12 = = 0.8531 NOTES
2
∴ f(x12) = 0.00215, i.e., positive.
Next interval is (0.8526, 0.8531).
0.8526 + 0.8531
and x13 = = 0.8529
2
∴ f(x13) = 0.0012, i.e., positive.
Since the two interations x12 and x13 are approximately same as
|(x12 – x13)| = 0.0002.
So, the desired root of the given equation is 0.8529 correct upto 3-decimal
places.
y A (a, f(a))
x1 x0 b x
O a
C B(b, f(b))
(x0,f(x0))
The equation of the chord joining the two points A(a, f(a)) and B(b, f(b)) is
given as,
a f (b) − bf (a )
or x0 =
f (b) − f (a )
which is the new approximation to the root. Now if f(a) and f(x0) are of opposite
signs, the root will lie between (a, x0) otherwise the root will lie in (x0, b). Now if
f(a) . f(x0) < 0, you will replace b by x0 and obtain the new approximation x1. If
f(x0) . f(b) < 0, than you will replace a by x0 and obtain the next approximation
x1 .
In general, if the root lies between (xi – 1, xi), i.e, f(xi – 1) . f(xi) < 0, then the
next approximation xi + 1 to the root ‘α’ is obtained as,
xi −1. f ( xi ) − xi f ( xi −1 )
xi + 1 =
f ( xi −1 ) − f ( xi )
Continue the same process till f(xi + 1) = 0 or | (xi + 1 – xi) | < ε, where ε is
a very small positive number known as desired accuracy. Usually the value of ε is
taken as,
1
ε = × 10−3
2
Example 4.3: Using Regula–Falsi method compute the smallest positive root of
the equation xex – 2 = 0 correct upto three decimal places.
Solution: f(x) = xex – 2
For initial guess
x: 0 1
Sign of f(x): –ve +ve
Since f(0) < 0 and f(1) > 0, the root for f(x) = xex – 2 lies between (0, 1).
i.e., f(0) = –2, f(1) = 0.7183
Since x0 = 0 and x1 = 1, the next approximation x2 is obtained as,
xi −1. f ( xi ) − xi . f ( xi −1 )
xi + 1 =
f ( xi ) − f ( xi −1 )
( xi − xi −1 )
xi + 1 = xi − . f ( xi ) ...(4.5)
f ( xi ) − f ( xi −1 )
Now, let α be the root of f(x) = 0 and let ei – 1, ei and ei + 1 be the errors
associated with (i – 1)th, ith and (i + 1)th iterations.
Then, xi + 1 = α + ei + 1, xi = α + ei and xi – 1 = α + ei – 1
By Equation (4.5) you have,
[α + ei ) − (α + ei − 1 )]
(α + ei + 1) = (α + ei ) − × f (α + ei )
f (α + ei ) − f (α + ei − 1 )
(ei − ei − 1 )
or ei + 1 = ei − × f (α + ei ) ...(4.6)
f (α + ei ) − f (α + ei − 1 )
By expanding f(α + ei) and f(α + ei – 1) using Taylor’s series upto third
term, you get,
(ei − ei − 1 )
ei + 1 = ei − × f (α + ei )
e2
f (α) + ei f ′(α) + i f ′′(α) + ... −
2!
e 2
( f (α) + ei − 1 f ′(α) + i − 1 f ′′(α) + ...
2!
After simplifying the denominator, you have,
e2
(ei − ei − 1 ) f (α) + ei f ′(α) + i f ′′(α) + ...
2!
ei + 1 = ei −
2 2
ei − ei −1
(ei − ei − 1 ) f ′(α) + f ′′(α)
2
ei2
ei f ′(α) + f ′′(α)
ei + 1 = ei − 2!
f ′(α) + i ei −1 ) f ′′(α)
( e +
2
ei2 f ′′(α )
ei +
2 f ′(α )
or ei + 1 = ei −
(ei + ei −1 ) f ′′(α)
1 + f ′(α)
2
–1
ei2 f ′′(α ) (ei + ei −1 ) f ′′(α )
= ei − ei + 1+
2 f ′(α ) f ′(α )
or ei + 1
2
Using Binomial expansion and retaining the first two terms, you get,
e2 f ′′(α) ei + ei −1 f ′′(α)
ei + 1 = ei − ei + i . 1− .
2 f ′(α) 2 f ′(α)
f ′′(α )
ei + 1 ~ ei ei − 1 .
f ′(α )
or ei + 1 = ei ei − 1 . M ...(4.8)
ei + 1
≤k
eip
or ei + 1 = eip . k ...(4.9)
or ei = eip−1. k
1/ p
e
or ei – 1 = i
k
Substitute the values of ei + 1 and ei – 1 in equation (4.8), you have,
1/ p
e
k . eip = ei × i .M
k
1 1
− 1 + 1+
or eip = M ×K P × ei p ...(4.10)
Comparing the powers of ei on both the sides of (4.10), you have,
1
p =1+ or p2 – p – 1 = 0
p
which on solving gives,
1± 5
p=
2
Taking only the positive sign, you have,
p = 1.1618
which is the rate of convergence of Regula–Falsi method.
In Secant method, you do not require the condition of f(a) . f(b) < 0. Thus, this
method is an improvement over Regula–Falsi method. Appart from this, Secant
method is just like Regula–Falsi method.
In this method also the graph of y = f(x) is approximated by a secant
line but at each iteration the two most recent approximations are used to find
142 Self-Instructional Material
out the next approximation. Thus, in this method, it is not required that the Solution of Algebraic and
Transcendental Equations
interval which you take for finding out the next approximation to the root must
contain the root.
The formula for Secant method is just the same as in case of Regula–Falsi
NOTES
method.
xi −1 f ( xi ) − xi . f ( xi −1 )
xi + 1 =
f ( xi ) − f ( xi −1 )
where (xi – 1, xi) is the most recent interval. Thus, the condition f(xi – 1) . f(xi) < 0,
is not required in this method.
Convergence of Secant method is 1.1618 as in case of Regula–Falsi method.
Example 4.4: Find a real root of the equation x log10x = 1.2 by Secant method,
correct upto three decimal places.
Solution: f(x) = x log10x – 1.2
Let the initial approximations be x0 = 1 and x1 = 2, then f(x0) = –1.2 and
f(2) = – 0.59794.
So, by Secant method, you have,
x1 − x0
x 2 = x1 − f ( x1 )
f ( x1 ) − f ( x0 )
2 −1
= 2− × (−1.2) = 3.9932
(−0.59794 + 1.2)
and f(x2 ) = 1.2012.
Next interval is (2, 3.9932).
x2 − x1
So, x 3 = x2 – . f ( x2 )
f ( x2 ) − f ( x1 )
3.9932 − 2
= 3.9932 – (1.2012)
1.2012 + 0.59794
= 2.6624
and f(x3 ) = – 0.0678
Next interval is (3.9932, 2.6624).
x3 − x2
So, x 4 = x3 − . f ( x3 ) = 2.7335
f ( x3 ) − f ( x2 )
and f(x4 ) = – 0.0062.
2.7407 − 2.7335
= 2.7407 − × (−0.00005)
(0.00005 + 0.0062)
= 2.7406
Since | x5 – x6 | = 0.0001, the value 2.7406 is correct upto three decimal
places and hence the desired root of x log10x – 1.2 is 2.7406.
Let f(x) = 0 be the equation. To find out the root of this equation by successive
approximations, rewrite f(x) = 0 in the form of ,
x = φ(x) ...(4.11)
The roots of f(x) = 0 are the same as the points of intersection of the straight
line y = x and the curve x = φ(x).
Let x = x0 be the initial approximation to the root α (see Figure 4.4). Then,
the first approximation x1 is given by,
x 1 = φ(x0)
x
y =
y
x
O x0 x2 x3 x 1
Figure 4.4 Iteration Method
1 1
Where as φ′2(x) = + × < 1 for (1, 2)
3 (5 − x) 2/3
x 1 = φ( x0 ) = 3 5 − (1.5) = 1.5183
x 2 = φ( x1 ) = 1.5156
x 4 = φ( x3 ) = 1.5160.
Since x3 and x4 are same. The root of the equation x3 + x – 5 = 0 is NOTES
1.5160.
(ii) f(x) = ex – 4x, f(0) = 1, i.e., positive
f(1) = e′ – 4 = –1.2817, i.e., negative.
Thus, the root lies between (0, 1).
Let,
x 0 = 0.5
ex
The equation can be rewritten as x = ,
4
ex ex
then, φ(x) = , φ′( x) = and hence | φ′( x ) | < 1 for (0, 1).
4 4
Thus, you apply Iteration method, to get,
ex
xn + 1 = φ(xn) with φ(x) =
4
and x 0 = 0.5
e0.5
then, x 1 = φ(x0) = = 0.4122
4
e0.4122
x 2 = φ(x1) = = 0.3775
4
e0.3775
x 3 = φ(x2) = = 0.3647
4
e0.3647
x 4 = φ(x3) = = 0.3600
4
x 5 = φ(x4) = 0.3583
x 6 = φ(x5) = 0.3577
x 7 = φ(x6) = 0.3575
x 8 = φ(x7) = 0.3574
Since x7 and x8 are the same. The root of equation ex = 4x is 0.3574
correct upto three decimal places.
ei2
(α + ei + 1) = φ(α ) + ei φ′(α) + φ′(α ) + ...
2!
Since ei is small, neglecting the second and higher order terms of ei, you
have,
(α + ei + 1) = φ(α) + ei . φ′(α)
or ei + 1 = ei φ′(α ) + (φ(α ) − α ) ...(4.14)
Now, since α is the root of f(x) = 0, therefore α = φ(α) and Equation
(4.14) reduces to,
ei + 1 = eiφ′(α)
ei + 1
or = φ′(α)
ei
Let y = f(x) or f(x) = 0 be the equation and let x = x0 be the initial approximation
to the root. NOTES
Let x = x1 be the exact root, then f(x1) = 0 and x1 = x0 + h.
Where, h is a small positive number known as step size.
⇒ f(x0 + h) = 0 (as f(x1) = 0)
Expanding by Taylor’s theorem, you have,
h2 h3
f ( x0 ) + hf ′( x0 ) + f ′′( x0 ) + f ′′′( x0 ) + ... = 0
2! 3!
Since h is a very small number, so, you can neglect the terms involving h2,
3
h , ... .
So, after neglecting the higher order terms, you have,
f ( x0 ) + hf ′( x0 ) ~ 0
f ( x0 )
or h =–
f ′( x0 )
f ( x0 )
Thus, x 1 = x0 + h = x0 –
f ′( x0 )
f ( x0 )
so, x1 = x0 − gives the bett er appro ximation to the root of
f ′( x0 )
f(x) = 0 than x0.
Similarly, if x2 denotes a better approximation than x1, you have,
f ( x1 )
x 2 = x1 −
f ′( x1 )
Continue in this manner, you have
f ( xn )
xn + 1 = xn −
f ′( xn )
This is known as Newton–Raphson method.
4.9.1 Geometrical Interpretation
Let x0 be the initial approximation near to the root α of f(x) = 0. Then, the equation
of tangent drawn at point A(x0, f(x0)) on the curve y = f(x) is given as,
y – f(x0) = (x – x0) . f′(x0) ...(4.15).
This cuts the x-axis at point x1 and y = 0.
root α B
x
x2 x1 x0
y = f(x)
f (α + ε n )
or εn + 1 = ε n − ...(4.17)
f ′(α + ε n )
ε 2n f ′′(α )
f (α ) + ε n . f ′(α ) + + ...
εn + 1 = ε n − 2! ...(4.18)
ε n2
f ′(α ) + ε n . f ′′(α ) + f ′′′(α ) + ...
2!
Since α is the exact root of f(x) = 0, thus f (α) = 0. Then from Equation
(4.18) you have,
ε 2n f ′′(α )
f (α ) + ε n . f ′(α ) + + ...
εn + 1 = εn − 2!
ε2
f ′(α ) + ε n . f ′′(α ) + n f ′′′(α ) + ...
2!
ε 2n ε3
f ′′(α ) + n f ′′′(α) + ...
or εn + 1 = 2 3
ε2
f ′(α ) + ε n f ′′(α) + n f ′′′(α ) + ...
2!
Since εn is very-very small, neglecting the third and higher power of εn,
you have,
ε 2n f ′′(α)
εn + 1 =
2[ f ′(α) + ε n f ′′(α )]
ε 2n f ′′(α )
~ .
2 f ′(α)
That is, the error at (n + 1)th iteration is proportional to the square of the
error at nth iteration.
f ( xn ). f ′′( xn )
which gives that <1
[ f ′( xn )]2
Thus, Newton–Raphson method will converge if,
Solution: Let, x = 6
125
then, x 6 = 125
or x6 – 125 = 0
let f(x) = x6 – 125
Apply Newton–Raphson method, you have,
f ( xn )
xn + 1 = xn − f ′( x )
n
f ′(x) = 6x5
5(2.5)6 + 125
Let x0 = 2.5, then, x1 = = 2.2967
6(2.5)5
5(2.2967)6 + 125
x2 = = 2.2399
6(2.2967)5
5(2.2399)6 + 125
x3 = = 2.2361
6(2.2399)5
5(2.2361)6 + 125
x4 = = 2.23606
6(2.2361)5
~ 2.2361 (Correct up to four decimal places)
Since x3 and x4 are same, the root of f(x) = x6 – 125 is 2.2361
or 6
125 = 2.2361.
1
Example 4.8: Obtain the Iterative formula for a , 3 a and using
a
Newton–Raphson method.
Solution: (i) For a:
Let, x = a then, x2 – a = 0.
Also, let, f(x) = x2 – a
Then, according to Newton–Raphson method, you have,
f ( xn )
xn + 1 = xn – ...(1)
f ′( xn )
NOTES xn2 + a 1 a
or xn + 1 = = xn +
2 xn 2 xn
1 a
xn + 1 = xn +
2 xn
(ii) For 3 3
a : Let, x = a
⇒ 3
x –a =0
Let, f(x) = x3 – a
Then, f(xn) = xn3 − a and f′(xn) = 3 xn2 . Substituting these values in
Equation (1), you have,
xn3 − a
xn + 1 = xn −
3xn2
2 xn3 + a
or xn + 1 =
3xn2
1 a
or xn + 1 = 2 xn + 2 .
3 xn
1
Let, x=
a
1
or x2 – =0
a
1
Also let, f(x) = x2 –
a
1
then, f(xn ) = x2n – and f′(xn) = 2xn
a
1 1
or xn + 1 = xn + .
2 axn
1
which is the required Iterative formula for
using Newton–Raphson method.
a
Example 4.9: Use Newton–Raphson method to find the smallest positive root of
the equation tan x = x.
Solution: Given, tan x = x
or tan x – x = 0
Let, f(x) = tan x – x
∴ f ′(x) = sec2x – 1
and f (–1) = – 0.5574 i.e., negative.
and f(1) = +0.5574 i.e., positive.
This root lies between (–1, 1).
Let, x0 = 1
then Newton–Raphson formula is,
f ( xn )
xn + 1 = xn −
f ′( xn )
0.5574
=1–
+0.8508
= 0.3449
tan (0.3449) − 0.3449
Also, x 2 = 0.3449 –
sec 2 (0.3449) − 1
= 0.3449 – 0.0127
Since the error between x15 and x14 is | 0.0009 – 0.0013 | = 0.0004. The
smallest positive root is 0.0009 correct upto three decimal places.
4.9.3 Newton–Raphson Method for System of Non-Linear Equations
Let f1(x, y) = 0 and f2(x, y) = 0 be two non-linear equations.
Let (x0, y0) be the initial approximation to the solution of the system and let
(x1, y1) be the exact solution.
Then, x 1 = x0 + h and y1 = y0 + k
where h and k are very small numbers known as step size.
Then, you must have,
f1(x1, y1) = 0 and f2(x1, y1) = 0
Expanding using Taylor’s series, you have,
∂f1 ∂f
f1 ( x0 + h, y0 + k ) = f1 ( x0 , y0 ) + h + k 1 + ...
∂x0 ∂y0
∂f 2 ∂f
and f 2 ( x0 + h, y0 + k ) = f 2 ( x0 , y0 ) + h + k 2 + ...
∂x0 ∂y0
∂f 2 ∂f
and h + k 2 = –f2(x0, y0) ...(4.21)
∂x0 ∂y0
∂f1 ∂f ∂f ∂f
Thus, = 5.08; = –5.08; 2 = 5.08 and 2 = 5.08
∂x0 ∂y0 ∂x0 ∂y0
Newton–Raphson formula for set of non-linear equations is,
∂f1 ∂f
h + k 1 = –f1(x0, y0) ...(1)
∂x0 ∂y0
∂f 2 ∂f
and h + k 2 = –f2(x0, y0) ...(2)
∂x0 ∂y0
x 1 = x0 + h and y1 = y0 + k.
By putting the values in Equations (1) and (2), you have,
5.08h – 5.08k = –(–3)
5.08h + 5.08k = +3.0968
then, h = 0.2363 and k = 0.0095
Thus, x 1 = 2.54 + 0.2363 = 2.7763 and y1 = 2.5495
Let the next approximation be (x2, y2) such that x2 = x1 + h and
y2 = y1 + k.
So, f1(x1, y1) = –1.7921
f2(x1, y1) = –1.7922
∂f1 ∂f ∂f
= 5.5526; 1 = –5.099; 2 = 5.5526
∂x1 ∂y1 ∂x2
this is called the generalized Newton’s method when m =1, it reduces to,
f ( xn )
xn + 1 = xn −
f ′( xn )
Short-Answer Questions
1. Write intermediate value theorem in brief.
2. What is the sighest power of x in a quintic equation expressed in terms of x.
3. What are methods to find solution to algebraic and transcendental equation?
4. What is bisection method?
5. What is the equation used for computing by Newton –Raphson method?
Long-Answer Questions
1. Find a root of the following equations, using the Bisection method correct
upto three decimal places:
(i) x log10x = 1.2 (ii) 2x3 + x2 – 20x + 12 = 0
(iii) x3 + x – 1 = 0 (iv) cos x = 3x – 1
2. Using Regula–Falsi method find the real root of the following equations
correct to four decimal places:
(i) 2x – logx = 6 (ii) x6 – x4 – x3 – 1 = 0
(iii) x3 – 4x – 9 = 0 (iv) xex = 2
3. Evaluate 15 by Iteration method correct upto four decimal places.
4. Find the roots of the following equations correct to three decimal places
using Secant method:
1
(i) x – e–x = 0 (ii) x = + sin x
2
5. Apply Iteration method to find the negative root of the equation
x3 – 4x + 7 = 0 correct upto four decimal places.
6. Find the root of the equation 2x = cos x + 3 correct upto three decimal
places using (i) Iteration method and (ii) Regula–Falsi method.
1 1
7. Use the following formula for to find correct upto four decimal
a 5
162 Self-Instructional Material
Solution of Algebraic and
(3 − axn2 ) xn Transcendental Equations
places: xn + 1 =
2
Show that the convergence is of order two.
8. Apply Newton–Raphson method to find roots of the following equations: NOTES
(i) x = 1 + sin x (ii) x4 + x2 – 80 = 0
(iii) x3 + x2 – 100 = 0 (iv) ex cos x = 2.5
9. Solve the given non-linear equations using Newton–Raphson method.
(i) x2 – y2 + x = 5 (ii) x2 + y2 = 4
(iii) 2xy + y = 5 (iv) ex + y = 1
10. Find the cube root of 41 correct upto four decimal places.
11. Solve x4 – 5x3 + 20x2 – 40x + 60 = 0, by Newton–Rephson method to
find a pair of complex roots.
12. Use Newton–Raphson method to find the smallest root of the equation
ex sin x = 1.
13. Find the negative root of the equation x3 – 21x + 3500 = 0 correct upto
three decimal places using Newton–Raphson method.
14. Using Newton–Raphson method find a root of the following equations
correct upto three decimal places:
(i) ex = x3 + cos 25x which is near 4.5
(ii) x sin x + cos x = 0 which is near π
(iii) x4 – x – 10 = 0
15. Find a double root of the equation x3 – x2 – x + 1 = 0 using Newton–
Raphson method. correct up to three decimal places only.
Mott, J.L. Discrete Mathematics for Computer Scientists, 2nd edition. New
Delhi: Prentice-Hall of India Pvt. Ltd., 2007.
Bonini, Charles P., Warren H. Hausman, and Harold Bierman. 1986. Quantitative
Analysis for Business Decisions. Illinois: Richard D. Irwin.
Charnes, A., W.W. Cooper, and A. Henderson. 1953. An Introduction to Linear
Programming. New York: John Wiley & Sons.
Gupta, S.C., and V.K. Kapoor. 1997. Fundamentals of Mathematical Statistics.
9th Revised Edition, New Delhi: Sultan Chand & Sons.
TO ORDINARY
NOTES
DIFFERENTIAL EQUATIONS
Structure
5.0 Introduction
5.1 Unit Objectives
5.2 Picard’s Method of Successive Approximations
5.3 Taylor’s Series Method
5.4 Euler’s Method
5.5 Runge–Kutta Method
5.5.1 Runge–Kutta Second Order Method
5.5.2 Runge–Kutta Fourth Order Method
5.6 Predictor–Corrector Methods
5.6.1 Modified Euler’s Method
5.6.2 Milne’s Predictor–Corrector Method
5.7 Summary
5.8 Key Terms
5.9 Answers to ‘Check Your Progress’
5.10 Questions and Exercises
5.11 Further Reading
5.0 INTRODUCTION
In this unit, you will learn that an ordinary differential equation is an essential tool
for modelling many physical problems in science and engineering. Also many
scientific laws are more readily expressed in terms of rate of change.
There are analytical methods to solve differential equations but these methods
are applicable to a limited class of differential equations. Thus, numerical methods
are useful in solving all these kinds of differential equations.
The first order differential equation, is as follows:
dy
= f(x, y) given y(x0) = y0
dx
The solution of differential equation is the function that satisfies the differential
equation and also certain initial conditions on the function. Here, some numerical
methods to find the solution of differential equations will be discussed.
dy
Let, = f(x, y); y(x0) = y0 ...(5.1)
dx
It is an ordinary differential equation of first order.
Integrating Equation (5.1) from x0 to x with respect to x, we get,
y x
∫ dy = ∫ f ( x, y ) dx
y0 x0
x
or, y – y0 = ∫ f ( x, y ) dx
x0
or, y = y0 + ∫ f ( x, y ) dx ...(5.2)
x0
Let y1(x) be the first approximation of the solution of Equation (5.1), then
the second approximation y2(x) can be obtained as
x
y2(x) = y0 + ∫ f ( x, y1 ( x)) dx
x0
x
Similarly, yn + 1(x) = y0 + ∫ f ( x, yn ( x)) dx
x0
Example 5.1: Find the successive approximate solution of the initial value problem
dy
= xy + 1, with y(0) = 1, by Picard’s method.
dx
166 Self-Instructional Material
Numerical Solution to
dy
Solution. We have, = xy + 1 or y′ = xy + 1 with y(0) = 1 Ordinary Differential
dx Equations
x2
x
y x = 1 + ∫ ( x + 1) dx = 1 + x +
(1)
0
2
x
x2
= 1 + ∫
x 1 + x +
2
+ 1 dx
0
x
x3 x 2 x3 x 4
= 1 + ∫
2
x + x + + 1 dx = 1 + x + + +
0
2 2 3 8
x 2 x3 x 4 x5 x 6
=1+x+ + + + +
2 3 8 15 48
Example 5.2: Find four successive approximate solutions for the following initial
value problem: y′ = x + y, with y(0) = 1, using Picard’s method. Hence compute
y(0.1), y(0.2) and y(0.3), correct upto four significant digits.
Solution. We have, y′ = x + y with y(0) = 1
The first approximation is given by,
x
y x = y(0) + ∫ [ x + y(0)] dx
(1)
x2
x
= 1 + ∫ (1 + x) dx = 1 + x +
0
2
0
2 6
NOTES
x
x3 x3 x 4
y (x) = 1 + ∫ x + 1 + x + x + dx = 1 + x + x + +
(3) 2 2
0
6 3 24
x
x3 x 4
y x = 1 + ∫
2
and (4) 1 + 2 x + x + + dx
0
3 24
x3 x 4 x5
= 1 + x + x2 + + +
3 12 120
Hence, the values of y(0.1), y(0.2) and y(0.3) are ,
x2
Solution. Given that y′ = with y(0) = 0. Using Picard’s successive
1+ y 2
approximation method the first approximation is obtained as,
x2 x
y x = y(0) + ∫
(1)
2 dx
0
1 + y (0)
x2 x3
x
= 0+∫ . dx =
0
1+ 0 3
x
x2
y x = 0 + ∫
(3)
(2) 2
dx
0 1 + [ y x ]
x
x2
= ∫ 2
dx which is not integrable.
0 −1 x3
1 + tan
6
So, in computing y(0.5) and y(1.0) we will make use of y(2)x only.
−1 (0.5)
3
Hence, y (0.5) = tan
(2)
~ 0.0416
3
13
and y(2)(1.0) = tan–1 ~ 0.3218
3
Example 5.4: Find three successive approximate solutions of the following initial
value problem, by Picard’s method. Hence, compute y(0.2) and y(0.4) correct
upto four decimal places.
y′ = x2 – y2, y(0) = 1
Solution. The first three approximations are given as,
x3 x3
x
y x = 1 + ∫ ( x 2 − 1) dx = 1 +
(1)
− x or 1 − x +
0
3 3
2x
x3 x 2 x3 x 4
y x = 1 + ∫ x − 1 + x − dx = 1 − x +
(2) + −
0
3 2 3 12
2x
x 2 x3 x 4
and y x = 1 + ∫ x −1 + x −
(3) − + dx
0
2 3 12
NOTES x 2 x3 x 4 x5
= 1− x + + − +
2 6 12 60
Hence,
(0.2) 2 (0.2)3 (0.2) 4 (0.2)5
y(0.2) = 1 – (0.2) + + − + ~ 0.8212
2 6 12 60
dy
Consider, = f(x, y); y(x0) = y0
dx
Expanding this function using Taylor series around x = x0, we get:
( x − x0 )2
y(x) = y ( x0 ) + ( x − x0 ) y′( x0 ) + y′′( x0 ) + ... .
2!
where y′(x0), y′′(x0) are the first order and second order derivatives and so on.
Then, at x = x0 + h, the expansion becomes
h2
y(x0 + h) = y(x0) + h . y′(x0) + y′′( x0 ) + ... .
2
which gives the value of x at x0 + h, where h is the step size.
Note: Taylor series method can be applied to those functions which are not
complicated enough to calculate the derivatives.
Example 5.5: From the Taylor series solution of the initial value problem, y′ =
x2 – y2, y(0) = 0 up to seven terms and hence compute y(0.5) and y(1.0), correct
upto four decimal places.
Solution. Given that y′ = x2 – y2 with y(0) = 0.
Differentiating successively we get,
clearly, y′(0) = 0
Now, y″ = 2x – 2y.y′
= 2x – 2y . (x2 – y2)
Hence, y″(0) = 2.(0) – 2(0).(02 – 02) = 0
Similarly,
y′″ = 2 – 2[y.y″ + (y′)2]
170 Self-Instructional Material
so y′′′(0) = 2 Numerical Solution to
Ordinary Differential
iv Equations
yiv = –2(yy″′ + 3y′y″) and y (0) = 0
y v = −2( yy iv + 4 y′y′′ + 3( y′′) 2 ) and y v (0) = 0. NOTES
y0
O X0 X 1 X2 XN X
x0 x0 + h x0 + 2h x0 + nh
Example 5.9: Using Runge–Kutta method of order 4, compute y(0.2) and y(0.4)
dy
from = x2 + y2, y(0) = 1 taking h = 0.1.
dx
dy
Solution. Given that = x2 + y2 ≡ f(x, y)
dx
Taking h = 0.1, calculating y(0.2) in two steps as follows:
Step 1: x0 = 0, y0 = 1, h = 0.1
so, k 1 = hf(x0, y0) = 0.1(02 + 12) = 0.1
h k1 0.1 2 0.1 2
k 2 = hf 0
x + , y0 + = 0.1 + 1 + = 0.1105
2 2 2 2
h k2 2 0.1105
2
k 3 = hf x0 + , y0 + = 0.1 0.05 + 1 +
( ) = 0.1116
2 2 2
k 4 = hf(x0 + h, y0 + k3) = 0.1[(0.1)2 + (1 + 0.1116)2] = 0.1246
h k1
k 2 = hf x1 + , y1 +
2 2
2
0.1 0.1245
2
= 0.1 0.1 + + 1.1115 + = 0.1400
2 2
h k2
k 3 = hf x1 + , y1 +
2 2
2
0.1 0.1400
2
= 0.1 0.1 + + 1.1115 + = 0.1418
2 2
k 4 = hf ( x1 + h, y1 + k3 )
2 2
= 0.1[(0.1 + 0.1) + (1.1115 + 0.1418) ] = 0.1611
1
k = ( k1 + 2k 2 + 2k3 + k4 ) = 0.1415
6
giving y(0, 2) = y1 + k = 1.1115 + 0.1415 = 1.2530.
Now to calculate y(0.4) by taking.
x 0 = 0.2, y0 = 1.2530 and h = 0.1.
In two steps as follows:
Step 1: x0 = 0.2, y0 = 1.2530, h = 0.1
k 1 = hf(x0, y0) = 0.1[(0.2)2 + (1.2530)2] = 0.1610
h k1
k 2 = hf x0 + , y0 +
2 2
2
0.1 0.1610
2
= 0.1 0.2 + + 1.2530 + = 0.1841
2 2
h k2
k 3 = hf x0 + , y0 +
2 2
h k1 2 2
k 2 = hf x1 + , y1 + = 0.1[(0.35) + (1.5478) ] = 0.2518
2 2
h k2 2 2
k 3 = hf x1 + , y1 + = 0.1[(0.35) + (1.5656) ] = 0.2574
2 2
k 4 = hf(x1 + h, y1 + k3) = 0.1[(0.4)2 + (1.6971)2] = 0.3040
1
∴ k = ( k1 + 2k 2 + 2k3 + k 4 ) = 0.2565
6
giving y(0.4) = y1 + k = 1.4397 + 0.2565 = 1.6962.
Thus,
x = 0.1 y = 1.1115
x = 0.2 y = 1.2530
x = 0.3 y = 1.4397
and x = 0.4 y = 1.6962.
Example 5.10: Find by Runge–Kutta method of order 4 an approximate value of
dy
y at x = 0.8, given that y = 0.41 when x = 0.4 and = ( x + y) .
dx
Solution. Given x0 = 0.4, y0 = 0.41, h = 0.4
f(x, y) = ( x + y)
h k
k 2 = hf x0 + , y0 + 1
2 2
Self-Instructional Material 177
Numerical Solution to
Ordinary Differential 0.4 0.36
Equations = 0.4 0.4 + + 0.41 + = 0.4363
2 2
NOTES h k
k 3 = hf x0 + , y0 + 2
2 2
k 4 = hf ( x0 + h, y0 + k3 )
= 0.4 ( 0.8 ) + ( 0.8533) = 0.5143
1
k = ( k1 + 2k 2 + 2k3 + k4 ) = 0.4389
6
giving y(0.8) = y0 + k = 0.41 + 0.4389 = 0.8489.
dy
Example 5.11: Given that (y2 – 2x)/(y2 + x) and y = 1 at x = 0; find y for
dx
x = 0.1, 0.2, 0.3 and 0.4.
Solution. Runge–Kutta fourth order method is,
k 1 = hf(x0, y0)
h k1
k 2 = hf x0 + , y0 +
2 2
h k2
k 3 = hf x0 + , y0 +
2 2
and k 4 = hf(x0 + h, y0 + k3)
1
k = ( k1 + 2k 2 + 2k3 + k4 )
6
Thus, x 1 = x0 + h and y1 = y0 + k
y2 − 2x
Given that f(x, y) = , x0 = 0, y0 = 1 and h = 0.1
y2 + x
12 − 2 × 0
Step 1: k 1 = hf(x0, y0) = 0.1 2 = 0.1
1 +0
h k (1.05)2 − 2 × 0.05
k 2 = hf x0 + , y0 + 1 = 0.1 = 0.087
2 2 2
(1.05) + 0.05
178 Self-Instructional Material
Numerical Solution to
h k (1.0435)2 − 2 × 0.05 Ordinary Differential
k 3 = hf x0 + , y0 + 2 = 0.1 = 0.0868 Equations
2 2 2
(1.0435) + 0.05
1
k = ( k1 + 2k 2 + 2k3 + k4 ) = 0.0892
6
y1 = y0 + k = 1 + 0.0892 = 1.0892
giving y(0.1) = 1.0892 at x = 0.1
Step 2: x 1 = 0.1 y1 = 1.0892 and h = 0.1
(1.0892) 2 − 2 × 0.1
k 1 = hf(x1, y1) = 0.1 2 = 0.0767
(1.0892) + 0.1
h k (1.1275)2 − 2 × 0.15
k 2 = hf x1 + , y1 + 1 = 0.1 = 0.0683
2 2 2
(1.1275) + 0.15
h k (1.1234) 2 − 2 × 0.15
k 3 = hf x1 + , y1 + 2 = 0.1 = 0.0681
2 2 2
(1.1234) + 0.15
(1.1233)2 − 2 × 0.1
k 4 = hf(x1 + h, y1 + k3) = 0.1 2 = 0.0780
(1.1233) + 0.1
1
k = ( k1 + 2k 2 + 2k3 + k 4 ) = 0.0713
6
so, y(0.2) = y1 + k = 1.0892 + 0.0713 = 1.1605
Step 3: x2 = 0.1, y2 = 1.1605 and h = 0.1
(1.1605) 2 − 2 × 0.2
k 1 = hf(x2, y2) = 0.1 = 0.0612
(1.1605) 2 + 0.2
(1.1911) 2 − 2 × 0.25
k 2 = hf x2 + , y2 + 1 = 0.1
h k
= 0.0551
2 2 (1.1911) 2 + 0.25
h k (1.1880) 2 − 2 × 0.25
k 3 = hf x2 + , y2 + 2 = 0.1 = 0.0549
2 2 2
(1.1880) + 0.25
NOTES 1
k = ( k1 + 2k2 , 2k3 + k4 ) = 0.0548
6
so, y(0.3) = y2 + k = 1.1605 + 0.0548 = 1.2153
Step 4: x3 = 0.3, y3 = 1.2153 and h = 0.1
(1.2153) 2 − 2 × 0.3
k 1 = hf ( x3 , y3 ) = 0.1 2 = 0.0493
(1.2153) + 0.3
h k (1.24)2 − 2 × 0.35
k 2 = hf x3 + , y3 + 1 = 0.1 = 0.0444
2 2 2
(1.24) + 0.35
h k (1.2375)2 − 2 × 0.35
k 3 = hf x3 + , y3 + 2 = 0.1 = 0.0442
2 2 2
(1.2375) + 0.35
(1.2595) 2 − 2 × 0.4
k 4 = hf ( x3 + h, y3 + k ) = 0.1 = 0.0396
2
(1.2595) + 0.4
1
so, k = ( k1 + 2k 2 + 2k3 + k4 ) = 0.0444
6
giving y(0.4) = y3 + k = 1.2597.
Example 5.12: Using Runge–Kutta fourth order method find y for x = 0.1 and
0.2 given that y′ = xy + y2, y(0) = 1, taking h = 0.1. Given that y(0) = 1
Solution. Given that f(x, y) = xy + y2, x0 = 0, y0 = 1 and h = 0.1.
Step 1:
k 1 = hf(x0, y0) = 0.1[0 × 1 + 12] = 0.1
h k
k 2 = hf x0 + , y0 + 1 = 0.1[(0.05 × 1.05) + (1.05) 2 ] = 0.1155
2 2
h k
k 3 = hf x0 + , y0 + 2 = 0.1[(0.05 × 1.0578) + (1.0578) 2 ] = 0.1172
2 2
k 4 = hf(x0 + h, y0 + k3) = 0.1[(0.05 × 1.1172) + (1.1172)2] = 0.1304
1
k = ( k1 + 2k2 + 2k3 + k 4 ) = 0.1160
6
giving y(0.1) = y0 + k = 1 + 0.1160 = 1.1160
180 Self-Instructional Material
Step 2: Numerical Solution to
Ordinary Differential
x 1 = 0.1, y1 = 1.1160 and h = 1 Equations
k
k 3 = hf x1 + , y1 + 2 = 0.1 (0.15 ×1.1950) + (1.1950) = 0.1607
h 2
2 2
k 4 = hf(x1 + h, y1 + k) = 0.1[(0.2 × 1.2767) + (1.2767)2] = 0.1885
1
k = (k1 + 2k2 + 2k3 + k4 ) = 0.1602.
6
Giving y(0.2) = y1 + k = 1.1160 + 0.1602 = 1.2762
Hence, y(0.1) = 1.1160 at x = 0.1
and y(0.2) = 1.2762 at x = 0.2
The methods discussed so far are single step methods. Now you will learn about some
multi step methods which make use of the past information of the curve to extrapolate
the solution curve. Such methods are known as Predictor–Corrector methods.
5.6.1 Modified Euler’s Method
As in Euler’s method, the curve of solution in the interval X0X1 is approximated by
a target at S (See Figure 5.1) such that at S1, we have
y1 = y0 + hf(x0, y0) ...(5.3)
Now the slope of the curve of solution through S1 is,
dy
= f(x0 + h, y1)
dx S1
and the tangent at S1 to S1T1 is drawn meeting the ordinate through X2 in S2
(x0 + 2h, y2). We can find a better approximation y1(i) of y(x0 + h) by taking the
slopes of the curve as the mean of the slopes of the tangents at S and S1, thus,
Self-Instructional Material 181
Numerical Solution to
h
Ordinary Differential y11 = y0 + [ f ( x0 , y0 ) + f ( x0 + h, y1 )] ...(5.4)
Equations 2
Since the slope at S1 is not known, we compute the value of y1 from Euler’s
NOTES method (i.e. Equation (5.3) and substitute it to the R.H.S of Equation (5.4) to
obtain the first modified value of y1(1) .
h (1)
Similarly, y1(2) = y0 + [ f ( x0 , y0 ) + f ( x0 + h, y1 )]
2
repeat this step till two consecutive values of y agree.
Once y1 is obtained to a desired degree of accuracy, compute,
y2 = y1 + hf(x0 + h, y1)
and a better approximation of y2, i.e. y2(1) is obtained as,
h
y2(1) = y1 + [ f ( x0 + h, y1 ) + f ( x0 + 2h, y2 )]
2
repeat this step till y2 is obtained with desired accuracy.
The modified Euler’s method can be rewritten as:
ynp + 1 = yn + hf ( xn , yn ) ...(5.5)
h p
and ync + 1 = yn + [ f ( xn , yn ) + f ( xn + 1, yn + 1 )] ...(5.6)
2
Equation (5.5) is known as Predictor formula and Equation (5.6) is known
as Corrector formula.
dy
Example 5.13: Solve = 1 – y, y(0) = 0 in the range 0 ≤ x ≤ 0.3 using
dx
modified Euler’s method by choosing step size 0.1. Also compare the obtained
result with the exact solution.
Solution. Given that y′ = 1 – y, y(0) = 0.
Modified Euler’s method is
ynp + 1 = yn + hf(xn, yn)
h
and ync + 1 = yn + f ( xn , yn ) + f ( xn + 1 , ynp + 1 ) .
2
Taking h = 0.1, the calculations are as follows:
y1p = y0 + hf(x0, y0) = 0 + 0.1 × (1 – 0) = 0.1
h
y1c (1) = y0 + f ( x0 , y0 ) + f ( x0 , y0p )
2
0.1
= 0+ [(1 − 0) + (1 − 0.1)] = 0.0950.
2
182 Self-Instructional Material
We will repeat this step untill y1 becomes stationary. So, with x0 = 0, Numerical Solution to
Ordinary Differential
y0p = 0.0950 Equations
0.1
y1c (2) = 0 + [1 + (1 − 0.0950)] = 0.0953
2 NOTES
0.1
y1c (3) = 0 + [1 + (1 − 0.0953)] = 0.0952
2
0.1
y1c (4) = 0 + [1 + (1 − 0.0952)] = 0.0952
2
Since y1c (3) and y1c (4) are same, we will stop here for x1 = x0 + h = 0 +
0.1 = 0.1 and y1 = 0.0952 (i.e. at x = 0.1)
Similarly, we will calculate for x = 0.2 and x = 0.3 also. The calculations are
arranged as shown in Table 5.1.
Table 5.1 Calculations
at at x = 0.3 y = 2.0994
and x = 0.4 y = 2.1378 NOTES
dy
Let = f(x, y) ; y(x0) = y0 be an ordinary differential equation.
dx
To compute the value of y we use the two formulae predictor and corrector as,
(i) Milne’s predictor formula
4h
yn + 1 = yn – 3 +
3
( 2 f n − 2 − f n − 1 + 2 f n ) ; n = 3, 4, ...
(ii) Milne’s corrector formula
h
( f n − 1 + 4 f n + f n + 1 ); n = 1, 2, ...
yn + 1 = yn – 1 +
3
By knowing four consecutive values of y as yn – 3, yn – 2, yn – 1 and yn we
calculate yn + 1 using predictor method and then use yn + 1 in the corrector formula
to find the value of yn + 1.
Example 5.15: Apply Milne’s method, to find a solution of the differential equation
y′ = x2 – y in the range 0 ≤ x ≤ 1 for the boundary condition y = 0 at x = 0.
Solution. Using Picard’s method, we have,
x
y = y(0) + ∫0 f ( x, y ) dx, given that f(x, y) = x2 – y
x x 2 x3
so, y 1 = y0 + ∫0 ( x 2 − y ) dx = 0 + ∫0
x dx =
3
x x x3 x3 x 4
y 2 = y0 + ∫ f ( x, y1 ) dx = 0 + ∫ x 2 − dx = −
0 0 3 3 12
x x 2 x3 x 4
y 3 = y0 + ∫ f ( x, y2 ) dx = 0 + ∫ x − + dx
0 0 3 12
x3 x 4 x5
= − +
3 12 60
Now we will determine the starting values of y for Milne’s method by taking
h = 0.2.
4 × 0.2
= 0.0025 + [ 2 × 0.1406 – 0.2975 + 0.4987 × 2]
3
= 0.2641
Now applying corrector formula,
h
y5c = y3 + ( f3 + 4 f 4 + f5 )
3
Here, x5 = 1.0, y5 = 0.2641 and f5 = 0.7359
0.2
so, y5c (1) = 0.0625 + [0.2975 + 4 × 0.4987 + 0.7359]
3
= 0.2644
Again, with x5 = 1.0, y5 = 0.2644 and f5 = 0.7356
0.2
y5c (2) = 0.0625 + [0.2975 + 4 × 0.4987 + 0.7356]
3
= 0.2644.
h
and y4c = y2 + ( f 2 + 4 f3 + f 4 )
3
4 x3 7 x 4 6 x 5
2
so, y(x) ~ 1 + x + x + + +
3 6 5
∴ y(–0.1) = 0.9088, y(0.1) = 1.1115 and y(0.2) = 1.2529.
Thus,
x 0 = – 0.1 y0 = 0.9088 f(0) = 0.8359
x1 = 0 y1 = 1 f1 = 1
x 2 = 0.1 y2 = 1.1115 f2 = 1.2454
and x 3 = 0.2 y3 = 1.2529 f3 = 1.6098
Using Milne’s predictor formula,
4h
y4p = y0 + [2 f1 − f 2 + 2 f3 ]
3
4 × 0.1
= 0.9088 + [2 × 1 – 1.2454 + 2 × 1.6098]
3
= 1.4387
Now using Milne’s corrector method,
h
y4c = y2 + [ f 2 + 4 f3 + f 4 ]
3
with x 4 = 0.3, y4 = 1.4387 and f4 = 2.1598
0.1
y4c (1) = 1.1115 + × [1.2454 + 4 × 1.6098 + 2.1598]
3
= 1.4396 and hence f4(1) = 2.1624
again,
0.1
y4c (2) = 1.1115 + [1.2454 + 4 × 1.6098 + 2.1624]
3
= 1.4397 with f4(2) = 2.1628
5.7 SUMMARY
h
ync + 1 = yn + f n − 2 − 5 f n − 1 + 19 f n + 9 f np+ 1 .
24
• The numerical solutions of differential equations certainly differ from their
exact solutions. The total error at any stage say ‘r’ is the difference between
the computed value yr and the true value y(xr).
Total Error = Truncation Error + Round-off Error
• The main and the most important objective of numerical methods is to
minimize the error and to obtain most appropriate solutions. Truncation
error can be reduced in any method by taking small subintervals whereas
the round-off error cannot be controlled easily unless the computer has the
double precision arithmetic facility.
dy
1. The first order differential equation is = f (x,y) given y(xo) = yo.
dx
2. The first approximation by Picard’s method is,
x
x2
x
y x =1+ ∫
(1) ( x + 1) dx = 1 + x +
2
0
Short-Answer Questions
1. Write equation used in Picard’s method of successive approximations for
differentiate equation of first order?
2. To which type of function Taylor Series method can be applied?
3. Define Euler’s method.
4. Write equation used in Runge-Kutta second order method.
5. Write formulae the Milne’s Predictor and corrector methods.
Long-Answer Questions
dy
1. Using Taylor’s Series method, solve = x2 – y, y(0) = 1 at x = 0.1, 0.2,
dx
0.3 and 0.4. Compare the values with the exact solution.
2. Use Picard’s method to approximate the value of y when x = 0.1 given that
dy y − x
y = 1 when x = 0 and = .
dx y + x
3. Use Picard’s method to find two successive approximate solution of the
initial value problem,
dy y − x
= , y (0) = 1
dx y + x
d2y dy
12. Solve 2
+ x. + y = 0 , y(0) = 1 and y′(0) = 0 and obtain y for x = 0.1,
dx dx
0.2 and 0.3 by any method. Further apply Milne’s method to compute
y(0.4).
Mott, J.L. Discrete Mathematics for Computer Scientists, 2nd edition. New
Delhi: Prentice-Hall of India Pvt. Ltd., 2007.
UNIT 64PROBABILITY
DISTRIBUTION
NOTES
Structure
6.0 Introduction
6.1 Unit Objectives
6.2 Classical Approch to Probability
6.2.1 Sample Space
6.3 Random Variables
6.3.1 Types of Random Variables
6.3.2 Joint and Marginal Probability Density Function
6.4 Discrete Theoretical Distributions
6.4.1 Binomial Distribution
6.4.2 Moments
6.4.3 Poisson Distribution with Mean and Variance
6.4.4 Uniform Distribution
6.5 Summary
6.6 Key Terms
6.7 Answers to ‘Check Your Progress’
6.8 Questions and Exercises
6.9 Further Reading
6.0 INTRODUCTION
In this unit, you will learn about probability distribution. Probability distribution
refers to mathematical models of relative frequencies of a finite number of
observations of a variable by listing all the likely results of an experiment. It is a
symmetric arrangement of probabilities associated with mutually exclusive and
collectively exhaustive elementary events in an experiment.
This unit will also walk you through random variables and discrete theoretical
distributions. These are supported by a number of examples to enhance
understanding of the subject clearly.
The classical theory of probability is the theory based on the number of favourable
NOTES outcomes and the number of total outcomes. The probability is expressed as a
ratio of these two numbers. The term ‘favorable’ is not the subjective value given
to the outcomes, but is rather the classical terminology used to indicate that an
outcome belongs to a given event of interest.
Classical definition of probability: If the number of outcomes belonging
to an event E is NE , and the total number of outcomes is N, then the probability of
NE
event E is defined as pE =
N
For example, a standard pack of cards (without jokers) has 52 cards. If we
randomly draw a card from the pack, we can imagine about each card as a possible
outcome. Therefore, there are 52 total outcomes. Calculating all the outcome
events and their probabilities, we have the following possibilities:
• Out of the 52 cards, there are 13 clubs. Therefore, if the event of interest is
drawing a club, there are 13 favourable outcomes, and the probability of
13 1
this event becomes: = .
52 4
• There are 4 kings (one of each suit). The probability of drawing a king is:
4 1
= .
52 13
• What is the probability of drawing a king or a club? This example is slightly
more complicated. We cannot simply add together the number of outcomes
for each event separately (4 + 13 = 17) as this inadvertently counts one of
16
the outcomes twice (the king of clubs). The correct answer is: from
52
13 4 1
+ −
52 52 52
We have this from the probability equation, P(club) + P(king) – P(king of
clubs).
• Classical probability has limitations, because this definition of probability
implicitly defines all outcomes to be equiprobable and this can be only used
for conditions such as drawing cards, rolling dice, or pulling balls from urns.
We cannot calculate the probability where the outcomes are unequal
probabilities.
It is not that the classical theory of probability is not useful because of the
above described limitations. We can use this as an important guiding factor to
nA
The sequence in the limit that will converge to the same result every time, or
n
that it will converge at all. To understand this, let us consider an experiment consisting
of flipping a coin an infinite number of times. We want that the probability of heads
must come up. The result may appear as the following sequence:
HTHHTTHHHHTTTTHHHHHHHHTTTTTTTTHHHHHHHHHHHHH
HHHTTTTTTTTTTTTTTTT...
This shows that each run of k heads and k tails are being followed by
nA
another run of the same probability. For this example, the sequence oscillates
n
1 2
between, and which does not converge. These sequences may be unlikely,,
3 3
and can be right. The definition given above does not express convergence in the
required way, but it shows some kind of convergence in probability. The problem
of formulating exactly can be considered using axiomatic probability theory.
Empirical Probability Theory
The empirical approach to determine probabilities relies on data from actual
experiments to determine approximate probabilities instead of the assumption of
equal likeliness. Probabilities in these experiments are defined as the ratio of the
frequency of the possibility of an event, f(E), to the number of trials in the experiment,
n, written symbolically as P(E) = f(E)/n. For example, while flipping a coin, the
empirical probability of heads is the number of heads divided by the total number
of flips.
Self-Instructional Material 197
Probability Distribution The relationship between these empirical probabilities and the theoretical
probabilities is suggested by the Law of Large Numbers. The law states that as
the number of trials of an experiment increases, the empirical probability approaches
the theoretical probability. Hence, if we roll a die a number of times, each number
NOTES would come up approximately 1/6 of the time. The study of empirical probabilities
is known as statistics.
6.2.1 Sample Space
A sample space is the collection of all possible events or outcomes of an experiment.
For example, there are two possible outcomes of a toss of a fair coin: a head and
a tail. Then, the sample space for this experiment denoted by S would be:
S = [H, T]
So that the probability of the sample space equals 1, or
P[S] = P[H,T] =1
This is so because in the toss of the coin, either a head or a tail, must occur.
Similarly, when we roll a die, any of the six faces can come as a result of the roll
since there are a total of six faces. Hence, the sample space is S = [1, 2, 3, 4, 5,
6], and P[S] = 1, since one of the six faces must occur.
= E(X) + E( Y)
Thus, E(X + Y ) = E(X) + E( Y)
Corollary:E(X + Y + Z ... + ) = E(X) + E( Y) + E( Z) + ...
Notes: 1. Expected value of a constant is the constant itself.
i.e. E(C) = C where C is a constant
2. E(X) is also known as the mean of variable X.
3. If X is a random variable and f(X) is a function of X which takes the
values f(x1), f(x2), ... when X takes the values x1, x2, ... with probabilities P1, P2,
... respectively, then,
E(f(X)) = P1. f(x1) + P2. f(x2) + ... + Pn. f(xn)
and ΣPi = 1
Now if, f(X) = xr , then,
E(xr) = P1x1r + P2 x2r + ... + Pn x nr ... = ΣPi(x ri)
This is defined as the rth moment of the discrete probability distribution
about x = 0 and is denoted by µr′.
Also, µr = E(X – E(X))r
4. µ′ = E(X) = P1x1 + P2 x2 + ... + Pn x n = ΣPixi clearly.
5. The first moment about mean is zero. Let X be a random variable with
the following probability distribution,
X: x1 x2 ... xn
Y: p1 p2 ... pn
n
Then, E(X) = ΣXP(X) = ∑ xi Pi
i =1
where, σX = var( X ) and σY = var(Y ) are the standard deviations of variables NOTES
X and Y respectively.
Probability density function of continuous random variable: Consider a small
interval (x, x + dx) of length dx round the point x. Let f(x) be any continuous
function of x so that f(x) dx represents the probability that X falls in the infinitesimal
interval (x, x + dx). Then f(x) is called the probability density function of X. The
continuous curve y = f(x) is called the probability density curve.
Symbolically,
P(x ≤ X ≤ x + dx) = fX(x) or f(x).
The probability density function PDF of X, i.e. f(x) or fX(x) has the following
properties —
(i) f(x) ≥ 0 for all x
∞
(ii) ∫ f ( x)dx = 1
−∞
Distribution function: Let X be a random variable with PDF f(x), then the
function
∞
FX(x) = P(X ≤ x) = ∫ f (t )dt , – ∞ < x < ∞.
−∞
∫ x r f ( x)dx
r
E(x ) =
−∞
which is defined as µ′r , the rth moment (about origin) of the probability
distribution and if g(X) = [X – E(X)]r = [X – x ]r
∞
∫ (X − x)
r r
then, E(g(x)) = E[X – E(X)] = f ( x)dx
−∞
Which is denoted by µr, the rth moment about mean x .
Moment generating function: Let X be a random variable. Then moment
generating function M(t) is defined as
M(t) = E(etx)
t2 X 2
= E(1 + tX + + ...)
2!
t2 tr
= 1 + tE ( X ) + E ( X 2 ) + ... + E ( X r ) + ...
2! r!
t2 tr
= 1 + t. µ′1 + µ′2 + ... + µ′r
2! r!
tr
the coefficient of in the expression is the rth moment of X about origin (i.e.
r!
µ′r ).
This shows that moments are generated by moment generating function.
For a discrete random variable X which takes the values x1, x2, ... with
probabilities P1, P2, ... the moment generating function is defined as,
and for a continuous random variable X with probability density function f(x),
NOTES
–∞<x<∞
∞
∫e
tx
M(t) = E(etx) = . f ( x)dx .
−∞
t(x – x )
t 2 ( X − x )2
M x (t) = E(e ) = E 1 + t ( X − x ) + + ...
2!
t2 tr
= 1 + t. µ′1 + µ′2 + ... + µ′r + ...
2! r!
where µr is the rth moment about mean. It is also called the central
moment.
Characteristic function: Let X be a random variable. The function φ(t) is defined
as,
φ(t) = E(eitx)
is called the characteristic function of X.
Here, i is the imaginary quantity defined as −1 = i.
For a discrete random variable X with values x1, x2, ... and probabilities
P1, P2, ..., the characteristic function is defined as,
itx itx
φ(t) = E(eitx) = P1e 1 + P2e 2 + ...
= ∑ Pr eitxr
r
For a continuous random variable with PDF f(x) the characteristic function
φ(x) is defined as,
∞
φ(t) = E(e ) = itx
∫ f ( x). eitx dx
−∞
∫ ∫ f (t − r )e
i (t − r ) x
(ii) φ(x, k) = dtdr.
00
is real and non-negative for all real x and all k > 0.
Example 6.1: What is the expected value of the number of points that will be
obtained in a signal throw with an ordinary die? Find the variance also.
Solution: The random variable X assumes the values 1, 2, 3, 4, 5 and 6 with
1
probability P(X = xi) = in each case.
6
1 1 1 1 1 1
so, E(X) = ∑ xi Pi = 1. + 2. + 3. + 4. + 5. + 6.
i
6 6 6 6 6 6
= 3.5
1 2
and Var (X) = E(X2) – [E(X)]2 = (1 + 22 + ... + 62 ) – (3.5)2
6
35
=
12
Example 6.2: A random variable X has the following probability function :
X: 0 1 2 3 4 5 6
p( X ) : 0 k 2k 2k 3k k2 2k 2
⇒ 0 + k + 2k + 2k + 3k + k2 + 2k2 = 1
⇒ 8k + 3k2 = 1
− 8 ± 82 + 4 × 3 − 4 ± 19
or, 3k2 + 8k – 1 = 0 ⇒ k= =
6 3
(ii) P(X < 6) = P(X = 0) + P(X = 1) + ... + P(X = 5)
= 0 + k + 2k + 2k + 3k + k2 = 8k + k2
= 8k + 3k2 – 2k2 = 1 – 2k2.
P(X ≥ 6) = 1 – P (X < 6) = 1 – (1 – 2k2) = 2k2
P(0 < X < 4) = P(X = 1) + .... + P(X = 3) = k + 2k + 2k = 5k
(iii) Distribution of X.
X f X ( x) + P ( X ≤ x)
0 0
1 k
2 k + 2k = 3k
3 k + 2 k + 2 k = 5k
4 k + 2k + 2k + 3k = 8k
5 k + 2k + 2k + 3k + k 2 = 8k + k 2
6 k + 2k + 2k + 3k + k 2 + 2k 2 = 8k + 3k 2 = 1
x=k: 0 1 2 3
p( x = k ) = p(k ) : 1/ 8 3/8 3/8 1/ 8
3
1 3 3 1 13
So, E(X) = ∑Px
k =0
k k = 0. + 1. + 2. + 3. =
8 8 8 8 8
let X′ be the number of tails., then.
3
C x′
P(X′ = x′) = , x′ = 0, 1, 2, 3.
8
13
and thus, E(X) = E(X′) = .
8
Example 6.5: A continuous random variable X has a PDF
f (x) = 3x2 0 ≤ x ≤ 1, find a and b such that
(i) P{X ≤ a} = P{X > a}, and
(ii) P{X > b} = 0.05
1
Solution: (i) Since P(X ≤ a) = P(X > a) each must be equal to , because the
2
total probability is always 1.
1 1
∴ P(X ≤ a) = and P(X > a) =
2 2
a a
1 1 1
⇒ ∫ f ( x)dx = = ∫ 3 x dx =
2
P(X ≤ a) =
2 0
2 0
2
a
x3 1 1 1
⇒ 3 = ⇒ a3 – 0 = ⇒ a = 3
3 0 2 2 2
1
(ii) p(X > b) = 0.05 ⇒ ∫ f ( x)dx = 0.05
b
x 3 1
⇒ 3 = 0.05 ⇒ 1 – b3 = 0.05
3 b
1/ 3
19 19
b = 3
⇒b= .
20 20
NOTES Usually we use the notation X ~ B(n, P) to denote that the random variable X
follows binomial distribution with parameters n and P.
Notes: 1. Since n and P are independent constants in the binomial distribution ,
they are called the parameters of distribution.
n n
2. ∑ P (k ) = ∑ n Ck .Pk .Qn−k = (Q + P) n
=1
k =0 k =0
n n
= ∑ x. Cx .P n x
.Q n− x
+ ∑ x( x − 1).n Cx .P x .Qn− x
x =0 x =0
n
= nP + n(n – 1)P . ∑
n−2
2 Cx −2 .P x −2 . Q n− x
212 Self-Instructional Material
x=2
= nP + n(n – 1)P2.(Q – P)n–2 Probability Distribution
µ 23 n 2 P 2Q 2 (Q − P)2 (1 − 2 P)2
Hence, β1 = = =
µ32 n3 P3Q3 nPQ
Q−P 1 − 6PQ
γ1 = β1 = and γ2 = β2 – 3 =
nPQ nPQ
= ∑
x=0
n
Cx ( Pet ) x .Q n – x
M0(t) = (Q + Pet)n
n! = 2π. n n+1/ 2 e− n
n! n n+1/ 2 .e− n
so, =
(n − x )! (n − x) n− x +1/ 2 .e− n+ x .n x
Thus, Equation (6.1) becomes,
λ x −λ n n+1/ 2 .e − n
P(x) = .e .L im
x! n→∞ ( n − x ) n − x +1/ 2 .e − n + x .n x
λ x .e −λ n n+1/ 2 .e− n
= .L im
x !.e x n →∞ (n − x) n− x+1/ 2 .e − n .n x
λ x .e−λ n n − x+1/ 2
= .L im 1
x !.e x n→∞ n− x+
n − x +1/ 2 x 2
n 1 −
n
λ x .e−λ 1
= x
.L im 1
x !.e n→∞
x 2
n − x+
x
1 − . 1 −
n n
λ x .e−λ 1 λ x .e −λ
= . =
x !.e x e− x .1 x!
Thus,
λ x .e−λ
P(X = x) = P(x) = ; x = 0, 1, 2, ...
x!
When n → ∞, nP = λ and P → 0
Here, λ is known as the parameter of Poisson distribution.
Definition: A random variable X is said to follow a Poisson distribution if
it assumes only non-negative values and its probability mass function is
given by
λ x .e−λ
P(x, λ) = P(X = x) = ; x = 0, 1, 2, ...; λ > 0.
x!
λx
Let P(X = x) = e −λ. ; x = 0, 1, 2, ... ∞; λ > 0 be a Poisson distribution. Then
x!
the moment generating function is given by, NOTES
∞ −λ x
e .λ
M(t) = E(etx) = ∑ etx . x!
x =0
∞ t x
(λ. e )
= e −λ . ∑
x =0 x!
t t
= e−λ . eλe = eλ (e − 1)
− r Cr . p r (−q) x ; x = 0, 1, 2, ...
so p(x) =
0 ; otherwise
which is the (x + 1)th term in the expansion of Pr(1 – Q)–r, a binomial expansion
with a negative index. Hence, the distribution is called a negative binomial
distribution. Also,
∞ ∞
∑ P( x) = Pr ∑ −r
C x ( −Q ) x = Pr × (1 – Q)–r = 1
x =0 x =0
Therefore, P(x) represents the probability function and the discrete variable
which follows this probability function is called the negative binomial variable.
Example 6.6: A continuous random variable X has a probability distribution function
f(x) = 3x2, 0 ≤ x ≤ 1.
Find a and b such that
x 3 1
1
⇒ 3 = 0.05 ⇒ 1 – b3 =
3 b 20
1/ 3
b =
19
20
6.5 SUMMARY
Short-Answer Questions
1. What is classical definition of probability?
2. What is discrete random variable?
3. Define the binomial distribution.
4. Write the types of random variables.2
Long-Answer Questions
1. Prove that the expectation of the sum of the two random variables, is equal
to the sum of their expectation, i.e.
E(X + Y) = E(X) + E(Y)
where X and Y are two random variables.
2. How would you define the correlation coefficient ρXY or r(X,Y) for two
random variables X and Y?
3. Write a short note on probability density function of continuous random
variable.
4. List the various points about Moment Generating Function.
5. What is the expected value of the number of points that will be obtained in
a signal throw with an ordinary die? Find the variance also.
6. Describe the types of random variables along with their probabilities and
mathematical expectation.
7. Discuss joint and marginal probability density function.
8. Obtain the first four moments about origin of binomial distribution.
9. Describe uniform distribution and moments of uniform distribution.
Sancheti, D.C. and V.K. Kapoor. Business Mathematics. New Delhi: Sultan
Chand & Sons. NOTES
Levin, Richard I. 1991. Statistics for Management. New Jersey: Prentice-Hall.
Chance, William A. 1969. Statistical Methods for Decision Making. Illinois:
Richard D Irwin.
Meyer, Paul L. 1970. Introductory Probability and Statistical Applications.
Massachusetts: Addison-Wesley Publishing Co.
Johnson, R.A., Probability and Statistics for Engineers, New Delhi: PHI
Trivedi, K.S., Probability and Statistics with Reliability, New Delhi: PHI, 1994
Levin, Richard. I., and David. S. Rubin. 1997. Statistics for Management, 7th
edition, New Jersey: Prentice Hall International.
Gupta, S.C., and V.K. Kapoor. 1997. Fundamentals of Mathematical Statistics.
9th Revised Edition, New Delhi: Sultan Chand & Sons.
Freud, J.E., and F.J. William. 1997. Elementary Business Statistics – The
Modern Approach. 3rd edition, New Jersey: Prentice Hall International.
Goon, A.M., M.K. Gupta, and B. Das Gupta. 1983. Fundamentals of Statistics.
Vols. I & II, Calcutta: The World Press Pvt. Ltd.
Hogg, Robert. V., Allen T. Craig and Joseph W. McKean. 2005. Introduction to
Mathematical Statistics. New Delhi: Pearson Education.
7.0 INTRODUCTION
Approximation theory deals with two types of problems, (i) when a function is
given explicitly and you want to find a simpler type such as polynomial for
representation and (ii) the problem concerns fitting function for the given data and
finding the best function in a certain class that is used to represent the data.
The major aim of approximation of functions is to represent a function with
minimum error as this is a central problem in software development. As you know,
computers are essentially arithmetic devices, the most elaborate function they can
compute is rational function, a ratio of polynomials. There are two ways to
approximate a function using a polynomial:
(1) Using truncated Taylor’s series
(2) Using Chebyshev polynomials
In this unit, you will learn to approximate a function by using these two
techniques.
In addition, you will learn about the normal distribution of variables and
polynomials.
Let f(x) be a function which has derivatives upto (n + 1)th in an interval [a, b].
This function may be expressed near x = x0 in [a, b] as follows:
( x − x0 ) 2
f(x) = f ( x0 ) + ( x − x0 ) f ′( x0 ) + f ′′( x0 ) + ...
2!
( x − x0 ) n n ( x − x0 ) n + 1 ( n + 1)
+ ... + f ( x0 ) + ·f (k )
n! ( n + 1)!
( x − x0 ) n + 1 f n + 1 ( k )
The term is called the remainder term. Here, k is a
( n + 1)!
function of x and lies between x and x0, i.e. x ≤ k ≤ x0.
This remainder term gives the truncation error only if the first n-terms in
Taylor series are used to represent the function.
( x − x0 ) n + 1
Hence,the truncation error (T.E.) = f ((kn) + 1) .
( n + 1)!
| ( x − x0 )n + 1 |
≤ M
(n + 1)!
( x − x0 ) 2 NOTES
f(x) = f ( x0 ) + f ′( x0 ) ( x − x0 ) + f ′′( x0 ) + ...
2!
x2 x3
= sin (0) + cos (0) ( x ) – sin (0). – cos (0) + ...
2! 3!
Thus,
x3 x 5 x 7
sin x = x – + − + ...
3! 5! 7!
x3 x5 x7
or, sin x = x − + − + ... ...(7.1)
6 120 5040
Since sin x is required correctly upto three significant digits and the truncation
error after three terms of equation is as follows:
1
T.E. ≤ = 0.000198
5040
Thus, the series is truncated after three terms.
x3 x5
So, sin x = x − +
3! 5!
This is the required representation to obtain the value of sin x correctly upto
three significant digits.
The Chebyshev polynomials {Tn(x)} are orthogonal on (–1, 1) with respect to the
weight function w(x) = (1 – x2)–1/2. The Chebyshev polynomial is defined by the
following relation:
For x ∈ [–1, 1], define
Tn(x) = cos (n cos–1 (x)) for each n ≥ 0.
It is not obvious from this definition that Tn(x) is an nth degree polynomial in
x, but now it will be proved that it is. First note that
T0(x) = cos 0 = 1 and T1(x) = cos (cos–1(x)) = x.
For n ≥ 1, introduce the substitution θ = cos–1 x to change this
equation to
Returning to the variable x and solving for Tn + 1 ( x) you have, for each n ≥ 1,
Tn + 1 (θ) = 2cos(n arccos x). x − Tn − 1 ( x) = 2Tn ( x) x − Tn − 1 ( x).
x
–1 1
T2(x)
–1
x = T1 ( x)
1
x2 = [T2 ( x ) + T0 ( x )]
2
1
x3 = [3T1 ( x ) + T3 ( x )]
4
1
x4 = [3T0 ( x) + 4T2 ( x ) + T4 ( x )]
8
1
x5 = [10T1 ( x ) + 5T3 ( x ) + T5 ( x )]
16
1
and x6 = [10T0 ( x ) + 15T2 ( x ) + 6T4 ( x ) + T6 ( x )]
32
These expressions are useful in the economization of power series.
Thus, after omitting T5(x) you have, sin x – 0.8802 T1 ( x) − 0.03906 T3 ( x)
substituting T1 ( x) = x and T3 ( x) = 4 x3 − 3 x you have,
sin x = 0.9974x – 0.1562x3
which gives sin x correctly upto three significant digits with only two terms for any
value of x.
7.3.1 Chebychev Inequality
If X is a random variable having µ as expected value and σ2 as finite variance,
then for a positive real number k
Pr(X = – 1) = 1/(2k2),
Pr(X = 0) = 1 – 2k2,
Pr(X = 1) = 1/(2k2),
Pr(|X – µ| ≥ kσ) = 1/k2,
Equality holds exactly in case of a linear transformation whereas inequality
holds for distribution which is non-linear transformation.
Theorem is useful as it applies to random variables for any distribution and
these bounds can be computed by knowing only mean and variance.
Any observation, howsoever accurate it may be, is never more than few
standard deviations away from the mean. Chebyshev’s inequality gives following
bounds that apply to all distributions in which it is possible to define standard
deviation.
At least:
• 50% of these values lie within standard deviations = Ö2
• 75% of the values lie within standard deviations = 2
• 89% lie within standard deviations = 3
• 94% lie within standard deviations = 4
• 96% lie within standard deviations = 5
• 97% lie within standard deviations = 6
Standard deviations are always taken from the mean.
Generally:
Minimum (1 – 1/k2) × 100% lie with standard deviations = k.
x3 x5
Example 7.2: Represent sin x = x − + ... by using Chebyshev polynomial
3! 5!
to obtain three significant digits accuracy in the computation of sin x.
Solution: Using Chebyshev polynomial, sin x is approximated as:
−x x 2 x3 x 4
e = 1− x + − + + ....
2! 3! 4!
in terms of Chebyshev polynomials.
Solution: The Chebyshev polynomial representation of
x 2 x3
e− x = 1 − x + − + ....
2! 3!
−x 1 1
as e = T0 ( x ) − T1 ( x ) + [T2 ( x ) + T0 ( x)] − [3T1 ( x) + T3 ( x )]
4 24
1 1
+ [3T0 ( x ) + 4T2 ( x ) + T4 ( x )] − [10T1 ( x ) + 5T3 ( x ) + T5 ( x)] + ...
195 1920
Thus,
e− x = 1.26606 T0 ( x) − 1.13021 T1 ( x) + 0.27148 T2 ( x)
7.4 DISTRIBUTION
Binomial, Poisson, negative binomial and uniform distribution are some of the
discrete probability distributions. The random variables in these distributions assume
a finite or enumerably infinite number of values but in nature these are random
variables which take infinite number of values i.e. these variables can take any
value in an interval. Such variables and their probability distributions are known as
continuous probability distributions.
A random variable X is the said to be normally distributed if it has the following
probability density function:
2
1 x – µ
1 –
2 σ
f(x) = e , for – ∞ ≤ x ≤ ∞
σ 2π
where µ and σ > 0 are the parameters of distribution.
Normal Curve: A curve given by
2
1 x – µ
–
2 σ
yx = y0e
which is known as the normal curve when origin is taken at mean.
1 x2
–
2 σ2
then yx = y0e .
y0 yx
O
Standard Normal Variate : A normal variate with mean zero and standard
deviation unity, is called a standard normal variate.
That is; if X is a standard normal variate then E(X) = 0 and V(X) = 1.
232 Self-Instructional Material
Then, X ~ N (0, 1) Approximation Theory
X – µ 5−8
Solution: (i) P[X ≤ 5] = P ≤
σ 4
= P (Z ≤ – 0.75)
= P (Z ≥ 0.75)
0.5
[By Symmetry]
= 0.5 – P(0 ≤ Z ≤ 0.75)
[To use relevant table]
= 0.5 – 0.2734 [See Appendix for value of ‘2’]
= 0.2266.
5.8 10 − 8
(ii) P[5 ≤ X ≤ 10] = P ≤Z ≤
4 4
X = 40
Z = – 0.8 Z = 0 Z=2
X –µ
(i) When X = 26, Z = = – 0.8
σ
X –µ
and for X = 40, Z = =2
σ
∴ P[26 ≤ X ≤ 40] = P[– 0.8 ≤ Z ≤ 2]
= P[0 ≤ Z ≤ 0.8] + P[0 ≤ Z ≤ 2]
= 0.2881 + 0.4772 = 0.7653
(ii) P[| X – 3 | > 5] = 1 – P[| X – 3 | ≤ 5]
P[| X – 3 | ≤ 5] = P[25 ≤ X ≤ 35]
25 − 30 35 − 30
= P ≤Z≤
5 5
= 2.P (0 ≤ Z ≤ 1) = 0.
= 2 × 0.3413 = 0.6826.
so P[| X – 3| > 5] = 1 – P[| X – 3 | ≤ 5]
= 1 – 0.6826 = 0.3174.
7.4.1 Central Limit Theorem
Let X 1, X 2,.... X n be n independent random variables all of which have the
same distribution. Let the common expectation and variance be µ and σ2
respectively.
–µt
n t
( X1 + X 2 + ... + X n )
σ
.E e σ n
=e
n
–µt t
=e
σ . M ( X1 + X 2 + ... + X n ) .
σ n
n
n
. M x
–µt t
σ
=e t
σ
σ n
This is because the random variables are independent and have the same
MGF by using logarithms, you have:
− µt n t
log Mz(t) = + n log M x
σ σ n
− µt n µ′1t µ′2 t i
2
=
+ n log 1 + + + ...
σ σ n 2! σ n
236 Self-Instructional Material
Approximation Theory
− µt n µ′ t µ′ t t 2 1 µ′1t
2
= + n 1
+ 2
. 2 + ... – + ... + ...
σ σ n 2! n σ 2 σ n
t2
= + O (n – 1/ 2 ) [Q µ2′ – µ1′2 = σ2µ1′ = µ]
2
Hence, as n → ∞
t2 2
log(Mz)(t) → i.e. M2(t) = et /2
2
However, this is the M.G.F. of a standard normal random variable. Thus,
the random variable Z converges to N.
This follows that the limiting distribution of x as normal with mean µ and
σ2
variance .
n
Before starting the laws of large numbers, let us define an inequality named as
Kolmogorov’s inequality. The set of Kolmogorov’s inequalities was defined by
Kolmogorov in 1928.
Suppose X1, X2, ...., Xn is a set of independent random variables having
mean O and variances σ12, σ22, .... σn2.
Let C n 2 = σ12 + σ22 + .... + σn2.
Then the probability that all of the inequalities
1
| x1 + x2 + .... + xα| < λCn, α = 1, 2, ...n hold, is at least 1 – 2 .
λ
7.5.1 Weak Law of Large Numbers
x1 + x2 + .... + xn x + x2 + .... + xn
If you put, x = ,x= 1
n n
En
and λ = = α.
n2
where ∝ is an arbitrary positive number and En is the mathematical expectation of
the variate
En En
≥1− = 1 − η provided < η.
n2 α2 n α2
2
This is known as the Weak Law of Large Numbers. This can also be stated
as:
With the probability approaching unity or certainty as near as you please,
you can expect that the arithmetic mean of the values actaully assumed by n variates
will differ from the mean by their expectations by less than any given number
however small, provided the number of variates can be taken sufficiently large and
En
provided the condition → 0 as → ∞ is fulfilled.
n2
In other words, a sequence X1, X2, ... Xn of random variables is said to
satisfy the weak law of large numbers if
Sn Sn
lim P − E < ε = 1
n→∞ n n
βn
for any ε > 0 where Sn = X1 + X2 + ... + Xn. The law holds provided → 0 as
n2
n → ∞ where Bn = Var. (Sn) < ∞.
7.5.2 Strong Law of Large Numbers
Consider the sequence X1, X2,....Xn of independent random variables with
expectation xk or µ k = E(Xk) and variance σ2.
If Sn = X1 + X2 + .... + Xn and E(Sn) = mn then it can be said that the sequence S1,
S2, .... Sn obeys the strong law of large numbers if every pair ε > 0, δ > 0
corresponds to N such that there is a probability (1 – δ) or better that for every r
> 0 all r + 1 inequalities,
Sn – mn
<ε
n
n = N, N + 1, ... N + r will be satisfied.
Bn n 1
∴ lim 2 = lim = lim = 0.
n→∞ n x→∞ 2 x→∞ n
n
Hence, weak law of large numbers holds for the sequence {Xk} of
independent random variables.
Example 7.8: Examine whether the laws of large numbers holds for the sequence
{Xk} independent random variables which are defined as follows:
1
P(Xk = ± k–1/2) =
2
Solution: E(Xk) = ΣXk pk
1 1
–
1 – 1
= k × + (−k 2 ) × = 0
2
2 2
– 1/ 2 2 1 1
E(Xk2) = Σxk2 . pk = ( k ) × + ( − k – 1/ 2 ) 2 ×
2 2
1 –1 1 –1 –1
= ×k + k =k
2 2
Var. (Xk) = E(Xk ) – [E(Xk)]2 = k–1 – 0 = k–1
2
n n
Bn = ∑ Var (X k ) = ∑ k–1 = n . k–1
i =1 i =1
Bn n . k–1 k–1
∴ lim = = lim =0
n → ∞ n2 n2 n→∞ n
9 − 10 X – 10 11 − 10
= P < <
5 5 5
T3 ( x) = 4x3 – 3x
T4 ( x ) = 8x4 – 8x2 + 1
NOTES
5 3
T5 ( x) = 16x – 20x + 5x
T6 ( x) = 32x6 – 48x4 + 18x2 – 1
4. Binomial, poisson, negative binomial and uniform distribution are some of
the discrete probability distributions.
5. A normal variate with mean zero and standard deviation unity, is called a
standard normal variate.
6. The set of Kolmogorov’s inequalities was defined by Kolmogorov in 1928.
1
7. When number of trials is large and probability p is close to , normal
2
approximation can be used to for binomial as well as poission distribution.
8.0 INTRODUCTION
In this unit you will learn about statistical inferences. Population in statistics does
not necessarily mean human population alone. It means a complete set of objects
under study, living or non-living. Each individual or object is called unit or member
or element of that population. If there are a finite number of elements then it is
called finite population. If the number of elements is infinite then it is called infinite
population. The procedure of testing the validity of our hypothesis using sample
statistics is called Testing of Hypothesis.
The techniques of regression analysis are employed in statistics for modeling
and analysing several variables, when the focus is on the relationship between a
dependent variable and one or more independent variables. More particularly,
regression analysis helps us understand how the typical value of the dependent
variable changes when any one of the independent variables is varied, while the
other independent variables are held fixed. Generally, regression analysis estimates
the conditional expectation of the dependent variable given the independent
A universe is the complete group of items about which knowledge is sought. The
universe may be finite or infinite. Finite universe is one which has a definite and
certain number of items but when the number of items is uncertain and infinite, the
universe is said to be an infinite universe. Similarly the universe may be hypothetical
or existent. In the former case the universe in fact does not exist and we can only
imagine the items constituting it. Tossing of a coin or throwing of a dice are examples
of hypothetical universes. Existent universe is a universe of concrete objects, i.e.,
the universe where the items constituting it really exist. On the other hand, the term
sample refers to that part of the universe which is selected for the purpose of
investigation. The theory of sampling studies the relationships that exist between
the universe and the sample or samples drawn from it.
8.2.1 The Two Concepts: Parameter and Statistic
It would be appropriate to explain the meaning of two terms viz., parameter and
statistic. All the statistical measures based on all items of the universe are termed
as parameters whereas statistical measures worked out on the basis of sample
studies are termed as sample statistics. Thus, a sample mean or a sample standard
deviation is an example of statistic whereas the universe mean or universe standard
deviation is an example of a parameter.
246 Self-Instructional Material
The main problem of sampling theory is the problem of relationship between Statistical Inferences
In statistics, population does not mean human population alone. A complete set of
objects under study, living or non-living is called population or universe for example,
graduates in Pondy, Bajaj tube lights, Ceat car tyres etc. Each individual or object
is called a unit or member or element of that population. If there are a finite number
of elements then it is called finite population. If the number of elements is infinite
then it is called infinite population. Take certain examples to know the parameters
of a population like, average life span of a Bajaj bulb, coverage mileage of ceat
car tyre, percentage of defective products of a company etc. Suppose you want
to know the average mileage of Ceat car tyre. All the tyres manufactured by that
company will constitute the population. Each tyre is run and the average mileage is
calculated which is known as the population average. This average calculated for
the entire population is called a parameter. Similarly any constant calculated for
the whole population is called a parameter.
Suppose there are N tyres and all of them are moving. Let X1, X2, ... Xn be
the mileages of N tyres. Now these N numbers constitute a population. These N
values are the values of a single variable X namely mileage. This variable X may
have a probability density function f(X, θ) where θ is an unknown parameter. This
Pdf f(X, θ) describes the population completely when θ is known. Now, this f(X,
θ) is itself considered as a population. Since it is not possible to test each and
every tyre, we do not get the data from that complete population and hence the
parameter cannot be calculated. When it is unknown, we may assume that the
form of f(X, θ) is known. When the parameter θ is unknown an assumption of
possible value of θ is called hypothesis. To test the validity of the hypothesis we
make use of sample observations and statistics. The procedure of testing the validity
of our hypothesis using sample statistics is called Testing of Hypothesis.
| x −µ|
For testing a hypothesis a value z, where is calculated from the
σ/ n
given data x is sample mean µ is population mean σ is standard deviation z = test
statistics, 1.281 is the z value at 10% significance level, n is sample popular
= n /5
250 Self-Instructional Material
For example, to ensure that 90 per cent of bulbs have more than 975 hours Statistical Inferences
of life
| x − µ|
≤ 1.281
σ/ n
NOTES
| 975 − 1000|
i.e., n = n /5 = 1.281
125
n = 6.405
= 6.4
n = 41 (approximately)
0.475 of 0.475 of
area area
When the population mean is either lower or higher than some hypothesised value,
one-tailed test is considered to be appropriate where the rejection is only on the
left tail of the curve. This is known as left-tailed test (refer Figure 8.2).
NOTES
For example, what will happen if the acceptance region is made larger? α will
decrease. It will be more easily possible to accept H0 when H0 is false (Type II
error), i.e., it will lower the probability by making a Type I error, but raise that of
β, Type II error.
Solution: No. Each is concerned with a different type of error. But both are not
independent of each other.
The level of significance implies the probability of type I error. A 5 per cent
level implies that the probability of committing a type I error is 0.05. A 1 per cent
level implies 0.01 probability of committing type I error.
Lowering the significance level and hence the probability of type I error is
good, but unfortunately it would lead to the undesirable situation of committing
type II error.
Type I Error: Reject H0 when it is true.
Type II Error: Accept H0 when it is wrong, i.e., reject H1 when it is true.
In practice, type I error amounts to rejecting a lot when it is good and
type II error may be regarded as accepting a lot when it is bad.
Thus P(rejecting a lot when it is good) = α
P(accepting a lot when it is bad) = β
where α and β are referred to as Producer’s risk and consumer’s risk respectively.
Power of Test: A good test should accept H0 when it is true and reject H0 when
it is false. Power of test or (1 – β) measures how well the test is working and is
called the power of test.
No. of sample pts. in the sample leading
to accept H 0 when it is false
Power of rest = (1 – β) =
Total no. of sample pts.
Degree of freedom: For a fixed value of the mean the number of free choices is
called the degree of freedom.
8.4.2 Testing a Hypothesis
Step 1: Set up H0. The following steps facilitates the testing of a hypothesis.
8.5 REGRESSION
n n n
and ∑ xi yi = a ∑ xi + b ∑ xi2 ...(8.3)
i =1 i =1 i =1
n n
1 1
also σX2 =
n
∑ xi2 − x 2 ⇒ n0
∑ xi2 = σ2x + x 2
i =1 i =1
Cov (x, y )
⇒ b=
σ2X
Solving Equations (8.2) and (8.3), we get
n Σxi yi − Σxi . Σyi
bˆ =
nΣxi2 − (Σxi ) 2
and ˆ .
â = y − bx
where â and bˆ are the least square estimates of (a, b) therefore the least squared
estimated value yˆi of yi is given by
yˆi = aˆ + bˆ xi .
Since b is the slope of the line of regression of Y on X and since the line of
regression passes through the point ( x , y ) , its equation is
Cov (X , y )
Y − y = b( X − x ) = (X − x)
σ2X
σY
⇒ Y − y = r. (X − x)
σX
r 2 − 1 σ X . σY
tan θ = .
r σ2X + σY2
π
(ii) If r = 0, then tan θ = ∞ and θ = , i.e., if the two variables are not
2
related then the lines of regression becomes perpendicular to each other.
(iii) If r = ±1 ⇒ tan θ = 0 ⇒ θ = 0°. That is, the lines of regression will
concide.
8.5.2 Regression Coefficients
Regression coefficients of Y on X and X on Y are given as:
σY nΣxy − (Σx) (Σy )
byx = r . =
σX nΣx 2 − (Σx)2
1 − r 2 σ X σY
tan θ = . as r2 ≤ 1
r σ2X + σY2
Example 8.2: The equation of two regression lines, obtained in a correlation
analysis of 60 observations are
5x = 6y + 24 and 1000 y = 768x – 3608
what is the correlation coefficient? Show that the ratio of coefficient of variability
of x to that of Y is 5/24. What is the ratio of variances of x and y?
Solution: 5x = 6y + 24 ⇒ x = 6/5y + 24/5
bxy = 6 / 5 .
σX
r = 6/5 ...(1)
σY
768 σY 768
Similarly, byx = ⇒ r = ...(2)
1000 σX 1000
By multiplying Equations (1) and (2), we get
⇒ r2 = 0.9216 ⇒ r = 0.96
divide Equation (2) by (1)
σ2X σX
= 1.5625 ⇒ = 1.25 = 5 / 4
σY2 σY
Since regression lines pass through the points ( x , y )
5x = 6 y + 24 and 1000 y = 768 x − 3608
On solving these two, we have x = 6, y = 1
σX
Coefficient of variability of X =
X
σ
Coefficient of variability of Y = Y
Y
σX Y Y σ 1 5 5
ratio = × = X = =
X σY X σY 6 4 24
x: 0 1.0 2 3 4
y: 1 4 10 17 30 .
x y x2 x3 x4 xy x2 y
0 1 0 0 0 0 0
1 4 1 1 1 4 4
2 10 4 8 16 20 40
3 17 9 27 81 51 153
4 30 16 64 256 120 480
Till now we have considered the case where dependent variable is a function of
single independent variable, i.e., y = f(x) type, it is a case of simple regression.
Here, we will consider a case where the dependent variable is a function of two or
more independent variables, i.e., y = f(x1, x2, ... xn) type. We call this type of
regression as multiple regression.
Although there are many different formulae to express this type of regression,
the most widely used are linear equation of the form
y = a + bx + cz
where a, b and c are regression coefficients which are to be determined. Applying
the method of least squares to determine these coefficients we minimize the sum of
square of the deviations.
That is,
e = Σ( yi − a − bxi − czi ) 2
Differentiating with respect to a, b and c we have,
∂e
= 2Σ( yi − a − bxi − czi ) (−1) = 0.
∂a
∂e
= 2Σ( yi − a − bxi − czi ) (− xi ) = 0.
∂b
∂e
and = 2Σ( yi − a − bxi − czi ) (− zi ) = 0.
∂c
On simplification, we get
or Σyi = n . a + bΣxi + cΣzi
2
Σxi yi = aΣxi + + bΣxi + cΣxi zi
and Σzi yi = a Σ zi + bΣxi zi + cΣ zi2
This is a system of three linear equations in three unknowns. On solving, we
will have the values of a, b and c and hence we determine the equation of regression
equation or plane.
2(−49)
1393 −
r= 11
2
2
492
1414 − . 1865 −
11 11
= 0.918
The regression coefficients of X on Y
ΣX . ΣY
ΣXY −
σ n
= r. X = 2
= 00.85
σY 2 (ΣY )
ΣY −
n
Σy 18
y = = =1
n 18
2 2
Σx 2 Σx Σy 2 Σy
σX = − σY = −
n n n n
2 2
60 12 96 18
= − = 2.88 = −
18 18 18 18
= 4.33
Self-Instructional Material 263
Statistical Inferences
so, σX = 2.88 = 1.7
= 0.56.
Regression equation of y on x:
σY
y − y = r. .( x − x )
σX
2.08
y – 1 = 0.56 × (x – 0.67)
1.7
1.7
x – 0.67 = 0.56 ( y − 1) or x = 0.46 y + 0.21
2.08
Example 8.6: Find a non-linear regression for the following data:
x: 1 2 3 4 5 6 7 8 9
y: 2 6 7 8 10 11 11 10 9
n x y x2 x3 x4 xy x2 y
1 1 2 1 1 1 2 2
2 2 6 4 8 16 12 24 NOTES
3 3 7 9 27 81 21 63
4 4 8 16 64 256 32 128
5 5 10 25 125 625 50 250
6 6 11 36 216 1296 66 396
7 7 11 49 343 2401 77 539
8 8 10 64 512 4096 80 640
9 9 9 81 729 6561 81 729
n=9 Σx = 45 Σy = 74 Σx 2 Σx 3 Σx 4 Σxy Σx 2 y
= 284 = 2025 = 15333 = 421 = 2771
on substituting these values in the normal equations, we have
a = –0.923, b = 3.520 and c = –0.267
Hence, the required non-linear regression is;
x: 1 2 3 4
y: 0 1 2 3
z: 12 18 24 30
Solution: Using y = a + bx + cz
equation will be 84 = 4a + 10b + 6c
240 = 10a + 30b + 20c and
156 = 6a + 20b + 14c
⇒ a = 10, b = 2 and c = 4.
So, regression plane will be
y = 10 + 2 x + 4 z. Ans.
X Y Z Z2 X2 XY YZ XZ
1 12 0 0 1 12 0 0
2 18 1 1 4 36 18 2
3 24 2 4 9 72 48 6
4 30 3 9 16 120 90 12
ΣX = 10 ΣY = 84 ΣZ = 6 14 30 240 156 20
266 Self-Instructional Material
8.5.5 Goodness of Fit Statistical Inferences
1.
(2) (3) (4) (5)
(1)
Five possible forms, which Scatter diagram may assume has been depicted in the above
five diagrams. First diagram is indicative of perfect positive relationship, Second shows
perfect negative relationship, Third shows no relationship, Fourth shows positive
relationship and Fifth shows negative relationship between the two variables under
consideration.
variable. Some formal expression of the relationship between the two variables is
necessary for predictive purposes. For the purpose, one may simply take a ruler and
draw a straight line through the points in the scatter diagram and this way can determine
NOTES
the intercept and the slope of the said line and then the line can be defined as Yˆ = a + bX i ,
with the help of which we can predict Y for a given value of X. But there are shortcomings
in this approach. For example, if five different persons draw such a straight line in the
same scatter diagram, it is possible that there may be five different estimates of a and b,
specially when the dots are more dispersed in the diagram. Hence, the estimates cannot
be worked out only through this approach. A more systematic and statistical method is
required to estimate the constants of the predictive equation. The least squares method
is used to draw the best fit line.
Y-axis
Consumption Expenditure (’ 00 Rs)
120
100
80
60
40
20
X-axis
0 20 40 60 80 100 120
Figure 8.4 Scatter Diagram, Regression Line and Short Vertical Lines Representing ‘e’
Self-Instructional Material 269
Statistical Inferences In the above figure, the vertical deviations of the individual points from the
line are shown as the short vertical lines joining the points to the least squares line.
These deviations will be denoted by the symbol ‘e’. The value of ‘e’ varies from
one point to another. In some cases it is positive, while in others it is negative. If the
NOTES line drawn happens to be least squares line, then the values of ∑ ei is the least
possible. It is so, because of this feature the method is known as Least Squares
Method.
Why we insist on minimizing the sum of squared deviations is a question that
needs explanation. If we denote the deviations from the actual value Y to the
n
estimated value Yˆ as (Y – Yˆ ) or ei, it is logical that we want the Σ(Y – Yˆ ) or ∑ ei , to
i =1
n
be as small as possible. However, mere examining Σ(Y – Yˆ ) or ∑ ei , is inappropriate,
i =1
since any ei can be positive or negative. Large positive values and large negative
values could cancel one another. But large values of ei regardless of their sign,
n
indicate a poor prediction. Even if we ignore the signs while working out | ei | ,
i =1
the difficulties may continue. Hence, the standard procedure is to eliminate the
effect of signs by squaring each observation. Squaring each term accomplishes
two purposes viz., (i) it magnifies (or penalizes) the larger errors, and (ii) it cancels
the effect of the positive and negative values (since a negative error when squared
becomes positive). The choice of minimizing the squared sum of errors rather than
the sum of the absolute values implies that there are many small errors rather than
a few large errors. Hence, in obtaining the regression line, we follow the approach
that the sum of the squared deviations be minimum and on this basis work out the
values of its constants viz., ‘a’ and ‘b’ also known as the intercept and the slope of
the line. This is done with the help of the following two normal equations:2
ΣY = na + bΣX
ΣXY = aΣX + bΣX2
In the above two equations, ‘a’ and ‘b’ are unknowns and all other values
viz., ∑X, ∑Y, ∑X2, ∑XY, are the sum of the products and cross products to be
calculated from the sample data, and ‘n’ means the number of observations in the
sample.
The following examples explain the Least squares method.
2. If we proceed centering each variable, i.e., setting its origin at its mean, then the two
equations will be as under:
∑Y = na + b∑X
∑XY = a∑X + b∑X2
But since ∑Y and ∑X will be zero, the first equation and the first term of the second
equation will disappear and we shall simply have the following equations:
∑XY = b∑X 2
b = ∑XY/∑X2
The value of ‘a’ can then be worked out as:_ _
a = Y – bX
270 Self-Instructional Material
Example 8.10: Fit a regression line Yˆ = a + bX i by the method of Least squares Statistical Inferences
When we talk about prediction or estimation, we usually imply that if the relationship
Yi = a + bXi + ei exists, then the regression equation, Yˆ = a + bX i provides a base
for making estimates of the value for Y which will be associated with particular
values of X. In Example 6.1, we worked out the regression equation for the income
and consumption data as,
Yˆ = 14.000 + 0.616Xi
On the basis of this equation we can make a point estimate of Y for any
given value of X. Suppose we wish to estimate the consumption expenditure of
individuals with income of Rs 10,000. We substitute X = 100 for the same in our
equation and get an estimate of consumption expenditure as follows:
Yˆ = 14.000 + 0.616 (100) = 75.60
ΣY 2 − aˆΣY − bˆΣXY
σ̂ = (n − 2)
For example, In a study of how wheat yield depends on the fertilizer input, ten
experimental plots were selected for study. The yields (Y) at different levels of
fertilizer input (X) were observed as given in table.
Fertilizer Yield (Y)
plot no. input (X) kg/acre X2 MT/acre XY Yˆ
(10).(18,175) − (565).(285)
= = 0.27
(10).(39,575) − (565) 2
ΣYi bˆΣX i
â = – = 13.245.
n n
Therefore, estimated line of regression is
Yˆi = 13.245 + 0.27Xi
Coefficient of determination R2 is given by
Σ(Y − Yˆ )2 11.717
R2 = 1 – 2
=1−
Σ(Y − Y ) 573.0
= 0.98 ~ 1
The standard error of estimate is given by
Σ(Y − Yˆ )2 11.717
σ̂ = = = 1.21
n−2 8
We can also plot the raw data Y and X and the estimated line on a graph.
8.6 SUMMARY
population.
• A very important aspect of sampling theory is the study of test of significance
which gives us a ground to decide the deviation between the observed
NOTES
sample statistic and the hypothetical parameter value or the deviation between
two independent sample statistics.
• A statistical hypothesis is a statement about a population parameter. There
are two types of statistical hypothesis which are null and alternative hypothesis.
• The Level of Confidence or level of significance expresses the degree of
credibility that may be attached to the result. The probability that we associate
with the interval estimate of a population parameter indicating how confident
we are that the interval estimate will include the true population parameter is
called confidence level. Higher probe means more confidence.
• Level of Significance is defined as the probability of rejecting the null
hypothesis when it is time. For example, 5% level of significance means
that there are about 5% changer or out of 100, you may reject the null
hypothesis 5 times, when it is time. The probability that a random value
of the test statistic t belongs to the critical region is known as the level
of significance.
• Regression analysis refers to the techniques used for modeling and analyzing
several variables, when the focus is on the relationship between a dependent
variable and one or more independent variables.
• If the variables in a bi-variate distribution are related, we will find that the
points in the scatter diagram will cluster round some curve called the ‘curve
of regression’.
• If the curve is a straight line it is called line of regression and there is said to
be linear regression between the variables, otherwise regression is said to
be curvilinear.
• Although there are many different formulae to express this type of regression,
the most widely used are linear equation of the form
y = a + bx + cz
where a, b and c are regression coefficients which are to be determined.
Applying the method of least squares to determine these coefficients we
minimize the sum of square of the deviations.
held fixed.
5. The coefficient of determination, R2, is a widely used measure of the goodness
of fit of a regression line. It measures the extent or strength of the association
that exists between the two variables, X and Y. Since this value is based on NOTES
the sample data points, it is known as the ‘Sample Coefficient of
determination’. It’s value ranges from 0 (poor fit) to 1 (good fit) or (perfect
fit).
6. In case the dependent variable is a function of two or more independent
variables, i.e., y = f(x1, x2, ... xn) type, this type of regression is called as
multiple regression.
7. The line of regression is the line which gives the best estimate to the value of
one variable for any specific value of the other variable. Thus, the line of
regression is the line of ‘best fit’ and is obtained by the principles of least
square.
Short-Answer Questions
1. What is meant by testing of hypothesis?
2. What is a null hypothesis?
3. What are type I and type II errors?
4. What are one sided test and two sided test?
5. When a sample is called small?
6. What do you understand by regression analysis?
7. What is non-linear regression?
8. Explain briefly coefficients of determination.
9. What do you understand by Standard Error of Estimate?
Long-Answer Questions
1. A company claims that 5% of its products are defective. In a sample of
400 items 320 are good. Test whether the claim is valid.
2. From rural schools 700 children were randomly chosen and examined.
200 of them had eye defects. From urban schools out of 900 children 180
had similar defects. Examine whether more percentage of rural school
children have eye defects.
3. A group of children was tested for measuring mathematical ability. They
were given tuition for a month’s time and second test of equal ability was
given.