UC-NRLF
$B 531 D3D
£Mnbuvob nDatbematical XTvacts
INTERPOLATION
AND
NUMERICAL INTEGRATION
A COURSE IN
INTERPOLATION AND
NUMERICAL INTEGRATION
for the Mathematical Laboratory
by
DAVID GIBB, M.A., B.Sc,
Lecturer in Mathematics in the University of Edinburgh
Xonbon:
G. BELL & SONS, LTD., YORK HOUSE, PORTUGAL STREET
1915
Digitized by the Internet Archive
in 2008 with funding from
IVIicrosoft Corporation
http://www.archive.org/details/courseininterpolOOgibbriGh
PREFACE
nPHE present work is intended for the use of students who
are learning the practice of Interpolation and Numerical
Integration. The advantages of a practical kno wedge of this
part of Mathematics ai-e so obvious that it is needless to insist
on them here : and these subjects form an important part of the
course in the modern Mathematical Laboratory. There are,
however, so many claims on the time of students that the extent
of this course, as of all others, must be kept within narrow
limits : and it has therefore been necessary to restrict the treat-
ment to the most central and indispensable theorems. A large
number of numerical illustrations and examples has been given
of a kind likely to occur in the applications of Mathematics.
My thanks are due to Professor Whittaker and to my
colleague, Mr E. M. Horsburgh, M.A., B.Sc, Assoc. M.Inst. C.E.,
for their valuable criticisms and suggestions.
D. G.
The Mathematical Laboratory,
University of Edinburgh,
July 1915.
325S65
CONTENTS
PAGE
Preliminary Note on Computation - - - i
CHAPTER I.
Some Theorems in the Calculus of P^inite
Differences.
Section
1. Differences ------- 3
2. Properties of the Operator A - - - - 5
3. Expansion of /(.v + z/tc) for Integral Values of;? - 7
4. Relations between the Differences and the Differential
Coefficients of a Function . . . . g
Miscellaneous Examples - - - - - 12
CHAPTER II.
Formulae of Interpolation.
5. Introduction ------- 14
6. Lagrange's Formula of Interpolation - - -15
7. Periodic Interpolation - - - - - 18
8. Notation ------- 18
9. Newton's Formula of Interpolation - - - 19
10. Stirling's Formula of Interpolation - - - 24
11. Gauss' and Bessel's Formulae of Interpolation - - 26
12. Evaluation of the Derivatives of the Function - - 29
Miscellaneous Examples - - - - -33
CHAPTER III.
Construction and Use of Mathematical Tables.
13. Interpolation by First Differences - - - - 3^
14. Case of Irregular Differences - - - - 3^
15. Construction of New Tables - - - - 39
16. Tabulation of a Polynomial - - - - 39
VUl CONTENTS
Section page
17. Calculation of the Fundamental Values for a Table of
Logarithms - - - - - - 40
18. Amplification of the Table by Subtabulation - - 43
19. Radix Method of Calculating Logarithms - - 46
20. Inverse Interpolation - - - - - 51
21. Various Applications of Inverse Interpolation - - 58
Miscellaneous Examples - - - - - 61
CHAPTER IV.
Numerical Integration.
22. Introduction - - - - - - - 63
23. Evaluation of Integrals by the Integration of Infinite
Series Term by Term - - - - - 64
24. A Formula of Integration based on Newton's Formula - 65
25. Case when the Upper Limit is not a Tabulated Value - 70
26. A Formula based on Bessel's Expansion - - 73
27. The Euler-Maclaurin Expansion - - - - 76
28. An Exceptional Case - - - - - 78
29. The Formulae of Woolhouse and Lubbock - - 79
30. Formulae of Approximate Quadrature - - - 81
31. Comparison of the Accuracy of these Special Formulae 85
32. Gauss' Method of Quadrature - - - - 86
Miscellaneous Examples - - - - - 89
PRELIMINARY NOTE ON COMPUTATION.
DEFORE introducing the formulae required in interpolation
and numerical integration, it may be advisable to mention
a few points which will facilitate the work of computing.
Where computation is performed to any considei'able extent,
computer's desks will be found useful. Those used in the mathe-
matical laboratory of the University of Edinburgh are 3' 0" wide,
r 9" from front to back, and 2' 6|" high. They contain a locker,
in which computing paper can be kept without being folded, and
a cupboard for books, and are fitted with a strong adjustable
l)ook-rest. Thus the computer can command a large space and
utilize it for books, papers, drawing-board, arithmometer, or
instruments. Each desk is supplied with a copy of Barlow's tables
(which give the square, square root, cube, cube root, and the
reciprocal of all numbers up to 10,000), a copy of Crelle's multi-
plication table (which gives at sight the product of any two
numbei's each less than 1000), and with tables giving the values
of the trigonometric functions and logarithms. These may,
of course, be supplemented by a slide rule, or any of the various
calculating machines * now in use, and such books of tables as
bear particularly on the subject in question.
Success in computation depends partly on the proper choice of
a formula and partly on a neat and methodical arrangement of the
work. For the latter computing paper is essential. A convenient
size for such a paper is 26" by 16"; this should be divided
by faint ruling into ^" squares, each of which is capable of
* For descriptions of these see Modern Instruments and Methods of Cal-
cidation, edited by E. M. Horsburgh. London : G. Bell & Sons. 1914.
e. 1
,2 INTERPOLATION AND NUMERICAL INTEGRATION
holding two digits. It will be found conducive towards accuracy
and speed if, instead of taking down a number one digit at a time,
the computer takes it down two digits at a time, e.g., instead of
taking down the individual digits 2, 0, 4, 7, 6, 3, it will be found
better to group them together, as 20, 47, 63. Every computation
should be performed with ink in preference to pencil ; this not
only ensures a much more lasting record of the work but also
prevents eye-strain and fatigue.
It need hardly be emphasized that seven-place accuracy cannot
be obtained by the use of four-figure tables. A fairly safe rule is
to make use of tables containing one digit more than the accuracy
requii-ed ; to use more accurate tables than these is simply to
increase the amount of labour without increasing the efficiency.
For the same reason contracted methods of multiplication and
division should be adopted when machines and tables are not at
hand.
CHAPTER I.
SOME THEOREMS IN THE CALCULUS OF
FINITE DIFFERENCES.
1. DiflFerences.
Let \f{x) denote the increment of a function f{x) correspond-
ing to a given constant increment w of the variable a-, so that
^f{x)=f{x + w)-f{x).
The expression l^/{x), which is usually called the First Differ-
ence of the function f{x) corresponding to the constant increment
w of the variable x, will, in general, be another function of x, and
as such will have a first difference which will be obtained by
operating on A/ (a-) with \ Calling the result of this operation
A^f{x), we have
AV(x) = A{A/(x-)}
= A{f{x + w)-/{x)}
= {/{x + 2w) -f{x + w)}- {/{x + ic) -f{x)]
=f{x + '2w)-2f{x + w)+f{x).
This is termed the Second Difference of f{x) corresponding to the
constant increment w of x. By repeating the process we obtain
in turn ^^f{x), ^^f{x), ... i\'\f{x), which are called the 3rd,
4th, ... nth differences of /{x) corresponding to this constant
increment 7V.
Xote. — It is always possible by means of a suitable trans-
formation to make the value of w unity.
For since w is the increment of x in f{x), then unity is the
increment of — . Replacing x in f(x) by wy, we obtain a new
w
function F {y), which is such that
^F{y) = F{y + \)-F{y).
INTERPOLATION AND NUMERICAL INTEGRATION [ch.
Example 1. — Tabulate and difference the successive values of f[x):
from a; = 110 to a; = 118.
A3 A*
X
7?
A
A2
no
1331000
36631
111
1367631
37297
666
112
1404928
37969
672
113
1442897
38647
678
114
1481544
39331
684
115
1520875
40021
690
116
1560896
40717
696
117
1601613
41419
702
118
1643032
In the above table the column headed x^ was obtained from Barlow's
Tables. The column A gives the first differences, being obtained by sub-
tracting each entry of the column x^ from the entry following it, e.g. if
x=114, a;+l = 115, A/(a:) = (a;-i- l)3-a;3= 1520875- 1481544 = 39331. This
result is set on a line midway between the numbers corresponding to a;=114
and a;= 115. In the same way the 2nd differences, A-, are obtained from the
1st differences, the results being set on a line midway between the latter,
and so on.
It will be noticed that the 3rd differences of ar* are constant and that its
4th differences are zero. It will be shown later that the nth differences of a
function of the nih. degree are constant and that the differences of higher
order are all zero.
Example 2. — To determine the successive differences of sin (ax + 6).
Let /(a;) = sin(aa:-l-t).
Then A/{a:) = sin {a x+w-vh) -sin [ax-vh]
= 2 sin ^a IV cos (ax + b + ^aio)
: 2 sin ha w sin (ax + b + haw + ir)
A-'/{x) = A{ 2 sin I a w sin {ax + 6 + ^ a ?t' + tt)}
= 2 sin ^ as
= {2 sin haw}- cos {ax + b + ^2aw + v)
ia w { sin {a x + io + b + ^aw + ir) -sin (a x + b + ^ a w + w)}
■{2smhaiv}"am (ax + b-\- ^2aiv + 2 ir).
1, 2] THEOREMS ON DIFFERENCES
Similarly, it can be shown that
A»/(a;) = { 2 sin ^awf sin{a x + h + ^3aw+^Tr)
and, more generally, that
A"f{x) = { 2 sin law}" sin {ax + b + ^n aw + 7nr).
This can be proved by induction in the following way :-
A"+i/(a;) = A [{ 2 sin | a w}" ain {ax +b + ^naw + n ir)']
: { 2 sin i a m; }" [sin (ax+w + b + ^naw + mr)
- sin {ax + b + ^naw + mr)]
= { 2 sin I a iv }"+i cos{ax + b + ^{n + l)aw + mr)
= { 2 sin i a «}«+i sin {a x + b + ^ {n+ I) a lu + {n+ I) n).
But the theorem has been proved for n—l, 2. Hence it is true always.
Example 3. — Find the first differences of the functions
cos {ax + b) ; tan— ^ ax ; «''•"=.
Example 4. — Find the nth differences of the functions
— ; a-' ; cos^(a.v + 6).
X
2. Properties of the Operator ^.
The operator A obeys certain of the fundamental laws of
algebra, viz. —
(i) The Law of Distribution.
For d {fix) + ^ {x)) = {/(x + IV) + c/, (x- + w)) - {/(.«) + c/) ix))
= [fix + w) -f{x)] + {c^{x + w)-<l> {x)]
= A/(a;) + Ac/,(x).
(ii) The Law of Com mutation if the coefficients of the various
terms are constants.
For, if a is a constant,
A [a .f{x)] = af{:c + w) - af{x)
= a{f{x^-w)-f{x)]
= a.\f{x).
(iii) The Law of Indices.
For A'- . \^f{x)
== {A . A . A ... (?• factors) .AAA ... (s factors)} /(a;)
= {A . A . A ... (r + s factors) }/(x-)
The indices r and a" are, of course, integers : the operator
A" was defined in the previous section only for integral values of ii.
6 INTERPOLATION AND NUMERICAL INTEGRATION [ch. i
Example 1. — If f{x) is a polynomial in x of the nih. degree, then its nth
difference is constant and its differences of higher order than n are zero.
Let f(x) — anX'^+an-\x^-'^+ ... +aix + ao
in which the coefficients a„,an-i, . . ai, ao are constants.
Then Af{x) = {a„ {x + lo)" + «"-i (x + w)"-^ + ... + ai (x + ir) + a^, }
— {a„ X" + a,i-i x"-'^ + ... + ai a; + ao }
= 7t an x"-^ ir + h„_.2 X"-- + ... +bi x + h„
where 6„_2, ..., ^i, 60 are constants.
We thus see that A/(x) is a polynomial in x of degree (71 -1). Repeating
the operation, we have
A^f{x) = {n a„ [x + »')"~^ "' + t„_.j {x + *")"-- + ... + ?>i (x -^ »•) + ?>,j}
-{nanX^'-'^ ii; + h„-2X"--+ ... +h.^x-\-h(^}
= n(,n-\)a„x"'-->if- + Cn~'iX'^-^-\- ... +Cia; + Co
where c„_3, ..., q, c^ are constants.
Hence A-/(x') is a polynomial in x of degree (/i-2). It is evident that
by continued repetition of this process the degree of the function will be
reduced to zero, and that we shall obtain
A«/(a;) = »i (vi - 1 ) (?j - 2) . . . 3 . 2 . 1 a„ »•" ,
i.e., the ?ith difference of the function is a constant. Obviously the next and
all differences of higher order will be zero. This result is of great import-
ance ; probably no other theorem will be applied as frequently as this. For
instance, it enables us to complete a mathematical table in which there are
gaps. This is illustrated in the next example.
Example 2. — If
/(a:) = 0-5563, /(x f l)-0-p682, /(x + 3) = 0-5911, /(a; + 5) =^0-6128,
obtain /(a; + 2) aad/(a; + 4) on the assumption that /(x) is a polynomial of
the third degree.
It has been shown in § 1 that
^""/{oc) =/(.T + 2 w) - 2f[x + w) +f(x).
By repeating the process we get
A'f{x)=f{x + 4.ir)-4:f(x + 3w) + 6/ix + 2,r)-4/(x + >r)+f(x).
Bat the given function is of the third degree. Hence for it we have
A-'/(.x-) = 0. Noting that w is unity we obtain
f{x + 4) -■i/(x + 3) + &/{x + 2)-i/{x+l)+/{x) = 0
or, denoting /(x + Ji.) by/,,,
/4-4/; + 6/,-4/.-f-/o = 0.
Similarly /o - 4/, + 6/ - 4/, 4-/1 - 0.
2, 3] THEOREMS ON DIFFERENCES 7
Substituting the given values in these two equations, we get
/4 + 6/3 = 4 -0809
/4+ y;- 1-1819
whence /, = 0-5798 , /j = 0-6021 ,
i.e., /(X- + 2) = 0-5798, /(x-r 4) = 0-6021.
Example 3. — If when a; = 0, 1, 3, 4, 6, the corresponding values of
f {x) are 1, 1, 79, '253, 1291, 6nd approximate values for f{x) when x = 2
and when .1=0. Under what conditions will the result be accurate ?
3. Expansion of f{x + mv) for Integral Values of n.
We shall now establish certain formulae which express
J{x-\-mv) (where n is an integer) in terms of differences of
various orders. These formulae will be shown in Chapter II, to
be valid under certain conditions even when n is not an integer.
It will be found convenient to introduce another operator E
defined by the equation
E=i+x y
This new operator will obviously obey the same laws as A.
Now EJ (x) = (1 + ^)/{x) = f{x) + \f{x) ^/{x + w).
Thus the effect of the operator E on the function f{.v) is to increase
its argument x by tlie given constant quantity w. Repeating the
process we get
E"J{x)=f{x + nw)
and since the operator sS. may under the conditions of § 2 be treated
as an algebraical quantity, this last equation may be written
f{x + nw) = {l + ^Yf{x)
=f{x)+,fi,^f{x) + „c,^''f{x)+ ... +^'\f{x) (1)
where ,/.v denotes the coefficient of a'' in the binomial expansion of
(1+a)". The function /(a; + 7iM/-) is here expressed in terms of
f{x) and its successive differences.
It is evident that, since the effect of the operator E on f{x) is
to increase x by w, the result of operating on f{x) with E~'^, where
E~^ is defined by the equation EE'^^^1, will be to subtract w
from X,
i.e., E-^/{x)=^/{x-w).
8 INTERPOLATION AND NUMERICAL INTEGRATION [cu. i
Hence the above formula (1 ) may be transformed as follows :
J {x + nw) =/{x) + ,.ci ^f^x) + {l + ^) E-' {,c, A- fix) + „c, ^y(x-)
+ ... +S"/{x)}
=/(•«'•) + ,fii ^J (a) + (1 + Aj {„c, A- f{x - iv)
+ ^^Cs^'/{x-7J))+ ... +1"/{X-W)}
+ (l+A)^-M„+iC,Ay(x-«;)+ ... }
+ (l+A){„^,c,Ay(x-2u;)+ ... }
+ ,,+1^4 AV'(« - 2^) + „+oCg A'f{x - 2m;) + ...
Operating next with (1 + A)^-' on the 7th and subsequent
terms, then on those beginning with the 9th term, and so on, we
obtain
f{x + nw) =/{x) + „ci ^/(x) + „c., /\'/{x -w) + „+iC3 Ay (a; - w)
+ „+,c, A V(«-' - 2*«) + „+.c, Ay (X -2w)+... (2)
the formula ending at the term .X-"~^/(x - n - I w).
In the same way it may be shown that
f{x + n w) ^f{x) + „c, Ay'(A' - w) + ,,^^c. A-/(^; - lo)
+ „^,c, AV(-*; - '^^) + ,,+.^4 A-* ;\x - 210)
+ „^.,C5 A'/'(x - otv) + ... {3}
the formula ending at the term A-'"y"(;f - n w).
Changing x into x + w and rt into n-\ in (3) we have
J\x + 11 w) =/{x i-w) + „_iCi \f{x) + ,fi., A-y"(a;) + „c. A"/ (a; - w)
+ „+jC, AVXa; -w) + „+,c, A' fix -'2w) + ... (4)
the formula ending at the term A-"/(x -n-2 w).
J, 4] THEOREMS ON DIFFERENCES
Adding (2) and (3) we obtain
f{x + n lo) =f{x) + ,fi, ,^ ^ ' + — ,fi, A-7 (,x - w)
A\f{x-w) + ^'/(x-2w)
+ ^„^,c,i\Y{x-2w) + ... (5).
Similarly, addition of (2) and (4) gives
Formulae (1), (2), (3), (4), (5), (6) may be compared with the expansions
introduced in Chapter II. under the names Newton's, Stirling's, etc.,
Formulae of Interpolation.
Again, from the definition of the operator E, we see that
Ay\x) = {E-\r/{x).
Expanding the operator in this last relation we obtain
A"/{x) = E''/{x) - ,fi, ^-VX*-) + .^. E^'-'fi^) ••• + ( - 1)V(^-)
=f{x -\-nw) - „Ci/(x- + n-lw)+ „Co/(.r + w - 2 w)
+ ... +{-\rf{x).
In this result we have the 7ith difference of the function f{x)
expressed explicitly in terms of f{x) and its successive values.
4. Relations between the Differences and the Dif-
ferential Coefficients of a Function.
The second example of g 1 may ha\ e reminded the reader of a
somewhat similar result in the diflerential calculus, viz.,
-— sin (a x + h) = a" sin (« x i-b + hn tt).
dx
This result may be derived from the example in question. For if
we put iv = .\x, where Ax' denotes a small inci-ement in x, and,
10 INTERPOLATION AND NUMERICAL INTEGRATION [ch. i
after dividing both sides of the equation by (A x)", make A x tend
to zero, we get
d'\f {^) T (-2 sin ia . Ax^" . , . , . , ..
^d^ = lu['""tx~ ) sm{a*' + 6 + l7.(a.Aa; + 7r)}.
Aa;->0
= a" sin {ax + b + J n tt).
This naturally suggests that there may be some relations connecting
the successive differences of fix) and the successive differential
coefficients of the same function. That such a relation does exist
we now proceed to show.
Applying Taylors Theorem* to the function /(a; + w) we obtain
f{:x + w) =/{x) + wf (a;) + ^y"' (.r-) + ... + ^./■'"" {^^ + • ■ •
or, since f{x + id) -J\x) = \f{x),
A/{x) = w/'{x)+ '^^/"{x)+ ... + ^f-'(x) + ...
which is a relation of the kind sought for.
Again, by Taylor's Theorem,
./(^- + "')=/ (-»•)+ wj'{x)+ — _/'(y/)+ ... + --/'""(.'/;)+...
'J> ! in\
f(x+2to)=/{x) + 2w/'{x)+^^^/"{x)+ .. +l^'/<™)(a;) + ...
and generally
/{x + n 2v) =/ (*) + n to/' (x) + ^' /" (x) + ...+ ^' f"" (o:) + ...
Hence
'^"/(•^-') =/(''' + '* '^') - »Ci /{-^ -t- w - 1 tv) + „c,/{x + 11-2 iv)
- ... +{-\)y(x)
= 2 -^ {"'" - „Ci {n - 1)'" + ..Co [n - 2)'" - . . . }/'"" (x).
* It would be inopportune at thi.s stage to introduce a discussion of the
questions of convergence whicli arise in connection with the application of
Taylor's Theorem. We shall assume a convergence sufficient for our
purpose.
4] THEOREMS ON DIFFERENCES 11
But
^"' ~ «Ci (w - 1 )'" + ,fi2 ('«■ - 2)'" + ...
= 0 if m<7i*
= nl if ni = 71
= i ri (w + 1 ) ! if m = ?i + 1
= — ?i(3« + 1) (?i + -2) ! ifm = 7i + 2
etc.
Making use of these results in the expi'ession for A"/(x) we find that
15 7r (?i + 1) „ . ,„
+ ^^ 'w"^'J^'^+-{x) + ...
Corollary 1. — Since the rzth differential coefficient of a function
_/'(x') of the nth degree is a constant, and the derivatives of higher
ordei- ai-e zero, we see again from this expression for A"/{x) that
the ?ith difference of a polynomial of degree ?*. is a constant, and the
differences of higher order zero.
Corollary 2. — By reversing the series for Z\"/"(.^'), we obtain an
expansion for the nth differential coefficient of /{x) in terms of its
differences, viz.
lon{n + 2) {n + 3) .... . ^
6! - "-^^'
CoruUary 3. — If the ?«th differences of a function are constant,
then the function is a polynomial of the nth degree.
For since the nth differences are constant, the higher differences
are zero, and therefore the nth. differential coefficient is constant,
Chrystal's Algebra, Fart ii., Chap, xxvii., § 9.
12 INTERPOLATION AND NUMERICAL INTEGRATION [ch. i
and the differential coefficients of higher order zero. Thus the
function, when expanded by Taylor's Theorem, takes the form
/(,«) =f(X) ^(x- x)f (X) + '•' ^ f^\r (X) + ...
.'-^^■,/ (.r, "
n !
and this is of the ?ith degree in x.
Corollary i. — Since e^ = 1 + j-: + — + ... -\ . + •■•
we see that the series which represents /{x + w), viz.
f{x) + ivf'(^x)+ y,/"(a^) + ...
may be symbolically represented as
d
e"~^f{x).
Hence Ef{x) = e'^f{x).
MISCELLANEOUS EXAMPLES.
(Asstime in the followiny exercises that w is unity, unless otherwise indicated).
1. Find the first differences of the functions
tail a X ; log x ; ( a-^ - a* ) / (« - 1 )•
2. Prove that
sin 0
A tan a; 0 =
eosx d cos(a,-+ \)d'
3. Prove that
1 1 e
^ X"'"^^ taii^;^ .
2^ tan—
4. Find the u'-'' differences of the functions
x+\ ^ . ,
—- _ — . x™+" ; sin ax cos ox.
X- \Zx '
5. Find the »ith difference of each of the factorials a;(a; + a) (x + 2a)... and
1
a;(a; + a)(x + 2a).
there being m factors in each factorial and the increment
4] THEOREMS ON DIFFERENCES 13
of X being a: and show how to develop any polynomial in x in a series of
factorials.
In particular show that
x^ =x^lx{x-\)^^x{x-\){x-'l)-Vx{x-Y){x~1){x-'>,).
6. If f(x) = a?,^ + 6 2^- prove that
AV(a;)-3A/(x') + 2/(x) = 0.
7. Prove that
A X" - 2A--^ x" + 3 A3 a;" - = (a- ^ 1 )" - (.r - 2)" .
8. Establish the formula for tJ e Jith difference of the product of two
functions
A« {f.x) . 4> {X) } = A«/(.T) • 0(-i-) + n A«-i/(.r + 1 ) . A 0 ( .r)
7i(n-l)
+ ^yy- A«- V(a' + 2) . A'-^,^(a-) + . . .
9. Show that
0 (E) a\f{x) = a^ <p (a E)f[x)
and find the value of
AV(-«) + — (a - cos 6) A/(.r) + ( 1 _ ^ cos 6 + ^ )f{x)
a \ a a- /
cos h X
where f{x) = — -- — .
10. The record of exports for the year 1813 was destroyed by fire. Make
an estimate of their value given that the values of the exp(rts for the years
1810, 18 11,. 1812, 1814, 1815, 1816 were 48, 33, 42, 45, 52, and 42 million
pounds respectively.
11. It is asserted that a quantity, which varies from day to ('ay, is a
rational and integral function of the day of the month, of less than the fifth
degree, and that its values on the first seven days of the month are 30, 30,
28, 25, 22, 20, 20. Examine whether these assertions are consistent. If so,
assume them to be true, and find (1) the degree of the function, (2) its value
on the sixteenth day of the month.
[CH. II
CHAPTER II.
FORMULAE OF INTERPOLATION.
5. Introduction.
The problem of determining the value of a function for
some intermediate value of the independent variable or argu-
ment, when the value of the function is given for a set of discrete
values of the argument, is called Interpolation. When the function
obeys some definite law the approximation to the true value may be
made as accurate as may be desired : but if, as often happens,
the rigorous analytical formulae which occur are very com-
plicated, they will be of little value from the point of view of
computation. It will then be necessary to use a method analogous
to that which has to be employed when the properties of the
function are entirely unknown, save for the values (usually obtained
as a result of experiment or observation) which it has for a certain
limited number of values of the argument. It is to this case that
our attention will be directed.
When the observations are given for values of the argument
equally distant from each other, and are in arithmetical progression
so that they can be plotted on a sti'aight line, the possibility of
interpolating is obvious. But, in general, the observed vahies
do not follow any apparent law, and unless some assumptions
are made regarding the form of the function, the problem is
indeterminate. It is, however, customary to obtain the entry for
values of the argument not far distant from each other ; and when
this is done, in most practical problems it is found that all the
differences at some particular order become sensibly zero. In such
cases the observations may be represented with considerable
accuracy by a polynomial in x of the form
y = flo + «i a; + «,; »■■ + ... + a„ a?".
This is the equation of a parabola of the nth order, on which
account the method now to be described is sometimes called
Parabolic Interpolation.
5, 6] FORMULAE OF INTERPOLATION 15
In some cases it is advisable to calculate the function not
directly, but in combination with some other function. The
calculation of the values of the probability integral e-i^" dt
furnishes an example of this kind, since, when the argument becomes
large, the integral and its differences become inconveniently small
and irregular. In this case it will be found convenient to compute
the function ei^' e'^^'dt. At other times some function of the
given function is calculated ; for example, if the observations
formed a geometrical pi'ogression, then we should prefer to make
the interpolation on their logarithms rather than on the quantities
themselves.
Wallis (1616-1703) may be said to have laid down the principle of inter-
polation in his approximation to the area of the circle j/ = (,r- t-)^ '-. In this
work he first assumed that the value of the integral I (x - .i^)^ ' dx might be
Jo
taken as the geometric mean of the values of
x")dx, i.e. of 1 and ^.
Subsequently he saw that this was not exactly true, and that the value of
I (.r -a;-)^ - dx must obey the law expressed by the series
I (x - x')" dx and } {x ■
le saw that this was not
; must obey the law expr
j {x - xY dx, {x- X-) dx, ... , I (x- x~)" d.i
this was equivalent to interpolating the value of the integral under con-
sideration.
6. Lagrange's Formula of Interpolation.
We shall first show how to find a polynomial of degree (n - 1)
at most, which has the same values as some arbitrary function /(oc)
for n distinct values of x, say x = a^, «o, ... a„. This polynomial
will be represented by the parabola of order {n - 1) already men*
tioned, and may be expressed in the form
/(x) = Cg + C^X + CoX-+ ... + c„_i x"~\
But it will be found more convenient to write it, as is always
possible, in the form
fix) = 2 A^{x-a-^) ... (x - a^_,) (x - a^^-^) ... (x - a„) (1 )
16 INTERPOLATION AND NUMERICAL INTEGRATION [ch. ii
When ,r is put equal to o, all the terms on the right-hand side of
(1) vanish except the tei'm containing A^: at the same time /(x)
becomes f{a^). Hence
^r =J ("r) / («r - «i) • • ■ («. - «,._,) («, - «,+,) ... (a, - rt,,)
and equation (1) becomes
fir) = y (a^ - CTi) ■ • ■ (■« - «,_,) (a? - a,+i) . . . {x - a„) _^
^"^ rfi K-cii) ... («,-a,._j)K-o,^,) ...(«,.-«„) ■^*"'-^ ^-^
This formula expresses /(x) for all values of x in terms of its
values at the n points «„ o.., . . . a„. If /(a?) is a polynomial of
degree less than n, the formula gives its accurate expression. But
\if{x) is not a polynomial of degree less than n, then it is no longer
possible (in default of other information) to determine f{x) from a
knowledge of its values at the points flj, a„_, ... a„ ; for if any one
such function were known we could obtain another by adding to it
any function which is zero at the points a^, a.., ... a,„ e.g., we could
add the function {x - a^) (x - a.,) ... (x - a„) (f) (x) where (f> (x) does
not become infinite when x = a^, o.,, . . . a„. If, however, /(x) is a
function of which nothing is known except that it takes certain
definite values corresponding to certain definite values a^, a.., ... a„
of the variable, then the simplest assumption which we can make
regarding its form is, that it is the function of degree (n - 1) which
sa,tisfies the given conditions ; and this function is given by
equation (2), The expression on the right-hand side of (2) may
therefore be taken as the expression for /{x) for all values of x in
the interval considered. It may be written in the form
fix) = Z -, , , , J (a.)
,.^1 (x- a,.) (b (a^)
where (/> (x) = {x - a^) {x - a^) ... (x - a„).
This is known as Layrange^s Formula of Interpolation. The
advantage which it has over other formulae for the same purpose
is that it is in a form suitable for logarithmic computation.
This latter form of Lagrange's formula may be obtained directly from the
formula previously given. It can also be obtained very readily by the
Method of Partial Fractions,
/(£l
<t> {X)
<t> {x) = [x - a,) [x - ff.,) ... (.T - a„ ).
For consider the proper fractional function , , . where
6] FORMULAE OF INTERPOLATION 17
Decomposing it into partial fractions we have
<p [x] X - ttj X - a., X - a„
But, by the general theory of partial fractions, the coefficient A,- is given by
f(a,.\
^'■- <p'{ar)
Hence
/{■r)_ /K) , /(0.3) I , /( an )
(p(x) (X - a^) (p' (a-^) {x - ao) (p' (ao) '" (x - a„ ) 0'(a„ )
which gives
0 (X)
^ ^-'"^ " S (*•-«,.)</.'(«,.) -^ ( ^-^ )•
Examplt /.—Construct a polynomial of the third degree in x which
shall have the values -5, 1, 1, 7 when the values of x are -2, - 1, 1, 2
respectively, and determine its value when x is %.
Let f{x) be the required function. Substituting the given values in
Lagrange's Interpolation Formula we have
(.r+l)(.r-l)(a'-2) (■r + 2) (.r - 1) (.r - 2)
/(•'■) - (_2+l)( -2- 1)(- 2-2)^ ^'^ (-l+2)(-l-l)(-l -2)^ ^^^
(X + 2) (.X- + 1) (.v - 2) , (X + 2) (.r + 1) (a- - 1) „
(I + 2) (1 + 1) (1 - 2) ^ (2 + 2) (2 + 1) (2 - I)
= tVU-+l) (-f-l) (.r-2) + i(,r + 2) (.r-1) (.r-2)
- 1 [X + 2) (.r + 1 ) (,r - 2) + /. (.r + 2) (.r + 1) (x - 1)
= ^-x\\.
From this it follows that /(:') = 2|. But in practice it is not necessary to
obtain the algebraic form of/(.r) if we merely require its value for some
value of X. Direct substitution in the formula will give the result much
more quickly. In the question under consideration
(Hl)(|-1) (i^-2) (1 + 2) (§-1) (§-2)
./ (:')-(_2+l)(-2-l)(-2-2)^ 5) + (_i^.O)(_i. i)(_i_2)
(f + 2) (1 + 1) (§-2) (5 + 2) (f + 1) (§-1) „
+ (1+2) (1 + 1) (1-2) + (2 + 2) (2+1) (2+1) ( + ''
2J. 7_ J. 3 5. I 345
~ 96 48'48~»6
'8*
Example ^. — Find the minimum value of the polynomial of the second
degree which has the following values
.r 0 , 1 , 5
/(x) 11, 5, 21.
G. 2
18 INTERPOLATION AND NUMERICAL INTEGRATION [ch. ii
Substituting these values in Lagrange's Interpolation Formula we have
., , (.r- l)(a--5) x(x-5) x{x~l)
f(x) =-V-(2a:-6) - |(2a;-5) + t^(2.r-11.
For a minimum turning value /' {x) = 0
i.e., 44(2.r-6)-25(2,T-5) + 21(2a;-l) = 0
SO.r- 160 = 0
,. _2
Hence the minimum value of /(,i-) corresponds to the value Q of the argument.
The corresponding value of f{x) is 3.
Example 3.—lif(x) is a function of x which assumes the values 20, 35, 56, 84,
when X has the values 0, 1, 2, 3 respectively, express f{x) in terms of x, on
the supposition that A*f(x) = {). Also find an approximate value of x if
/(.t) = 25, third differences being neglected.
Example 4.— If four consecutive values of a function corresponding to the
values 1, 2, 3, 4 of the variable be 601, 777, 902, 999 respectively, find the
value of the function when the argument has the value 1-5.
7. Periodic Interpolation.
When there is reason to suppose that the function is periodic, it is better
to assume for /(a;) an expansion in terms of periodic functions, say,
f{x) = «() + aj cos X + ttg cos 2x+ ... + a„ cos n x
+ 61 sin .1- + 60 sin 2 .T + ... +6„sinn.T,
the number of unknown coeiEeients corresponding to the number of observa-
tions given. When the observations are given for values of the argument
spaced at equal intervals, the problem becemes one in Fourier's analysis and
will be found treated in the tract on that subject. But when the values of x
are not in arithmetical progression we may use a formula which is analogous
to Lagrange's formula of interpolation, and which Gauss has shown to be
equivalent to the above expansion ior f{x).
If/(,r) is given for x-Xf^, .r,, Xo, ... x^n the corresponding formula is
ff^\ - '^ sm^{x-X(t) ... sin^(x-Xr-:)am^{x-Xr+i)... sin i (.i- - .r2„) _
'■-0 Sin J (xr - .Ti) ... sin 4 {.r,. - Xr-i) sin | {xr - av+i) ... sin ^ (.r,. - X2„)
8. Notation.
Hitherto the points a^, a„, ... a„ have been taken to be any
points whatever. But in the most important class of cases it is
possible to obtain the entry for values of the argument equally
distant from each other. In what follows we shall assume that the
r. 9] FORMULAE OF INTERPOLATION 19
arguments are spaced at equal intervals, and shall use the notation
of the following scheme which for many purposes * is more con-
venient than that used in Chapter I.
A6
Argument.
Entry.
A
A'
A^
A^
A5
a — o ?r
./■-;i
^Ai
a - 2iv
./'-2
IJ- 8f- 2
8V-
2
l^f-^
S./"-^
/.8V_
V 8\f-
a - ir
./'-I
i"-s/'-l
s-/-
1 /^ SV-
1 sy-i
/^./■-i
8/-i
f, sv-
i 3V-
i /^^V-i
8V-k
a
A
^s/;
8Vo
/* ^Vo
sv;
l^8V.
i\f\
Vi
H-o~A,
•^Vi
/-sVi
^A
a + n-
A
1^ s./'i
ov;
/^ ^v;
8'A
i^A
3/5
/^svi
sy,
a + 2 w
A
i-A
1^8 A
8A
8%
a + 3 IV
A
whe^i-e
8Ai
=A-A-
1
/^./i = h
(A + Ai)
s./l
=A-A
/xS/,= l
m + 8A
i)
etc.
etc.
8%
9. Newton's Formula of Interpolation.
If we assume that the differences n.^ - a^, a.-a^, ... in Lagrange's
interpolation formula are all equal to iv, and that a^ in the formula
corresponds to a in the scheme, then the 'quantities /(a^) take the
following forms : —
./■(«j=./;=/o+3.4
f(ao)-A=A'+8A =A+-^8j] + 8V,
/(«.o=y;.^3/i . "^ 8v. . ^±z}ltzl} ,j^ , ...
* Especially in connexion with central-diflference formulae.
20 INTERPOLATION AND NUMERICAL INTEGRATION [ch. ii
Lagrange's Interpolation Formula then becomes, if the number
of given points is denoted by k,
^ {x - g.,) {x - flg) (a; - g,) ■ . ■ {x - a,.)
(rti - a.,) (ai - a,) («!-«,)... (a, - a, ) ' "
(a? - CTi) (ft7 - CTo) {x -a^) ...{x- a,.) , . , . ,
+ / w V, ^ ; ; C/o + (^Jh )
{a.2 - «]) {a., - tto) (a„ - a^) ... (Oo - a^.) ■'
, (..-aOU--)(.^ -«.)...(--«.) ^j-^,^j-_,^^,.^
{as - a,) (a„ - n..) (a. - rt^) . . . (03 - %)
(x - CTi) ■ ■ ■ (.r - a^) (x - g^+o) ... (x - a^)
(g,.+i - gj) ... (gv+i - gj (g^+i - g,+o) ... (g,.+i -
where g^ = g, g., = g + ?<;, . . ,, g^,. = g + (A; - 1 ) iv.
This is an expression of degree (^ - 1) in a?. In it the coefficient
oi /q is unity when x = a^, a.2, ... a^,, and hence is unity always.
Again, the coefficient of S/1 is 0 at a^, 1 at gg, 2 at gg, ... r at
g,.+i, ..., and as this coefficient must be a polynomial in x of degree
less than k, it must be no other than . In fact the terms in
8/^ in the above series simply constitute the expansion of
S/. by Lagrange's formula. We shall denote bv w, so
w '^ "^ w ■^
that w is a real quantity, not necessarily an integer. The
coefficient of 8f, is thus n. Similarly the coefficient of S-/J is 0
r(r-l)
at gj, 0 at g^, 1 at gs, ... — - — — at g^+i . . . , and in the same way
we see that this coefficient is — — - for any value of n. The
other coefficients can be determined in the same way, and we
thus finally have Lagrange's formula transformed into the equation
/(g + «M;)=/o + riS/i + S-/ + ' ^/^ + ...
9] FORMULAE OF INTERPOLATION 21
This is called Neicton's Formula of Interpolation* In it the
diifei-euces are those which occur along a line parallel to the upper
side of the triangle in the scheme. It may on this account be
termed a Forward Difference formula of interpolation, and will be
found very useful when the value of the function near the beginning
of a limited series of observations is required.
It must be understood that the function furnished by Newton's
formula is the polynomial of lowest degree which has the given
values at the k given points a, a+v:, a+'2io, ... , a-\-{k-\)w:
and the series on the right is a terminating series. If the given
values at the k given points are the values of some given function
^ {x) at the points, we cannot conclude that <^ {x) is represented by
Newton's formula unless 4> {^) is the polynomial of lowest degree
which passes through the points.
Suppose, now, that the values of a function are given not
merely for a finite set of values of its argument, but for the infinite
set
x = a + riv {r = Q,\,'2,'c>, ... ad inf.).
Then we can, by Newton's formula, construct a polynomial of
degree {k - \) which has the given values for the first k of these
values of the argument, where k is any positive integer : and it
may happen that, as k increases indefinitely, the series on the right
in Newton's formula becomes (for some range of values of a;) a con-
vergent series, having for its sum a function f{oi). The function
/(*■) will now have the given values for each of the values a + rw
(r = 0, 1, 2, 3, ... ad inf.) of the argument. There are, however, an
infinite number of functions which satisfy this condition and /(a-)
is that one of them which is the limit, when ^ -> x , of the polynomial
of lowest degree which has k of the given values.
An expression for the diff'erence between the value of the
function and the sum of the first (»" + l) terms of the limiting
form of this polynomial will now be obtained. Let
n{n-\) ... (ti- r ■\-\) . , „
t\a + mv)=S, + n^J. + ... + '—^^ 67^ + R.
Then "
r/ ^ r s/ m{m-\) ... {m-r+\)
f{a + m9v) -J^-niSji^ - ... — — 6/^ - H
Cf. § 3, Formula (1)
INTERPOLATION AND NUMERICAL INTEGRATION [en.
will be zero for ?n = 0, 1, 2, ... r and n: its derivate of order (r+ 1)
will therefore be sero for some value n of 7n lying between the
least and the greatest of the quantities 0, r, n, and hence the
residue is given by
E = „C,+,to'-+' /"■+'' (a + n' n-).
It may be remarked that Newton discovered the binomial coefficients
when studying the expansion in series of a function by interpolations. He
first of all considered the expansions of the expressions (1 -a'-)", (l-.i-)^,
(1 -.r-)'-^, ... and deduced from them the expansions of (1— ,<-)'^, (1 -x^fi , ... .
Then, by analogy, he obtained the expression for the general term in the
expansion of a binomial, and thus established the binomial theorem.
determine the value oi/{3-5).
X fix)
1 0-208460
Example 1. — From the following table of values of the function f{x)
A3
•29242
29029
28789
28523
0-237702
0-266731
0-295520
0-324043
-213
-240
-266
27
In this example
a=l, n = ^, w=\
/y =/(!) = 0 208460.
Ph
Hence /(3-5)=^ + i 5/, - I S'-^j
= 0-208460
14621
27-2
= 0-223106.
E.ramplc J.— Corresponding values of x and/(,t:) are given in the fuUowiiig
table :
4 0
4-1
4-2
4-3
4-4
4-5
/(-■)
54-5982
60-3403
66-6863
73-6998
81-4509
90-0171
Obtain the value of /(4-05).
9] FORMULAE OF INTERPOLATION 23
A certain amount of discretion is necessary in the application
of Newton's formula. Let us assume, for instance, that the
f!;iven observations are finite — say seven — in number, and
correspond to the values a - 3 w, a-2w, ... a + Ziv of the
argument as in the scheme, § 8. If a + n w lies between a - 3 w
and a-'lw, then Newton's formula fory(a + nw) will involve the
quantities /_:„ S/_ s, S-/_._;, S\f_ ^, 8^/_i, 8^/_ ^ and 8^/1, all of which
occur in the scheme. Thus a good approximation to the value of
/(a + nw) may be obtained by means of the formula. But if
a + 7iw lies between a + 2w and a + 3w, then in the Newton's formula
for /{a + n w) the only known terms are the first two depending
on /a and 5/g. In general, these are quite insufficient to determine
the value of f{a + n w) to the necessary degree of accuracy. It is
therefore expedient in this case to find a formula involving a
larger number of known quantities. A suitable one may be
obtained by replacing iv in the original formula by - w. Then,
since
we must replace t>j\ by
Similarly, since
/,-/-! or -S/_i.
we must replace S'-yj by
./; - i/Li +/L, or 671i .
In the same way we can determine the functions which must
take the places of S*./#, S^/j, .... Newton's formula then becomes
n(n-\) ^, , n(7i-l) hi-2) ^. ,
J{a-mv)=/,-ndJ^.^+ \^ ^ <^-./-i - -^ 37 ■^V-s+ •••
This expansion involves diiFerences along a line parallel to the
lower side in the scheme, and it may therefore be regarded as
a formula for Backward Interpolation. When a + n w lies between
a + 2 w and a + 3 1<;, the formula will depend on the quantities
fi-, 5/s, S-/^, S^/|, S^/i, 8f/^, and S^/o) all of which are known: it
is therefore very convenient when the value of fix) is required for
values of X near the end of the series of observations.
24 INTERPOLATION AND NUMERICAL INTEGRATION [ch. ii
Example 3. — From the following table of values oi f{x) determine the
values of /( 104 -25).
X
/(.f)
A
A--^
101
1030198
30906
102
1061104
31518
612
103
1092622
32136
618
104
1124758
32760
624
105 1157518
In this example a — 105, n — 0'75, "■ = 1
/y=/(105) = 1157518
/(104-25)
0-75 (0-75 - 1) 0-75 (0-75- 1) (0-75 - 2)
-/o-0-75o/.^ + ^^— ^5y_,- --^ 5y-«+...
=/o-0-75 5/'_, - 0-09375 5-/_, - 0-0390625 sy. + ...
= 1157518-0-75 (3-2760) -0-09375 (6-24) - 0-0390'625 (6)
= 1157518-24570
58-5
0-234375
= 1132889-265625.
Example 4- — Using the data of Example 2, find the value of/(4-45).
10. Stirling's. Formula of Interpolation.
As Newton's formula involves differences -which lie along a
slanting line it will frequently not be suitable for calculating the
value of the function for a value of the argument near the middle
of a series of observations. By a suitable transformation of this
formula, however, we may derive another involving differences
which lie on a horizontal line : such a formula is termed a Central-
Difference Formula.
Since
6J\ - 8/_^ = S^Jl
we have
Hence
8jl =ijlSj[, + }.8'-/o.
), lUj
FORMULAE OF INTERPOLATION
In the same way the higher differences which occur in Newton's
formula may be determined in terms of the differences which occur
on the same horizontal line Sis/,,. Substituting these in Newton's
formula we obtain
/{a + n w) - /; + 7i{ji 3 /; + }, 3;/;} + ^^^^^ — - { S-/o + /^ 5"_/o + -h ^Vo}
n {n - 1 ) {to - 2)
3"!
{fx oVo + § 8\/l + II S% + i S%} + . . .
=Jo + n II oy,3 + — 6-/o + -^^-^ — jx b-J, + ^ SVo
n {ii- - 1) (vr — 4j
i^ 8%
This is known as Stirling's Formula of Interpolation. It should
not be used for determining the value of a function for values of the
argument near the beginning or the end of a series of observations,
as it would not then be possible to obtain a sufficiently high order
of differences.
Example i.— Using the data of §9, Example 1, determine the value of
/(3-5).
.r
fU-)
A
A-
A3
1
0-208460
29242
2
0-237702
29029
-213
-27
3
0-266731
(28909)
28789
-240
(-27)
-26
4
0 "295520
(28656)
28523
-266
(-26)
- 26
5
0-3-24043
28-231
-292
6
0-35-2274
In this case
a
= 3, IV = 1,
11 -
--\
fo
= /(3) = 0-26673L
/(3-5)
= /o + i^5/o
+ i-
5Vo - tV 5Vo
= 0-266731
14455
2-
30
= 0 '281158.
Examj}/e ^^ — Using the data of § 9, Example 2, Hnd the value of /(4-21).
26 INTERPOLATION AND NUMERICAL INTEGRATION [cir. ii
If n is greater than 0'5, it is better to use backward differences.
The corresponding formula may be obtained by merely changing
the sign of w in Stirling's Forward Difference Formula, viz. :- —
/{a - n iv) = /; -nix of, + ^^ 8\tl - ^^ ^'^ ~ /x S\fl
71- (n- - 1 ) ^ , , 71 (71- - 1 ) («- - 4) .. .
4 !
Of course, either formula may be applied whether it is greater or
less than 0-5 : thus, when a check is necessar)^, it will be found
better to compute/(.r) by the formula not already used.
Example 5.— Use this Bickward Difference Formula to check Example \.
Here a = 4, c- = 1, n — ^
/o=/(4) = 0-295520
/(3-5) -/, - i M 5/o + I 5- A + -h 5Vo
= 0-295520 - 14328
33
2
= 0-281157.
The agreement with the former result is perfect when we take into
consideration the fact that the last figure is forced in both cases.
Example 4.— Using the data of § 9, Example 2, find the value of/(4-29).
11. Gauss' and Bessel's Formulae of Interpolation.
Another central difference formula may be obtained by trans-
forming Stirling's formula.
Since /^S/, = 8y,- |8V;
,xSi^j^ = S^f^-i8% etc.,
we see that Stirling's formula may be written in the form
J {a +71 IV) =j, + 71 8/^ + ^^, sv^ + ^ — Yr — - ^A
(n+l)7i{7l-\) (71-2)
^ {71 +2) {71+ I) 71 (71- I) {71 -2) ^^^. ^
5 ! • I ■■■
In this result, which is frequently known as Gauss Formula,
the /x-terms do not occur. The differences which occur in this
10, 11] FORMULAE OF INTERPOLATION 27
formula are those which lie on two adjacent horizontal lines of
the scheme in § 8.
Again, since t\ +/o ^ - /Vl
hav(
^f, = \^^f,-\^f, etc.
w€ nave
and hence
f\a + n w)
(^t + 2)(7^+ 1) 7i (« - 1) (?t - 2) ,.
+ g-^ 0/^+ ...
, ■. . 71 (lb -\) ^, 71 (n — \) (n - h) ^.. .
{n + l)n{n-l){n--2)
+ - ^ ^-^ — /^^A
+ ;-; 6'yj^ + . . .
This is frequently known as BesseVs Formula. In it the differ-
ences are those wdiich occur on a horizontal line mid-way between
the lines on which /J, and/j stand.
The formula most suitable for the construction of mathematical
tables, when differences of higher order than the fourth are
neglected, is obtained from the first of these on putting
When this substitution is made we get
./ (a + 71 iv) = /o + n 6/, + \^ ^ ' 6% \ ' ^J\
W (1 - 71-) (2 - 7l) ,, .
28 INTERPOLATION AND NUMERICAL INTEGRATION [cii. ii
Bessel's formula will be found particularly suitable when n is
equal to |-, for in this case the formula becomes
/(»+j„j.,4+ i^.ov,. (i±Mli^ll(izlVsv>...
in which only even differences occur.
As with Stirling's formula, none of these results should be
employed when the value of the function is required for values
of the argument near the beginning or the end of a series of
observations. Stirling's formula will be found more convenient
when 71 is very small, or when the number of observations given
is odd : Bessel's when n is nearly equal to a half, or the number
of observations even.
The backward difference formula corresponding to Bessel's
forward difference formula may be obtained in the same way as in
Newton's and Stii-ling's cases. But, except as regards sign, the
coefficients in the two formulae will be the same, and therefore the
one cannot be used as a check on the other. When a check is
required to a result computed from Bessel's formula, it is better
to compute it again using Stirling's formula.
Example /. —The values of a function for the values -2, -1, u, 1, 2, of
the argument are respectively 0-12569, 0-17882, 0-23004, 0-27974, 0-32823.
Using Gauss' formula, show that when the argument has the value 0-4, the
value of the function is 0-25008.
X fix) 1 X' A3
-2 0-1-2569
5313
- 1 0-17S82 - 191
5122 39
0 0-23004 -152
4^ ,^~ 31
1 0-27974 - 121
4849
2 0-328-23
In this example a = 0, n = O'i, ic = 1
/;, =/(0) = 0-23004.
Hence /(0-4) = /„ + 0 4 5/^ - 0-12 o'Vo 0-056 5^;..
= 0-23004
1988
18-2
= 0-25008.
11, 12] FORMULAE OF I?^TERPOLATION 29
Example 2. — Using the data of §9, Example 2, find the values of/(4"21)
and /(4 -29).
12. Evaluation of the Derivatives of the Function.
When, as often happens, the values of the differential co-
efficients of a function which is given by a series of observations
are required, use may be made of the symbolic equality established
in § 4, Cor. 4.
For, since
e"'~^f{:c) = Ef{x)
we have
w — .J\x) = log E .f{x).
The method will be evident from the following example : —
Examph 1. — The values of a function corresponding to the values I'Ol,
1-02, 1-03, 104, 1 05 of the argument are 1030301, 1-061208, 1092727,
ri24864, 1*157625 respectively. Find the value of the first differential
coefficient of the function when the argument has the value 1 "03.
Here ,r = O'Ol. Let /o = /(I -01).
d
Then *'• j^/(l-03) = log £". /(1'03)
= log E.f,
= log E.E^f,
= log(l+A).(l + A)2./„
/ A2 A^ X
= 1-^ - T+ T -...). (1 +2 A + A--'j./o.
Neglecting differences of higher order tlian the 4th, this gives
0-01-^^,/(l-03) = (A + |A2+iA^'-AA^)/„
= {j\-iE+^E'-j\E')/,
= i(/3-/i)-TV(/4-/«)
= 0-042437 -0-010610
= 0-0318-27
/'(l-03) = 3-1827.
Examplt 2. — The values of a function corresponding to the values 1-5, 2,
2-5, 3, 3-5, 4, of the argument are 4-625, 2-000, -0125, -1, 0-125, 4000,
respectively. Find the value of the first differential coefficient of the function
at the points where the argument has the values 2-5 and 3-0.
30 INTERPOLATION AND NUMERICAL INTEGRATION [en. ii
The method is perfectly general. The same result may be
obtained by diiierentiating any of the formulae of interpolation
already established. For example, if we differentiate with respect
to n the formula of Stirling,
f{a + n w) = /; + nix 8jl + J^ S\f„ + ^ '^~ /^ ^Vo
7r{n--l)
+ 4!"" -^"^ -
we have
wf {a + nw) = tJ. 8/0 + n 8' f, + ^{Sn--l)ix S\fl
+ ^\,{2rv'-n)S\f,+ ...
w-f" {a f n w) = S-_/q -\-nii S-'/q + J-.- (6 rr - 1 ) S-*/', + . . .
Putting ?i = 0 these become
wf (a) = IX 8/1 -}ix 8% + gV H- ^Vo - • • •
w"-/" (a) = S\/l - ^\ 8\fl + ^ 8% -...
These formulae are of use in determining the turning values of
a function, the values of the argument corresponding to given
values of the differential coefficients, the points of inflexion on a
curve of which isolated points are given, and also in determining
whether a series of observations is periodic.
Example 3. — The consecutive daily observations of a function being
0-099833, 0-208460, 0-314566, 0-416871, 0-514136, 0-605186, 0-688921,
0 '764329, show that the function is periodic and determine its period.
From the given observations we obtain the table : —
X
1
y=f{x)
0-099833
A
108627
A"
A"
A*
2
0-208460
106106
-2521
-1280
3
0-314566
102305
-3801
-1239
41
4
0-416871
97-265
-5040
-1175
64
5
0-514136
91050
- 6215
-1100
75
6
0-605186
83735
-7315
-1012
88
7
0-688921
75408
-8327
8
0-7643-29
12] FORMULAE OF INTERPOLATION 31
Taking a = 3, /« -/(a) =./"(3), we have, since tv =\,
/"(3) = svo- w5*/;-.
= -3801-iV41
= -3804
1 ^, 3804
•■• /(3) •/"(^) = - 314566
= - 0-0121.
In the same way we find that
(i) When a = 4, /, -/(4)
/"(4)= -5045, j;^/"(4)= -0-0121
(ii) When a = 5, /„ = /{5)
f"(o)= -6-221, jr^f"(b)= -0-0121
(iii) When a = 6, ./'o =./'(6)
/"(6)= -732-2, ^-:^-/"(6)= -0-0121.
1
Hence in all cases 77T\/"('t') - -0-0121. Putting y = J {x), this gives the
differential equation
1 (Py
- -r.,:= -0-0121
y dx-
d-y
whence -^., + 0-0121 3/ = 0
y - A cos 0*11 a- + 5sin011.T.
This result shows that y or f(x) is a periodic function of x, and that its
period is 2 tt/ Oil, or 57-12 days.
Example 4. — The following table represents one of the Bessel func-
tions, / (x) : —
1-0
0-765198
1-1
0-719622
1-2
0-671133
1-3
0-620086
1-4
0-566855
1-5
0-5118-28
1-6
0-455402
dy d-y
Find the values of -j-, and -j-,, for x — 1 -3, and hence, by substituting (with
four-place accuracy) in the equation
ah^ \ dy / 71^ \
d^^ + lc d^ + {^-lc^)y = ^'
determine the value of n.
INTERPOLATION AND NUMERICAL INTEGRATION [ch. ii
Forming the table of differences we have
Here
Hence
0-765198
0-7 19622
0-671133
0-620086
0-566855
0-511828
0-455402
A
- 45576
- 48489
-51047
(-52139)
- 53231
- 55027
- 56426
2913
■2558
2184
1796
1399
A^*
355
374
(381)
388
397
w = 0-1, rt = 1-3, /o =/(a) =/(l-3
0-l/'(l-3) = -0-052139 -i(0-0003Sl)
/'(1-3) = -0-52203
•01/" (1-3) = -0-002184 -iV (0-000014)
/"(I -3) = -0-2185.
Substituting these in
we have
t4 + - T^ + (1 - — )y = 0,
dx-^ X d X \ x^ } " '
•0-2185 + j-:^( -0-5220) + (l - Tf^.) (0-6201) = 0
-0-2185 - 0-4016 + (l - TfWi) (0-6-201) = 0
-0-6201 + (l - TyW,) (0-6201) = 0
1 + 1
(1-3)-^
ti = 0.
Hence the Bessel function in question is /(,(»').
Example 'j. — Show that the following observations are these of a periodic
function, and find its period : —
X (years)
,/■(.»■)
1
0 198669
2
0-295520
3
0-389418
4
0-479425
5
0-564642
6
0-644217
7
0-717356
8
0-783327
12] FORMULAE OF INTERPOLATION 33
Example tf. — Observations were taken of the cooling of a boiler, and the
results are given in the following table : -
t (time in hours) 1-2 19 26 3-3 40
e (temperature °F) 182-3 1752 ITl'l 1643 16U0
(le
Determine as carefully as you can the value of -^ wlien < = 2-6 hours.
MISCELLANEOUS EXAMPLES.
1. Find a polynomial of the third degree in x which assumes the values
37, II, L3, 29 when x is equal to 2, 1, -1. -2 respectively. Also find the
minimum value of the function.
2. A steam electric generator on three long trials, each with a different
point of cut off on steady load, is found to use the following amounts of steam
per hour for the following amounts of power : —
Lb. of steam per hour 40-20 66o0 10800
Indicated horse-power 210 480 706
Kilowatts produced 114 290 435
Find the indicated horsepower and the weight of steam used per hour when
330 kilowatts are being produced.
3. Using first differences only, find the value of /(4817'36) if
.<■ 4816 4817 4818
f{x) 69.3974063 694046108 694118146
To what extent may the result be relied on ?
4. Interpolate for .>; = 35-2967 from the following table and state the
accuracy of the result :^
X
/(-f)
35-28
162021-40139
35-29
162043-22198
35-30
16-2065-04777
5. If the sun's apparent longitude at O. M.N. on 1914 Aug. lat was
128° 2-2' 25"-8 ; on Aug. 3rd 130" 17' 15" -2 ; on Aug. 5th 13-2" 12' 7"-8 ; find
the longitude at G.M.N, on Aug. 2nd.
6. At the following draughts in sea water a particular vessel has the
following displacements : —
Draught in feet - - 15 12 9 6
Displacement in tons - '2098 1512 1018 536
What are the probable displacements when the draughts are 11 and 13 feet
respectively ?
7. The following table gives the expectation of life for males of the ages
mentioned : —
Present age .... 5 10 15 25 35
Expectation of life in years - 50-87 47-60 43-41 35*68 28-64
Find the expectations of 3 men whose ages are 18, 20, 30 years respectively.
34 INTERPOLATION AND NUMERICAL INTEGRATION [ch. ii
8. The HM table of the Institute of Actuaries gives the following values
of h :-
a:
l^
10
100000
12
99113
14
98496
16
97942
18
97245
20
96223
Obtain the values of Ix when x has the values 11, 13, 15, 17, 19.
9. The sun's apparent right ascension at mean noon is given in the follow-
ing table : —
June 1
4
34
0-69
6
4
54
31-75
11
5
15
10-64
16
5
35
55-34
21
5
56
43-18
26
6
17
30-99
Determine the apparent right ascension for 1914 June 15d. 16h.
10. Let the number of recruits in the United States army whose height
did not exceed x inches be denoted by y. Determine the number whose
height exceeded 66 inches, but did not exceed 66| inches, if
62
618
63
1855
64
3802
65
6821
66
10296
67
14350
68
17981
69
21114
70
23189
71
24674
m's apparent
declination 3 at mean noon ^
314, May 1
+ 14°
54' 19" -6
3
+ 15°
30' 25" -0
5
+ 16°
5' 28" -4
7
+ 16°
39' 27" -6
9
+ 17°
12' 20" -3
11
+ 17°
44' 4" -3
Determine the hourly difference in ? for 1914 May 5 and 7.
MISC. Exs.] FORMULAE OF INTERPOLATION 35
12. In the following table s denotes the space described at the time t by
a point in a mechanism. Tabulate the values of the velocity and the acceler-
ation, and obtain their values when t = OAo.
t (time in sees.) 0 O'l 0-:2 03 0-4 0-5 O'O O'T O'S 0-9 1-0 1-1
« (space in feet) 0 O-OS-' 0-12 0-2S1 0-J2j 0-804 1-200 1-725 2-229 2-753 3-2S2 3-Sll
13. There is a piece of mechanism whose weight is 200 lbs. The following
values of s in feet show the distance of its centre of gravity from some point
in its straight path at the time t seconds from some era of reckoning. Find
its acceleration at the time < = 2-05, and the force in pounds which is giving
this acceleration to it.
t 200 2-02 2-04 2-06 2-OS 2-10
s 0-3090 0-4931 0-6799 0-S701 10643 1-2631
14. The following values of p (lbs. per sq. foot) and d"C are for saturated
steam. Calculate with as much accuracy as the numbers will allow the value
dp
of ^for ^=115T.
e 105 110 115 120 125
p 2524 2994 3534 " 4152 4854
15. From the following table of the secular variations of the magnetic
declination determine the chief period : —
Year. Declination.
1540-5 - 7° -463
1580-5 - 8° -500
1620-5 - 6°-250
1660-5 - 0°-500
1700-5 + 8° -228
1740-5 +15° -755
1780-5 +20° -832
1820-5 +22° -380
1860-5 +19° -364
1900-5 +14° -833
[CH. Ill
CHAPTER III
THE CONSTRUCTION AND USE OF
MATHEMATICAL TABLES
13. Interpolation by First Differences.
All mathematical tables are alike in that they give numerical
values of the function for certain values of the argument, which
are so chosen that intermediate values of the function may
be derived by interpolation. The interval of the argument
varies in the different mathematical tables, its choice in
each case being determined by the rapidity of variation of the
function. In most of the published elementary tables, which are
intended for the use of students, the interval of the argument is
so chosen that the second difference is smaller than one unit in
the last place of decimals retained : on this account it is not
necessary to consider second differences in the interpolation, and
the function can be calculated for any value of the argument
intermediate between two of the tabulated values by a simple
proportion in the following manner : —
Consider a function f{x) which is tabulated for the values
.To, Xq + w, Xq + 2w,.... of the argument. Then, since the function
increases uniformly, when digits beyond a certain place of decimals
are neglected, we have
/(xq + w') -/(Xq) _ «;'
f{Xo + w)-f{Xa) w
where 0<w'<w.
From this it follows that
in'
/{x, + iv') =/{x,) + — {/{X, + W) -f{Xo) },
so that, if the first difference /(a-Q + w) -/(x'o) is known, the value
of the function for any value of x between Xq and Xq + w may be
found. This equation is known as the " Rule of Proportional
Parts." Geometrically, what we do is to I'eplace the curve through
two consecutive points by the chord joining them.
13] MATHEMATICAL TABLES 37
If the first differences are only sensibly constant, the result
obtained will be subject to a certain error. If tv = n iv where
0</i< 1, then the error may be put in the approximate form
n{n-\) ,
wj- (x-^rn w)
where 0<n'<l.
Example, 1. — In the following table the entry is the cube root of the
corresponding argument. Find the cube root of 1619"25.
Argument. p]ntry. A
1619 11-7421858
24171
1620 11-7446029
24161
1621 11-7470190
If we desire accuracy up to the fifth decimal place, then we may regard the
first differences as constant.
Here ,r„^1619 »'-=l !/y = 0-25.
Hence ^^1619-25 = 11 74219 + -^^ (000-242)
= 11-74280.
Example 2. — If a; = 1000 show that, when second differences are
neglected, the error committed in computing by this method the value of
log 10(1000 +/'•), where 0 -= ir < 1, will have no effect on the sixth decimal
place of the logarithm.
Here log {x + w) = log x + w { log {x + 1 ) - lug x }
and the error committed is
E-~^''^\-;-^/"ix+.')
where
0<m;'< 1 and /(.r) = logio.i-.
Now
J [x+w) = ^2 logjo(.c + M;')
-1 "' t
I (x + iv')' I
where m — modulus of common logarithms = 0-43429.
But the maximum value of w{l - »•) is ^ since
Hence the error committed is less than
J 0 43429
' ■ (1000+H-')- '
1
16 . 10«
0-0000000625,
38 INTERPOLATION AND NUMERICAL INTEGRATION [cii. iii
which shows that the error committed does not affect the sixth decimal place
in the logarithm of the number (1000+ id).
Example 3. — The differences between the six-figure logarithms, to base 10,
of the successive integers between 5000 and 5100 are sensibly constant.
Prove this and determine the value of the constant difference.
14. Case of Irregular Differences.
When the second diiFeiences cannot be neglected, however,
this method cannot be applied. This happens frequently when
the rate of change of the function is large compared with that of
the argument, as is the case with the logarithmic sines and
tangents of small angles. We must then use the formulae of
Chapter II., including differences to the order required.
Example 7.— To obtain log sin 0° 16' 24" -5.
From Schron's tables we have
A A- A»
iysinO°16'10" = 7-6723450
44543
L sin 0° 16' 20" -7 -6767993 -452
44091 9
L sin 0^^ 16' 30" = 7 6812084 - 443
43648
i/sinO 16'40" = 7-6855732
Using Gauss' formula
n(n-\) (n+\)n{,i-\)
f[a+nii-)=^j„ + nbjy+ — ^y^ 5-./o + 5T~ ^".4
we obtain
L sin0° 16' 24" -5 = 7 -6767993 - 1
19841
56
= 7-6787889
This is the answer correct to 7 places, the value correct to 10 places being
found from Vega to be 7-6787889383.
The value obtained by the " Rule of I'roportional Parts" is 7-6787834.
This result is correct only to 4 places of decimals, and shows the necessity of
interpolating by some such method as the above the value of the logarithmic
sine of a small angle. This may be shown to be also true in the case of the
logarithmic tangents of small angles.
Example, i?. — Compute the value of log tan 0° 15' 35".
14-16] MATHEMATICAL TABLES 39
15. Construction of New Tables.
In spite of the large number of different mathematical tables
that have been constructed to assist the computer, further
research in the various branches of mathematics and the kindred
sciences shows that these are not sufficient for all purposes. A
research student often requires new tables for the question which
he is discussing, and these he must construct for himself.
The first question to be settled is the degree of accuracy to
which the tables are to be computed, i.e., the number of decimal
places to be retained. This varies in different cases ; but it is
always desirable to extend the original calculations beyond the
number of places which are ultimately to be retained, in order to
avoid the cumulative effect of the neglected digits.
The general plan of the calculation is to compute first with great
accuracy the values of the function for a series of widely-spaced
values of the argument. For this purpose infinite series are generally
employed. The values of the function for a more closely-spaced
series of values of the argument are then derived from these by
means of the formulae of Chapter II. : this second process is called
Subtahulation. We shall consider these processes in §§ 17, 18.
16. Tabulation of a Polynomial.
The simplest case to be considered is the tabulation of a polynomial.
If the degree of f{x) is n, it will not be necessary to compute its value
directly for more than (w-l- 1) values of the variable, as may be seen from the
following examples.
Example 1. —Let J{x) = .c^ + 3.v- + 2x-\.
When .I'-O, 1, 2, 3, we have f{x)= -1, 5, 23, 59. We thus obtain the
following table : —
Ax) A A-i
-1
12
18
23
36
3 59
Since /(a) is of the third degree in x, its third differences will be constant ;
in this case they will be equal to 6. Extending the column A'^ by inserting
6'8 we can then extend the column A- by simple additions or subtractions.
A repetition of the same process will enable us to extend the column A, and
40 INTERPOLATION AND NUMERICAL INTEGRATION [CH. in
thereafter to determine the values of /{x) for all integral values of ;r. The
resulting table is
X /(x) A A^ A3
-3 -7 -12
12
23
36
60 6
4 119 30
90 6
5 209 36
By this simple method we are enabled to determine the values of a poly-
nomial f{x) for all positive and negative integral values of .r. We can then
use any of the formulae of interpolation already established to obtain the
values of f(x) for fractional values of .r. The method may be extended to
functions which are not rational or integral, but in this case it will be
necessary to make the interval between successive values of x so small that
the dififarenoei of some definite order are constant throughout the range of
values under consideration.
Example 2. — Form the table of values of J'{x)-x*^ ox^ - T.c-^-.c - 4 for the
values 0, 1, 2, 3, 4 of the argument ,r, and then deduce the values of /(,()
corresponding to the values - 1, - 2, - 3, - 4 of ,r.
17. Calculation of the Fundamental Values for a
Table of Logarithms.
We shall now illustrate the process of computing tables of
functions other than jaolynoniials by describing the calculation
of a table of logarithms.
The first stage in this consists in determining the logarithms
of all the prime numbers between certain limits. It is evident
that when the logarithms of the primes are known, the logarithms
of all composite numbers can be derived from theui by mere use of
the formula
log a h = log a + log h.
17J MATHEMATICAL TABLES 41
[n order to obtain the logarithms of the primes we take two
equations whose roots are integers and which differ only in their
constant term. Such a pair are
X" - 3 u; + 2 = 0
ic- - 3 a; - 2 = 0.
Now x'-3x + 2 = (x-iy (x + 2)
x-^-Sx-2 = {x+l f {x - 2).
Hence
log {x-\) + h log {X + 2)- log {x+l)-h log {X - 2)
_ a''-3a;+2
2^^o
x" -
1 +
-3
X-
o
■J
0
ilog
OCT'
-3
X
1 -
o
3x
which, by use of the ordinary logarithmic series, may be written,
if the logarithms are all to the base e,
2 I 2 1^ r 2 V'
X' - 3 X '^ Ix''^ - 3 x) ■' Ur -3x1
This is called Bordas Formula.
Putting x = b, 6, 7, 8, we obtain the four equations
log 2 - f log 3 + Mog 7 = 0-0181838220854376...
4- log 2 + logo- log 7 = 0-0101013536587597...
- 2 log 2 + 2 log 3 - i log 5 =0 0062 11 2599992784...
- ^ log 3 + i log 5 + log 7 = 0-0040983836020895...
Solving these four equations we obtain the logarithms (to the
base e) of the prime numbers 2, 3, 5, 7, viz. : —
log 2 = 0-6931471805599453...
log 3 = 1-0986122886681096...
log5 = 1-6094379124341003...
log 7 = 1-9459101490553133...
Further substitution for. x in Borda's formula will give immedi-
ately the logarithms of the other prime numbers. For exauiple, on
putting x = 9 the left-hand side becomes
\o" 8 + .1, lo" 1 1 - loir 10 - J log 7.
42 INTERPOLATION AND NUMERICAL INTEGRATION [ch. in
Of these log 7 has been found ;
logs -3 log 2 = 2-0794415416798359... ;
and log 1 0 = log 2 + log 5 = 2-3025850929940456 ....
Hence we obtain log 11 = 2-3978952727983705....
It is to be noted that these logarithms are to the base e. But
as we have now found log^ 10, we can at once derive its reciprocal
logio 6 = 0-43429448 190325 18..., and all the above logarithms
can now be converted to the base 10 by multiplying them by this
"modulus."
We thus obtain the table
logio 2 = 0-3010299956639812
logjo 3 = 0-4771212547196624
logio 5 = 0-6989700043360188
log,o 7 = 0-8450980400142568
logioll = 1-0413926851582250
etc.
For the calculation of the logarithms of the larger primes Haros" Formula
will be found more useful. In this case consider the equations
.f*-25.i-- = 0
.c-"- 25 .1-2 + 144 = 0.
From these we deduce
(■f + 5)(.u-5).r- __ x^-2ox^ _ (x-^ - 25 .r'-' + 72) - 72
^^^ (.f + 3) (.1- - 3) (.f + 4) {x - 4) ~ ^°S x^ - 25 .x-'- + 144 " ^°2 (x^ - 25 x^ + 72) + 72
whence, if the logarithms are to the base 10, we have
log(.r + 5) = {log(.i- + 3) + log(.c-3) + log(.c + 4) + log(,v-4)}-{log(.c--5) + 21og.i'}
^, r ?2 /__ 72 x3 \
'-''"\.T^-25.r= + 72 + ■^Kx^-lbx'^l-l) +•• /
where m = log^j e.
Example 1. — To obtain the logarithm of the prime number 43 to the
base 10.
Let .v + 5 = 43.
Then x = 3S
log 43 = log 41 + log 35 + log 42 + log 34 - lug 33 - 2 log 38
- 0-8685889638065036 [0-000035137'24O'2040 + 0-UUOO000O00OO0145[.
But the logarithms of 41, 35, 42, 34, 33, 38 are assumed to have been
previously calculated. Hence
log 43 - 1 -6127838567197355 - 1 -5185139398778875
1 -544068044350-2756 3- 159567 1932336'203
I -6232492903979005 00000305198190724
1-531478917042-2551
-1-6334684555795865
Example ^.— Show that logjo l;31=i 2 -1172712956557643.
17-18] MATHEMATICAL TABLES 43
18. Amplification of the Table by Subtabulation.
Having now obtained the logarithms of the prime numbers
within a certain range, we derive from them, by mere addition, the
logarithms of the composite numbers which lie in this range, and
thus obtain the logarithms of all the integers within the range.
We shall now show how, by subtabulation, the extent of this table
may be increased manyfold. It will be supposed that 7-place
accuracy is required.
For example, suppose the logarithms of the integers 31, 32, 33,
34, 35, 36 have been obtained directly ; we shall show how the
logarithms of all the integers between 330 and 340 can be derived
from them.
Neglecting the characteristics of the logarithms, which can
always be written down by inspection, we form the following
table : —
X-
log X
A
A-
A^
A*
31
•49136169
1378829
32
•50514998
1336396
-42433
2535
33
•51851394
1296498
-39898
2312 (■
-223
-208)
34
•53147892
1258912
- 37586
2120
-192
35
•54406804
1223446
- 35466
36
•55630250
Putting a = 33, n=^
iV '"=^ "1
Gauss' form
ula
/{a + u IV) =y;
+
+
, . «(1-
'' "-^i 2 !
n (1 - ,v) (2
4!
(1-
3
r'^'A
we
obtain
/(33-1) =/; + 0^1 6/, - 0-045 6-/o - ^-0165 6;/^ + 0^0081 /x 87;
= •51851394-38
129650 1
1795
= •51982800.
44
INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
Thus to 7 decimal places log 331 = 2-5198280. But it is not
necessary to calculate the values for all the numbers between
330 and 340 in this way. It is, in fact, simpler to calculate
log 332 from the formula
y (a + ^% w) =/ {a + T-V ic) +[J{a + j\ w) -/(a) } + 8.,
where 8.,=/{a+ f^w) - '2/{a + -^w)+f{a)
= 0-01 87o + 0-001 sy, - 0-001 At sy,
and then to calculate log 333 from the formula
J{a^ j;"^ w) =f(a + ^oW) + {/{a + ^o «') -/(« + to "^) } + ^3
where S3 =/(« + f '„ w) - 2f{a + ^^ 10) +/ {a + -iy ^)
= 001 SVo + 0-002 8V^ - 00017 {x SV, .
The complete table of these quantities S is as follows : —
8„ = 0-01 SVo + 0-001 8\f, - 0-0010 /x 8'J]
S3 = 0-01 Sy + 0-002 Sy - 0-0017 /j. 8\f^
8, = 0-01 s-/; + 0-003 8' /; - 0-0019 /x sy;
Sg = 0-01 6-y; + 0-004 sv; - 0-0021 /x sy^
s, = 0-01 sy + 0-005 sv; - 0-0020 iM sy
8, = 0-01 Sy + 0-006 S-;/' - 0 002 1 /x 8\/[
Sg = 0 -01 sv; + 0 -00 7 s v; - o -oo 1 8 /x sv;
Sg - 0-01 sv; + 0-008 sv; - 0-0017 /x sv;
8j„ = 0 -0 1 sVo + 0 -009 svi - 0 -00 1 2 /x sv';
The numerical work is arranged as in the following table, in which
the logarithms required are found in the last column.
-39898-0 + 2
312-0-2
08-0
129649-8
2-51851394
1
To
+ 1795-4-
38-1 -
1-7 =
- 1755-6
131405-4
2-519827994
To
- 399-0 +
2-3 +
0-2==
- 396-5
131008-9
2-521138083
fV
- 399-0 +
46 +
0-4 =
- 394-0
130614-9
2-522444232
A
- 399-0 +
6-9 +
0-4 =
- 391-7
130223-2
2523746464
5
TO
- 399-0 +
9-2 +
0-4 =
- 389-4
129833-8
2-525044802
To
- 399-0 +
11-6 +
0-4 =
- 387-0
129446-8
2-526339270
To
- 399-0 +
13-9 +
0-4 =
- 384-7
129062-1
2-527629891
To
- 399-0 +
16-2 +
0-4 =
- 382-4
128679-7
2-528916688
10
- 399-0 +
18-5 +
0-4 =
- 3801
128299-6
2-530199684
To
- 399-0 +
20-8 +
0-2 =
- 378-0
127921-6
2-531478900
18] MATHEMATICAL TABLES 45
The method followed in constructing this table will best be shown by
considering two examples.
In the line tV
+ 1795 '4 is the numerical value obtained for 0'045 5Vo
- 38-1 „ ,, „ -0-0165 8%^
- 1-7 ,, „ „ 0-008 m5^/i
1755 "6 is the sum of the three preceding numbers.
129649'8 is the numerical value obtained for 0 01 dj\.
131405-4 is the sum of the two preceding numbers, and therefore
represents
0 -01 d/^ + 0 -045 5-Jl - 0 0165 8\/\ + O'OOS /x d\f^ .
2-51851394 is the numerical value of /;,
2-519S27994, which is obtained as the sum of the two preceding
numbers, represents
/o + O-Ol 5./; + 0045 5Vo- 0-0165 o V) + 0 -008 m 5 Vl •
In the line ^j^
- 399-0 is the numerical value obtained for 0-01 S"/o
+ 2-3 „ „ „ 0-001 5V;
+ 0-2 ,, ,, „ -0-00lfx5{f\
-396-5, which is obtained by summing the three preceding numbers,
represents
0-01 8"fo + 0-001 5 VI - 0-001 M 5V\ or 5.,.
131405-4, which has previously been obtained, represents
131008-9, which is obtained by summing the numbers immediately
above it and to the left of it, represents
/(« + tV "•)-/(«) + 5,.-
Now 2-519827794, as we have already seen, represents f(a+i\w).
2-521138083, which is obtained by adding the numbers immediately
above it and to the left of it, represents
f(a + r^w).
The results in the last column are the logarithms of the numbers
330, 331, ... 340. They cannot, of course, be relied on up to the
last figure, but will be found correct up to seven decimal places.
This accuracy may be obtained by taking one decimal more in the
first difference than in the tabular function, one decimal more in
the second difference than in the first, and so on, when these can
be obtained. A check on the accuracy is obtained by carrying
on the process until log 340 is computed. Comparison of this
result with the value for log 34 will show whether the computa-
tion has been accurately performed.
46 INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
We may remark that it is an improvement on the foregoing-
process if the differences required for the subtabulation are
computed independently from the expressions in series, instead of
being derived by forming a table of differences.
Example i.— Given log 24 = 1 -38021 12
log 25 = 1 -3979400
log 26 = 1 •4149733
log 27 = 1 4313638
log 28 = 1-4471580
log 29 = 1-4623980
obtain the logarithms of the numbers 261, 262, ... 269.
19. Radix Method of Calculating Log^arithms.
When the logarithm of a number is required to a large number
of decimal places, its determination either by means of an infinite
series or by interpolation becomes very laborious and liable to
error. A much simpler and more useful method in such a case is
that known as the Radix Method.
Let N be the given number whose logarithm is required to n
places of decimals. Then since
ir=10^.a.{l+iVo)
where N^ is a decimal, x an integer, and a = 0, 1, 2, ... 9,
logio N=^x + logjo a + logio (1 + i\'o)-
Thus, to determine logjoiV, it will be sufficient to compute log,oa
and logio(l + N^). The former has been obtained above. We shall
now show how to evaluate the latter.
Let r■^ denote the decimal consisting of the first m figures in
iVo) ^"d let A^i be the difference of Nf^ and r^ Then
Now divide iVj by (1+?-,) until there are m figures in the
quotient r„. If we call the remainder N.j, then
N,-{l+r,)ro = N„.
Similarly, if the remainder, after dividing No b}' (l+ri)(l+n)
until there are m figures in the quotient r., is denoted by N^, we
have
N,-{l+r,){\+r.;)r, = N,,
and in general,
iV^_l-(l+r,)(l+r,)...(l+r^_Or^ = ^"^.
19] MATHEMATICAL TABLES 47
These equations may be written in the form
JV, -JV. ={l+r,){l+r,)-{l+r,)
JV. -N, ={l+r,){l+r^(l+r,)-{\+r,){l+r,)
N,_,-N^_, = {\ +n) (1 +r,) ... (1 +.v_0
-(l+r,)(]+n)...(l+7V,)
iVV:-i\^^ =(l+r,)(l+r,)...(l+rv)
-(l+r,)(l+r,)...(l+r^_i)-
Adding, we obtain
(1+iV^o) -A\ =(l+rO(l+r,)...(l+r,)
whence
1 + i^o = (1 + n) (1 + r„) . . . (1 + ?v) + A^j,
If, therefore, the remainder A^^, has no significant figures in its
first n decimal places,
log(l+iV'o) =log(l+r,) + Iog(l+r,)+ ... +]og(l+r^,)-
The quantities (l+r,), (1+?^), ... are called the radices. As
they will be of the form [ 1 + Ynf) where k and I are integers, and
as k is always small compared with 10' we may compute log (1 + r^)
by means of the logarithmic expansion
log {^+r,) = r„-hrj' + \rj; - { r^/+ ...
The values of the logarithms of these i-adices when k = 0, 1, 2, ...
999, and ^ = 3, 6, 9, 12, 15, 18, 21, 24, as well as of the numbers
1, 2, ... 9, have been computed to twenty-four places of decimals.*
The inverse process of determining the antilogarithm corre-
sponding to a given logarithm is equally simple. When the
logarithm is not in the table, we have merely to take out the
next lower and subtract it from the given one. Then take out
the next lower and subtract it from this remainder; and so on
until the remainder consists of zeros. The product of the radices
corresponding to the logarithms taken out is the number required.
* Gray, Tables for the formation of Logarithms and Antilogarithms to
twenty-four or any less number of places. London, C. & E. Layton (1900).
48 INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
We shall illustrate these two methods in the following
examples : —
Example. 1. — To determine logj||43 to 24 places of decimals.
Here.r = l, a = 3.
Taking m = 3, as in Gray's tables, we have
\ + N^ = \- 433 333 333 333 333 333 333 333
l+r, = l-433
A'l = 0 • 000 333 333 333 333 333 333 333
r.3 = 0- 000232
A'2 = 0- 000 000 877 333 333 333 333 333
( 1 + ri ) r., - 0 ■ 000 332 456 /. e. , 0 • 000 333 333 - 0 • 000 000 877
{\+r^)(\+r.^ = [\+r,) + (\+r,)r.,
= 1-433 332 456.
The above will be sufficient to indicate the general theory. The process
can, however, be simplified by arranging the work in the following manner : —
1 • 433 333 333 333 333 333 333 333
1-433
DifiF. ^ 0- 000 333 333 333 333 333 333 333
1 • 433 ) 0 • 000 333 333 333 333 333 333 333 ( 0 • 000 232
286 6
46 73
4299
3743
2 866
1-433333333
0-000000877
Remainder.
Diff. = 1-433 332 456 = New Divisor.
1 • 433 332 456 ) 0- 000 000 877 333 333 333 333 333 ( 0- (OOOf 612 '
859 999 473 6
17 333 859 73
14 333 324 56
3 000535173
2 866 664 912
133 870261 = Remainder.
1-433333 3.33 333 333333
0-000 000 000133 870 261
Diff. = 1 - 433 333 333 199 463 072 = New Divisor.
Should the original dividend not be known beyond the 24th place of
decimals, we should now have to proceed by a method of contracted division.
In the present ease we know the figures after the 24th decimal, and may
therefore continue as above. But we shall adopt the contracted method.
* 0-(000)-612 is written in place of 0- 000 OCO 612
0- (000)^093 ,, „ 0-000000000 093 and soon.
I
19] MATHEMATICAL TABLES 49
The new divisor contains 19 figures ; the new dividend
•(000)3133 870261333 333
contains only 15 significant figures after the point. We must therefore cut
off 4 figures from the divisor. Thus we obtain
1 • 433 333 333 199,4^6 ) 0" (OOO)^ 133 870 'J61 333 333 ( 0- (000)3 093 *
128 999 999 987 951
4870 261345 382
4299 999 999 598
570 261 345 784 = Remainder.
1-433 333 333 333 333 333 333 333
0 • 000 000 000 000 570 26 1 345 784
Diff. = 1-433 333 333 332 763071 987 549 = New Divisor.
As the new dividend 0- (OOO)"* 570 261 345 784 contains only 12 significant
figures and the divisor "25 we must cut off 13 figures from the latter. It
then becomes 1-433 333 333 33. As these are the figures in the original
dividend there will be no further new divisors. This will usually happen
when the fourth division is reached. We may therefore determine the
remaining figures by continued contracted division.
1 • 433 333 333 33 ) 0 ■ (000)^ 570 261 345 784 ( 0 - (000)^ 397 856 752 873
430000000000
140261345 784
129 000000000
11261345 784
10033 333333
12-28 012451
1146 666 666
81345 785
71666 667
9679118
8 600 000
1079118
1003333
75785
71667
4118
2867
1251
1146
105
100
~5
4
50 INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
The remaining radices are thus
1- (000)^ 397
l-(000pS56
1 • {000)« 752
1 • (000)7 873
Using Gray's tables of the radices we have
logs =0-477
121254 719 662437 295 028
log 1 • 433 = 0 • 1 56 246 1 90 397 344 475 994 693
log 1-000232 =
100 744 633875 845 680796
log I -(000)- 612 =
265 788141593627 087
log 1- (000)3 093 =
40 389386 815124
log 1- (000)^397 =
172414 909316
log 1- (000)5 856 =
371756077
log 1- (000)6 752 =
326589
log 1- (000)' 873 =
379
Sum = 0-633468 4555795865-26 405089
Hence since .t'= 1 the required logarithm is
log 43 = 1 • 633 468 455 579 586 526 405 089.
Example 2. — To determine the antilogarithm of
1-120 619 882 724133 396 646 450.
As the integral part only affects the characteristic we shall neglect it for
the present.
Using Gray's tables we then obtain
0- 120 619 882 724 133 396 646 450
log 1-3-20 =0- 120 573 931 205 849 868 472 706
Diff.
log 1-000105
Diff.
log 1 -(000)2 812
Diff.
log 1 -(000)3 793
Diff
log 1 -(000)^ 443
Diff.
log 1 -(000)5 922
Diff.
log 1 -(000)6 831
Diff.
log 1 -(000)' 516
Diff
45 951518 283 528173 744
45 598 526 719080137 349
352 991564448 036 395
352 646976130 787551
344 588 317 248 844
344 395 5-24 012 7-26
192 793 236118
192 392 455 483
400 780635
400419512
361 123
360899
224
Multiplying these radices together we have
1 -(000)^ 443 X 1 -(000)' 922 x 1 -(OOO)* 831 x 1 -(000)' 51(
1 (000)^443 922 831 516.
19. 20] MATHEMATICAL TABLES 51
This product of the last four factors can always be written down by
inspection. The other multiplications may be arranged as follows : —
1 • 000 000 000 000 443 922 831 516
1-000000000 793
1 • 000 000 000 000 443 922 831 516
700 000000 000 310
90000000 000040
3 000000 000 000
1 • 000 000 000 793 443 922 831 866
1-000000 812
1 • 000 000 000 793 443 922 831 866
800 000000 634 755138
10 COO 000 007 934 439
2000000001586 888
1 • 000 000 812 793 444 567 108 331
1- 000 105
I - 000 000 812 793 444 567 108 331
100000 081279 344 456 711
5 000004 063 967 2-22 836
1 • 000 105 812 878 787 878 787 878
1-320
1 • 000 105 812 878 787 878 787 878
300 031 743 863 636 363 636 363
20 002 1 16 257 575 757 575 757
1 • 3-20 139 672 999 999 999 999 998
Since the cliaracteristic in the given logarithm is 1, there will be 2 figures
before the point in the corresponding antilogarithm. Thus the required
result is
13- -201 396 729 999 999 999 999 98.
Exami^le 5.— Given tt^:^ 3- 141 592 653 589 793238 462 643 show that
log TT = 0- 497 149 872 694 133 854 351 268.
Examjih 4- — Show that the antilogarithm of
0-434294 481903-2518-27 651129 is 2- 718 "281 828 459045 235 360290.
20. Inverse Interpolation.
Suppose the values of a function have been tabulated, and that
it is required to find the value of the argument which corresponds
to some given value of the function, intermediate between two of
the tabulated values. The process by which we obtain this value,
which may be called the "antif unction," is known as Inverse
Interpolation.
52 INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
Lety (a + n w) be the given value of the function. Then
^ . n in - 1 ) ^, .
J{a-^nw) =y,j + n 6j, + - ^ , - - 5;/j + . . .
7(,(n - I) ... (71 - 7- + I) , .
+ ^ — , oy,. +...
?• . Y
In this equation /(a + 71 w) is given, and as the vakies of
fotfi, ... have been tabulated, the values of of^, 8"-/], ... are also
known. Hence the only unknown quantity in the equation is 71,
and this is the numbei- required for the determination of the anti-
function (a + 7iw).
To obtain it we must solve this equation. It has been pointed
out already that the differences of the function under consideration
usually become sensibly constant at some definite order, say r.
The higher differences will therefore be zero, and this equation then
becomes one of degree r in ?t. It can, therefore, be solved by a
method of successive approximation, as follows : —
If we write the above equation in the form
V^^~^-A^
and neglect differences of the second and higher orders, we get as
a first approximation the value «i where
Let us next take second differences into account and put jIj for n
in the denominator of the expression on the right-hand side. Then
the next approximation 71.2 is given by
f—f
When third differences are taken into account, 7i.^ is substituted in
the denominator for w^, giving
^ .f^Ii
""" sy; -f 1 (n, - 1) sv; + \ {71., - 1) {n., - 2) 37: •
Proceeding in this way we obtain closer and closer approxima-
tions to the value of 71. The method, though rather laborious, has
the advantage that any error introduced in the computation of any
20] MATHEMATICAL TABLES 53
one appi-oxiniation will not cause the final result to be wrong,
being itself rectified in the higher approximations. Thus the sole
effect of the eri-or introduced is merely to make the approximation
less rapid, and therefore to increase the amount of labour involved
in the accurate determination of n. On this account the process is
a safe one to employ.
Other formulae for determining the antifunction may be
obtained by reverting any of the formulae of interpolation. We
shall exemplify this by reverting Stirling's formula
/(a + 71 tv) =Jl + „ ^ 8/o + — SV; + 3, H- SYn
?i-(w--l),, . 71(71" - I) (n- - 4:1 .. .
Neglecting differences of higher oi-der than the fifth, and arranging
in powers of «., we have
/(a + 71 w) =/„ + (/. 5y; - 1 /^ 5^0 + .V /^ sv;) ,i + (i sv; - a. sy^) ,,^
which may be written
F= An + Bri" + Cn^ + D71' + E7i^
where
F=f{a + nw)-fQ
A=ix8/l-lix?i'fo + rhl^^%
C = hl-^Vo-^\f^8Vo ,
^ = ASVo
^=Tk/-S7;
Reverting
this series we obtain successively
(i)
F
(ii)
F B /F\-
F B /FV- f,/BV C\fF\'
54 INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
F BVF B (F\^ fo^^V ^l/^l'T
" Z IT ^ T u / / ~ vi U /
" ^ ^ U/ y^\A) a]\a]
A more convenient form * may be obtained by putting
F n
— = Wj and -7- = ^' a-nd rearranging the terms.
We then obtain the formvila
n = Wi + /i 7ii +/^ ?^l- +/; ?z/ 4-/4 Wj''
where
/ = _ 5r + 2 {Brf - 5 (^r)^^ + 14 (Br)'
f,^Cr{ -I +5 Br- 21 {Brf}
f, = Dr{-l + 6Br) + 3{Crf
J\=-Er.
If the first derivate of the function /{a + 7iiv) is tabulated, the
values of the quantities A, B, C, D, E may be obtained by forming
the successive differences of wf {a + n tv). As fourth difFei-ences of
the derivate are of the same order as fifth differences of the
function we need not extend the table beyond this order of
differences.
* Van Orstrand, Phil. Mag. (6) 15 (1908), p. 630.
20] MATHEMATICAL TABLES 55
Argument.
Entry.
A
A-
A3
A^
a - 3w
«'./'-3
""^f'-n
a -'2 IV
wf'-2
wSf _^
n-B\r_,
w^/L.
a-w
n-f'-.
wS/'_,^
«.'S7"-.
w^/:^
w8Y'-.
a
"^./■'o
ic\f\^
wS^'o
w^'^
wS'f'o
a + w
Wf\
10 S/\
w o-f\
«;5y'3
wB^f\
a + 2 9v
Wf'r,
^v ¥\
w8\r.
a + 3w w/'s
Differentiating the equation
/{a + nw) =y; + An + Btr + Cnf + Dn* + Eiv"
we see that
A=wf\: B = \w\f:'; C = \vff;"; D=.}fW'f^'^
Now from Stirling's formula we have
n'" j,„ ., n hr - 1 ) ^ .,
wj" {a + nw) = IV J 'o + nwfxoj Q+ ^ iv &f ^ + ^-y— w fi 6' J „
4!
oo ., 3 ?i" - 1 „„ ,, 2 n" - 71
w-f ' {a + n w) = n-fx8j o + n?r6y „ + ^ tv ixd'J „ + — p^ — wdj q
iv\f"' (a + n 71-) = 7vS\f'o + mvixS'f\ + '~^r^ w^Y'o
?6-V<-" (a + n tu) = IV fjL Sy ',j + n tv 8\f'o
?.t'5/*^' {a + n ^v) = w Sy 0
Putting w = 0 in these equations we obtain
A = IV f (a)
5 = 1 10' f" (a) = i XV IX. 8/'o - yV w /x Sy'o
c = i ioy '" (a) = 1 .^ sy 0 - ^v xv sy 'o
^=Tk^^y'''(«)=Tk^5y'o.
56 INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
A third method of determining the antifunction is to determine
?ij as in the first method, and then, using shorter intervals, to
retabulate the function for values of n in the neighbourhood of n^
The second approximation is now obtained from this new table in
the same way as n■^ was obtained from the original table. One or
two repetitions of this process will give the value of n to a very
high degree of accuracy.
An approximation to the value of the error committed by computing ?i
by the first method will next be obtained.
Since
/-/o
n-1
f-fo
n-\
^f.+^jW'f" [x + n'w)
■-/o f, , »i - 1 o f"{x + n'w) \ -1
1^1
the numerical value of the error committed is
54 • 2! "■' 54 I
I 7i(?i- 1) _^ /" [x^-n! w) I
Example 1. — The following table gives the logarithms of the distance
of Venus from the earth. To find when the logarithm had the value
9-9351799.
Date. log. A A^ A^ A* A^
1914 Aug. 18 9-9617724
- 74079
20 9-9543645 - 1497
- 75576 - 46
22 9-9468069 -1543 3
-77119 -43 1
24 99.390950 (-77912) -1586 (-41) 4 (1)
- 78705 - 39 0
26 9-9312245 -1625 4
- 80330 - 35
28 9-9231915 -1660
- 81990
30 9-9149925
We shall use the second method described above.
20] MATHEMATICAL TABLES 57
Let
/•o= 9-9390950.
Then
'a= -0-0077912 + 0-0000007= -0-0077905
if =-0-0000793
C= -0-0000007
D= 0
E= 0
9-9.3ql799- 9-9390950 _^,_. -
"i- -0-0077905 -0-50253O
r= -64-5062
5r= 0 0051153
Cr= 0-0000452
2)r= 0
Er^ 0
/i=- 0-0050637
f\n^ = -0-002545
/;= - 0-0000440
fyi^^= -0-000011
/,= 0
Hi= 0-502535
./•4= 0
n = 0-499979.
Hence the time required is
1914 Aug. (24 d. + 0-499979 x 48h.)
or 1914 Aug. 24 d. 23 h.
59 niin. 56 s.
Example S.—To show that when .t = ]
1000 and w -= 1 the error in com-
puting w from the formula
log (a; + 10)
- log X
'''-log(x+l)
-logo;
is less than 0-000125.
Here
I £" I = — — — — -r-jy — where 0 •
I w{\-n-) m 1
2 (a- + w'f log (a- + 1 ) - log X I
But by Taylor's theorem
in
log (.T + 1 ) - log X = ^ , ^^,„ where 0 < iv" -=1.
I V)(\- w) x + iv" I
^ F - 1000"
< 0-000125.
Example 3. — The following table gives the values of /j(.r) : —
X J,(.r)
1-1 0-47090
1-2 0-498-29
1-3 0-52202
1-4 0-54195
1-5 0-55794
1-6 0-56990
1-7 0-57777
Find the value of x corresponding to the value 0-55302 of Ji(x).
58 INTERPOLATION AND NUMERICAL INTEGRATION [cii. iii
21. Various Applications of Inverse Interpolation.
The method of inverse interpolation may also be employed to
obtain the maximum and minimum values of a function which is
given as a series of observations, or as a graph. In this case we
should differentiate one of the formulae of interpolation, say,
/(a + n w) =/o + n 8J\ + "^''-^^ 8^, + ...
Differentiation with respect to n gives
_ 2 w - 1 ,, ^
7vJ {a + n w) = 6/^^ + — — 6-/i + . . .
For maxima or minima we have
f {a + n iv) = 0,
so that n is given by the equation
2 w - 1 ^„ , 3 ?r - 6 7i + 2 „„ .
0-8/1 + —^ 5v; + ^ svg+...
w^iich is solved in the same way as the corresponding equation in
last section.
When it is desired to obtain the value of ?;, for which the
gradient of the function has a particular value o, say, the equation
for determining n becomes
. ^ 2 n - U._
The same method may be employed to obtain the solution of
any equation involving one unknown, whether it be algebraic or
transcendental. For we may tabulate the values of the expression
occurring in the equation and obtain their siiccessive differences.
The solution of the equation is then the same as finding the values
of the argument corresponding to zero values of the function.
Examiplt 1. — Approximate to the roots of the equation
1
21] MATHEMATICAL TABLES 59
First form the table.
X
fU-)
A
A2
A3
-3
_ 2
1-875
-2-5
- 0 125
0-125
-1-750
0-75
-2
0
- 0-875
-1-000
0-75
-1-5
- 0-875
- 1-1-25
- 0-250
0-75
-1
_ 2
- 0-625
0-500
0-75
-0-5
- 2-625
0-625
1 -250
0-75
0
- 2
2-625
2-000
0-75
0-5
0-625
5-375
2-750
0-75
10
6
8-875
3-500
0-75
1-5
14-875
13-125
4-250
0-75
20
28
18-125
5-000
0-75
2-5
46-125
23-875
5-750
3-0
70
The positions of the roots will be found by examining the column ./'('^')-
We then see that when /(.r) = 0, the value of x- lies between 0 and 0-5. It
will also be seen that f(x) is zero when x= -2, and that there is a
possibility of another root in the interval -2-5<a'<-r5. We shall
consider the various approximations to the root for which 0< x < 0-5.
Using the formulae of the first method of § "20 we have
/-/o
0-(-2)
^ 2-6-25
= 0-7619
"5A + i(%-i)5Yi
_ 2
~ 2-625- 0119 X 2-75
= 0-8704
INTERPOLATION AND NUMERICAL INTEGRATION [cii. iii
yv^^
"' " 5 /; + i ("o - 1 ) s-./'i + h ("o - 1 ) («2 - '2) 5v;
2
^ 2 -625 -0 0648x2 -75 + 0 0244x0 -75
= 0-8113
■• ~ 5./; + h ( "a - 1 ) Kfi + U "3 - 1 ) (n. - 2) 5^/;
2
"^ 2-625 - 0-0944 x 2-75 + 0-0374x0-75
= 0-8356.
The first four approximations to the value of n are therefore
n, = 0-7619
Wg =: 0-8704
w, = 0-8113
714 = 0-8356
But the tabular interval is one-half. Hence the first four approximations to
this root are
.^1 = 0-3810
Xo = 0-4352
xs = 0-4057
Xi = 0-4178
The correct value of the root is 0-4142. Subtracting each of the values
obtained from 0-4142 we get
0-4142 -.T, = +0-0332
0-4142 -r. = -0-0-210
0-4142 -.r., = +0-0085
0-4142 -.(■, = -0-0036
These differences show that each step gives a much closer approximation
to the root, and that by continuing the process we can get a result which
shall be correct to any given number of decimal places.
Examples of the second method, in which only first differences are taken
into account, will be found in any text-book which discusses the solution of
numerical equations.
Example 2. — Approximate to the real roots of the equations
(i) e-^-a;3 = 0
(ii) e^ + a;2-4 = 0
(iii) 0-5.r"5-121ogioa- + 2sin2a- = 0-921.
I
MISC. EXS.] MATHEMATICAL TABLES 61
MISCELLANEOUS EXAMPLES.
1. If y = xt%n--^x, calculate a table of first differences near .v — O'-tS for
increments 0'0(X»2 in .v.
2. The following table gives the moon's right ascension
Date. R. A.
1914, March 5d. Oh. 5h. 3m. 7 -26 8.
10 h. 26 m. 58 -96 8.
20 h. 51m. 14-83 s.
6d. 6h. 6h. 15 m. 48 -96 8.
16 h. 40 m. 34 -76 8.
7d. 2h. 7h. 5 m. 25-36 s.
Complete the table for every hour from March 5d. 20 b. to March 6d. 6h.
3. Prove that in a table of logarithmic tangents to base 10, the difference
for one minute in the neighbourhood of 60° is 0-00029.
4. Determine the values of (1) log tan T 29' 33"
(2) log sin r 15' 1-2"
5. Show that logjo 1031 = 3" 013 -258 665 "283 516 546 909 664.
6. Show that the number whose logarithm to 20 places of decimals is
given by 0- 004 892 890 303 "239 039 78 is 1 -01 133.
7. From the following table determine the value of x for which /(x) has
the value 0-2-2389.
X fix)
1-5 0-51183
1-7 0-39798
1-9 0-28182
2-1 0-16661
2-3 " 0-05554
8. Approximate to the real roots of the equations
(i) .re' + 2.v-5 = 0.
(ii) 6^-e-* + 0-4.r-10 = 0
(iii) 2.r3'i-3.T-16 = 0
(iv) 2-42x- - 3-15 log^.f - 20-5-0.
9. From the series
berx = 1
2-2 , 4-2 1- 2'2 . 4^^ . 6- . 8- 2'- . 4- . 6- . 8^ . 10- . 12-
eompute a table of ber a; from .r:=00 to .r — 50 at intervals of 0-5, correctly
to 6 places of decimals.
Hence, by subtabulation, compute a table of ber.^- from a; = 2-0 to a: = 3-0
at intervals of O'l.
62 INTERPOLATION AND NUMERICAL INTEGRATION [ch. hi
10. The confluent hypergeometric function W^ ^^ (:;) is computed for large
values of z from its asymptotic expansion
I 71=1
while for small values of z it is computed from the formula
r(-2m) r(2m)
where
[ \ + m-k {\ + m-k){% + m-k) \
M,, ,n (^) = zi+m ,-fc|l + 1, (2m + l) ^ + 2! (2m + l)'(2m + 2) ^'+ -j _
Compute tables and draw graphs of W^ ,^ (:) for positive values of z in the
cases
y!;=-3-l, wi=+2-2
k^-O-l, m=:+0-2
^=+2-2, m=+0-2.
CHAPTER IV
NUMERICAL INTEGRATION
22. Introduction.
The method of evaluating the definite integral of a function
from a series of numerical values of the function is called
Mechanical Quadrature or Numerical Integration.
If a function is given by its graph, the simplest way of deter-
mining an area bounded by the curve, two given ordinates and a
given abscissa, is by means of a planimeter. This method would
naturally be used in the case of the indicator diagrams of steam,
gas, or oil engines, or the cards from certain hydraulic motors,
or the stress-strain diagrams drawn by various types of testing
machines. The accuracy of the result obtained would in such
cases be as good as is required, but it is to be remembered that
graphical methods are not susceptible of great refinement. If,
therefore, considerable accuracy is required, or indeed in general
whenever the function is specified by a table of numerical values, the
method of numerical integration is preferable to the use of the plani-
meter. The method of nujnerical integration must also be resorted to
when a function is specified by its analytical expression, but cannot
be integrated in terms of known functions by the methods of
the Integral Calculus ; e.g., the elliptic integrals belong to the class
of functions which cannot be evaluated by the elementary
methods of the integral calculus.
The formulae necessary for numerical integration are derived
from the formulae already established for interpolation. As in
the case of interpolation, the order of differences which must be
taken into account, will depend entirely on the rapidity with
which the differences decrease as the order increases. It need
hardly be said that unless the convergence of the series can be
ensured the process is of no avail. Consider, for example, the
64 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
determination of the area of a semicircle, the ordinates being
perpendicular to the diameter. No matter how many ordinates
are taken (i.e. no matter the order of differences employed), the
approximation to the area by means of many of the formulae to be
hereaftet established will not be a good one. But if we evaluate
the area included between an arc of a circle, less than a semicircle,
a diameter, and the two perpendiculars from the extremities of the
arc on this diameter, we can, from the same formulae, obtain an
approximation to the true value with any degree of accuracy we
please. The reas(m for the failure in the first instance is that the
differential coefficients of the ordinates at the limits of integration
are infinite, and as a result the series is not convergent. For this
same reason we must assume that the function represented by the
observations, and its differential coefficients up to the order of
differences employed, are continuous ; a condition which is not
, necessarily fulfilled in natural problems, as may be seen by an
examination of an ordinary barometric curve.
23. Evaluation of Integrals by the Integration of
Infinite Series Term by Term.
Before applying the methods of finite differences we shall give
an example showing how in many cases the value of the integral
of a function which is analytically known may be obtained by
expanding the function in a uniformly convergent infinite series,
and then integrating it term by term.
Example i. — Consider the complete elliptic integral
TT
dx
^<^--l.V
fc- sin" X )
= 1 dx (l - k- sin" x) -
Jo
Assuming that | A; | < 1 and expanding the term (I -k-sin-x)'^
by the binomial theorem, we have
TT
^(^■' I) = [ "c^ x(l + l k- sin- X + .57-77 ^' s^^' *
Jo ^ " -". ^ .
1.3.5
2^ 3 !
¥' sin" X +
23] NUMERICAL IXTEURATION 65
1 . 3.5 ... (2n-l) -
But
Hence
I .siir" 6 dd = — ' — ^ ^ —
J, 2. 4. 6... 2 71
If k is small this series is very suitable for obtaining the
value of F ; but it is evident that if k approaches unity a large
number of terms must be included to get a fair approxima-
tion to tlie value of /''.
Let k-0-\. Then to obtain F correct to 4 decimal places it is only neces-
sary to iaelude two terms of the series. For
/ ^\ ^ S 025 0140623 )
^V^l.")=^il + 10Cr+ 10000 +-1
^ ^ {1 + 0-0025 + 0-0000140G25 + .. }
= ^{1-0025}
= 1-5747.
Example ;?.— Find the value of F (k, ^) when ^• = sin 10"
Example 5. — Show that I Jl-lc^sin-x is equal to r54415 when k has
the value sin 15°.
2 fo-o 2
Example 4. — Show that the value of the integral — p^ I c'^'dx is
V TT Jo
0 '60386.
24. A Formula of Integration Based on Newton's
Formula.
We shall now show how integration may be performed by the
aid of interpolation-formulae. Integrating Newton's interpolation
formula
J (a + 71 w) =/o + n 8/ + !ii!Lzi) S2 /• + _ _
66 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
with respect to n, we have
1 /(a + ?i ?r) d7i = /o I d ?i + 0 /■ l/i o? w + — o'./i I w (?i - \)dn
+ ^Sf/;, [n{n~\){n--2)dn+...
He
+ ...
Jo
f /(a + « .r) dn =f^_, + 1 8y;_, - J. 8V; + ^ ^'fr+^ - r^o ^Vr-
J r-\
Adding these equations we have
/{a + n?c)dn= ./o +/i + ••• +/r--i
But
8/1 + 8 /;+...+ o/;_, = /; -ji +./; -./; + . . . +. /; -/-i
87; + 8v;+ ...8v; = oy;-8/; + 8y;-8y; + ... + s/;^,-8y;_,
=sy;+.-8yi
24] NUMERICAL INTEGRATION 67
Hence
J\a + n ?(■) d n = hf^ +/i +/, + . . . +f,_, + 1 /;
-7^(sv;.+.-sv;)+
Putting x^a + nn- this gives
1 f «+'•«■
—I ./'(x ) c£ X = If, +y; + /.; + . . . + /;_ J + 1 /;.
- A (^;+i - « A) + -v (sv:-+: - ^\t\)
-T^^(sv;.+.-sv;)+
Tliis expansion will give the value of the definite integral, pro-
vided it is possible to obtain values of f{x) outside the limits of
integration in order to obtain the differences.
When this result is written in the form
1 f"+'-'"
—J ./-(.r) dx=it\ +./; + ... +/-,) + h (./;.+. /o) - iV (s/;+i - 34)
it is easily seen that the coefficients are those of a; in the expansion
, 7- : = 1 + Jo X + A^ X- + A., x^+ ...
log(l+a;)
Clausen * has computed the first thirteen of these, viz. : —
1 ^ 33953
Ao=+— A,
2 ■ 3628800
1 . 8183
^'12 ' ' 1036800
1 , 3250433
A„= +— A,
J4 " 479001600
19 , 4671
'^'~ 720 ^'° ■ 788480
3 , 13695779093
As= +T-7rX ^11 =
160 '' 2615348736000
863 . 2224234463
A
60480 '- 475517952000
275
54192
* Jour, reine angew. Math., 6 (1830), p. 287-9.
68 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
In the same memoir it is shown that when the limits of
integration are a - h n- and a + (r+ i) ?t- the coefficients are tliose
in the expansion
1 + B.x" + B.gif + B.x^ + ...
J {\+x)\og{\+x)
These are also tabulated.
The formula obtained above must be subjected to another
transformation when values outside the limits of integration cannot
be obtained.
iSTow
(i) 8/;,+, = 5;/:.+s/._x
(ii) sv;^;-2sv;-+sv:-i=sy.
i.e. s=/.^,-2 8y,= -5y,_, + s*/^.
(iii) ^/r=^A.+^-^^fr+l + ^^fr-l~^fr-^
i.e. ^fr+% = ^Vr + 3 ^.^, - 3 sv;_, + i%_^
Substituting these values in the formula we have
^£^' f{x) d X = \f, +f, + /; + . . . + /,_: + hfr
- ^i^fr-i - sy.) - tIo (sy.-. + sy )
This form involves only diiferences which can be calculated
from values of the function which are within the limits of
integration, and is the form to be used when no values outside
those limits can be obtained.
Example l.—lo determine it from the equation
I
4 -J^ 1+.^
using
differences of J^.
Here
f{^) =
1
/o=/(0) =
1;
/:
11
= 0-1; r =
10
:/(l)=.0-5;
Inspection of the formula will show that it is not necessary to form a
complete table of differences. The highest order of differences retained will
be the 4th, and as d^/.^ and oVs are the only two of this order required, the
central portion of the table may be omitted as in the following scheme : —
0-0
0-1
0-2
0-3
0-4
0-5
0-6
0-7
0-S
0-9
1-0
fix)
1 -00000
0-99010
0 96154
0-91743
0-86207
0 80000
0-73529
0-67114
0-60976
0 •55-249
0-50000
NUMERICAL INTEGRATION
A A-
- 990
-1866
-2856
-4411
-5536
-6415
-6138
5727
-5249
From this table we have
^/„ = 0-50000
/i ^ 0-99010
/. = 0-96154
f, = 0-91743
/, = 0-86207
/g = 0-80000
/h = 0-73529
/- = 0-67114
/s = 0-60976
/g =^ 0-55-249
i/i„- 0-25000
Sum = 7-84982
Hence
1555
1125
311
430
-277
411
478
134
67
x0 05-249
-0-00990
Diff. =
-0-04259
5Vi -
4 0-00478
-0-01866
Sum =
-0-01388
0 00067
0 00311
Diff. =
-0 00244
5v; =
5V, =
-0-00067
0-00119
Sum ^
0-00052
10
Jo 1 + ^^
7-84982
19
72~0
0 04259)
-0-00-244)
A*
119
67
(0-00052)
70 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
= 7-84982
355
58
6 - 1
= 7-85400
IOtt
—T- = 7-85400
TT = 3-14160.
Example 2. — Determine log« 2 from the equation ^ loge 2 — I \ -
using differences of -^ .
25. Case when the Upper Limit is not a Tabulated
Value.
Ill the previous section it has been assumed that the quantity r
in the upper limit was an integer, so that the tabular value of f{x)
was known when x had the value a + r tr. When the upper limit
is a + r'u; where r<r'<r+l, and the lower limit coincides with
one of the values of the argument, so that y (a;) is known when
x = a, a + n; ... a + 7- n; a + r + 1 w, . . . the value of the integral
may be obtained as follows : —
First determine, as in last section, the value of
ra + ru.
j f{x)dx.
Then, dividing the range {a^-ru\ a + r'ic) into a convenient
number of equal intervals, interpolate the values of /(x) for each
of these sub-intervals. The value of
f"+'- »■
Ja + rw
)dx
may now be obtained by making use of these new values of f(x)
and their differences.
When the lower limit does not coincide with a tabular value of
the argument, a similar process will enable the value of the integral
between the given limit and the nearest tabular value of the
argument to be obtained.
25] NUMERICAL INTEGRATION 71
Another method which might be employed is to obtain the
value of I f{x)dx for different ranges of integration, and then
determine the value for the given range by interpolation.
J 20°
20'
Example 1. — Evaluate I cosx'rfx, given the values of cus a; when
J 20°
a: = 20°, 22°, 24°, 26°, 28°, 30°, 32°.
Forming the table of differences we have
X cos.f A A- A^ A*
20° 0-9396926
- 125087
22° 0-9271839 -11297
- 136384 166
24' 0-9135455 -11131 16
- 147515 182
•26° 0-8987940 -10949 9
- 158464 191
28° 0-88-29476 -10758 16
- 169222 207
30° 0 8660254 -10551
-179773
32 0-8480481
1*30"
j20^.'
Firat Method. — We shall first of all determine | cos x dx. Substituting
in the formula
i- P' "/{x) dx^hf, +J\ + . . . +/-, + 4/- - ^ (5/; _ . - 5/p - ^ (5^/._i + 8- A)
we have
1 P°" 7 f..ifin«^ftQ , r-169222-1 , T - 10758-] x,, T lOn
0-9135455
0-8987940
0-88-29476
0-4330127
= 4-5257896
(•30"
I cosa;(Za- = 01579798.
J 20"
Now by using Newton's formula for backward interpolation
n(n-l) „ . n(n-l){n-2) ^^ .
/(a - 71 ir) = /;, - n dj\ j + - .-^-j — 5-f_ ^
3:
determine the values of . /'(.(,), I.e. cos x, when x has the values 30° 20', 30° 40',
31°, 31° -20', 31° 40'. We thus obtain the following table : —
INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
30^0'
0-8660254
- 29235
30" 20'
0-8631019
- 295-28
30° 40'
0-8601491
-29818
sro'
0-S571673
-30109
3r20'
0-8541564
- 30397
31 = 40'
0-8511167
- 30686
32=0'
0-8480481
om this ta
ble
we obtain
1
•3]
" 20'
30
COS .f dx -
^0-4330127 -i\
0 0058178
0-8631019
0-8601491
0-8571673
0-4270782
-3-4405189
f
Zl" 20'
COS X dx
-0-0200161
•293
290
291
288
•289
r- 301091 _ , r---29n
L + 29235J ^*L-293j
Adding these results together we get
rsi" -20'
cos .i(/.t= 0-1779959
Second il/e//iO'/. — Obtain by the general me I hud thu values of I cos .i- cZa;
when the lower limit is 20 and the upper limits are 22', 24 , 26 , -28', 30% 3-2'
respectively. We thus obtain the following table : —
X
1 cos X dx
A
A^'
A3
20°
0 0000000
325865
22''
0-0325S65
321300
-4565
390
24-
0-0647165
316345
-4955
-385
26°
0-0963510
311005
-5340
-381
28°
0-1274515
305284
-5721
-370
30°
0-1579799
299193
-6091
32°
0-1878992
25] NUMERICAL INTPXiRATIOX 73
Using Newton's formula of backward interpolation we obtain the value of
the required integral, viz.
f3i<'20' ^Ui-1) A(4-i)(i-2)
cos.fdt—O-1878992-t (299193)+ ., , ( -6091) - 3 , ( " -^'O)
= 01779961.
The values obtained by the two methods agree to 6 places of decimals.
The result correct to 8 places is 0'17799.398, but as the last figure is always
forced, the one method can hardly be assumed better than the other. Probably
the latter method is to be preferred, as in the former the division of the range
between the given upper limit and the new upper limit into suitable intervals
may give rise to values of the argument consisting of many decimal figures
which would increase the amount of labour involved in the interpolation.
r 30" -20'
Example J. — Evaluate I sinxdx, given the values of sin.r when
J 20"
x = 2{)°, 22°, 24% 26% 28% 30% 32 .
f4-5 dx 1
Example 3. — Evaluate 1 . if the values of , , , are given for
integral values of x.
26. A Formula Based on Bessel's Expansion.
Another series for the definite integral is obtained by integrating
Bessel's formula of interpolation,
, „ . n(7i-]) ^, . V4(M - l)('/i - -i) ^, .
J{a+ n ?c) = II J ^ + (n-h) 8j^ + \, , /x ir/^ + — S"/,
n{ir-l}{H - 2)
-\ — n 0%/, + . . .
Putting x = a + nw we have by integration
— f(x)dx= t\a+niv)dn
= l\fA dn + 8jA {n-h)du + ii6-J^\ —^, dn
<
0
11 (71 - \) (n - h)
•i/i(?r-- l){n-2)
a ii
74 INTERPOLATION AND NUxMERICAL INTEGRATION [ca. iv
Similarly,
^ J a+K ' ' '
^^ J, ,+ (,-!,«,
Hence
w
+
or, since
SV; + 5\/; + . . . + SV;._i = 54 - 6y- + Sy- - Sy, + . . . + 8y;_ , - Sy;_ .
and therefore
we have
1 f ''+'•«■ , .
—J ^ /(x) d X = -jy; +yl +/;+••■ +/.-1 + h/r
+
But
26] NUMERICAL INTEGRATION
and, similarly,
Hence, finally, we obtain the formula
1 r«-^'-'"
— I ^ J (x) (/ X = ifo +j\ +/, +... +/-1 + y'r
75
f lOo dx
Example i.— Calculate — to 8 places of decimals.
J 100 •*-
1
Here /{x)
-, /o=100, r^o, w^l.
In order to obtain the dififereaces required, it will be necessary to form
/(■c) for all integral values of x from x- = 9S to ;c=107. These may be
obtained from Barlow's Tables. We thus obtain the following : —
.(-■
/(a-)
A
A-
A^
98
0-010204082
- 103072
99
0 010101010
- loioio
2062
-62
r — 1(A)
0-t)10(MH)00(»
( - 100010)^
- 99010
•2(3(30
( - 60)
-58
101
0-009900990
- 97(368
1942
-58
102
0-0(^9803922
- 95184
1884
-53
103
0-(W9708738-
- 93353
1831
-53
104
0(309615385
- 91575
1778
-51
^— 105
(J -009523810
(~ 90712)
- 89848
1727
(-49)
-47
106
0-009433962
- 88168
1680
107
0-009345794
.Substituting in the formula
f I'J.J dx
-- = ^ .'o + -'i + -/--i + J'^ +J4 + hA-^ (m S/, - p. 5/o) + ^Vo (/^ 5% - M 5Vo)
Jioo ^
we get
f 105 dx
I - rr u 048790940 - 0 -(3(;k3000775 f 0 -( )0(XM3(3(30<3
Jiou ^
^0018790165.
The value correct to 10 places is 0 0487901642.
76 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
Note. —It is to be remarked that the expression
l/o +J\ +/2 +/) +/4 + 4/5 = 0 -048790940
gives the answer correct to six decimal places ; this approximation is quite
sufficient for most practical purposes.
f3 dx
Example 2. — Approximate to the value of the integral 1 ,, Q_t>„2 1 o\ '
jo \ y^ -X" + —J
27- The Euler-Maclaurin Expansion.
In the foregoing the value of the integral is expressed as a
series involving differences. In many cases in which the function
is known analytically, it will be found more convenient to use what
is known as the Euler-Maclaurin expansion which employs
'f differential coefl&cients instead of differences. The rapidity of
' convergence is of the same order as that given by Bessel's formula.
It may be obtained as follows : —
Since
= ^[J{a+w)-j\a-w)]
W'
= h {/(«) + "^/' («) + y/" («) + •••
-fia) + wf (a) - ^/" (a) + . . . }
= w/'{a)+ '^/"'(a)+...
and similarly
/x Sy, = tvj ' {a + r to) + — y '" {a + r to) +
w''
Ij. B^f,. = vff" (a + r iv) + —./"*'■' (a + r w)
4
we have by substituting in the formula
If*'"
fix) dx=\j\ +j, +y; + . . . +./;._, + y.
26, 27] NUMP:RICAL INTEGRATION 77
the result
^j"^' " ./■ (^^ dx = 1 /; + /; + .. . + /;_, + 1 /;.
-TroiTo"'-M/.'^'-/o'") + ...
This is the EuIerMacIauriii* expansion. The coefficients of the
various items in it are proportional to the well-known Bernoullian
numbers. For if f we put f{x) = e^-" in the formula and at the
same time chansfe the limits to a ~ w and a + tv we iret
e-" dx = ifo +/] + y: + A w {fj -/;) + b «- (/;'" -./;'")
where A, B, ... are constants.
But
— I e" " dx ^ — I
W Ja-w W J-„.
dx
— (e'" - e-'").
w
Hence
— (e'" - e-"') = .i/o + /; + 1/, + ^ 7r (/,' -./;/) + B >r (./:'" - /;,'") + . . .
= I e-'" 4-1+1- e'" + A w (e"' - e'"') + B ir^ (e'" - e~'"} + ...
from which it is easy to deduce that
^; ::. :, = I - A iv - B >r* - ...
w to , , „ r. ■
T— . cot -T—. = \ - Aw- - Bijcr - ...
'1% 2 1
But the Bernoullian numbers are the quantities ^j, B.., ... in the
expansion
and therefore .4, ^, ... are proportional to B^., B.,, ...
* For a rigorous discussion, c/". Whittaker and Watson's Modern
Analysis, p. 128.
t Poisson, Mem. de I'Acad. des Sc, lb23.
1
1
78 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
Example i. — Calculate — to 8 places of decimals, and check the
J 95 '^
result by computing it as - loge(l - ygy) by the logarithmic series.
Here /o = tfV, ^"=^, r = ^-
Then
Jp5 ^ 2 ■ 95 ' 96 97 98 99 2 100
T^ ^ 100- 95-7 ^ - " V " '*»' 95V
w -O-OOOlOO0OO\i/ -0 000000010 ^
- 0-005-63158 1^ + 0-000110803/ ^'-"\ +0-000000012/
0 010416667
0-010309278
0-010204082
0-010101010
0-005000000
= 0-051294195
-0-000000900
= 0-051293295
Again
, /, 5 \= ^^. (0-05;^ (0-05)3 (005)^ (0-05)3 (o-05)«
-i°g^(i-ioo; o'^«+-T-+-^-+— r-^~5-^-6~~
= 0-050000000
0-001250000
0-000041667
0-000001563
0-000000063
0-000000003
= 0-051293296
The two results will be seen to agree to 8 places of decimals.
fl05 d.r
Example 2. — Evaluate I -7- to 8 places of decimals.
IT
Example 3. — Approximate to tlie value of 1 co^-xdx.
28- An Exceptional Case-
Legendre* remarked that if the odd differential coefficients above a
certain order take the same value at both limits, the formula fails to give an
accurate value of the definite integral. For instance, if the integral under
* Fonctions Elliptiques, Vol. 11., p. 57.
28, 29] NUMERICAL INTEGRATION 79
discussion were the elliptic integral of the second kind taken between limits
0 and 7r/2, viz.
rTT's
I J l—JS^ sin^ X dx,
Jo
all the odd differential coefficients become zero at both limits, and thus the
value of the integral would appear to be given by
i/o +/i +/2 +•• ■+./;-! +i/;-.
But the value given by this expression is not in close agreement with the
true value of the integral. The reason of this is that the numerical
coefficients introduced in the successive differentiations increase without
limit, so that each term takes the form oc x 0, which is, of course, inde-
terminate.
Better methods of evaluating the three kinds of elliptic integrals will be
given in another tract in the present series.
29. The Formulae of Woolhouse and Lubbock.
We shall next derive from the Euler-]\Iaclaui-in expansion a
formula which is of use when it is required to evaluate the sum of
a large number of terms.
If the interval between consecutive values of the argument in
the Euler-Maclaurin formula is w, we have
and if this interval is w j m the formula becomes
- f{x) dx= ./; + /■ , +/ , + . . . +./;. - 1 (/o +/,.)
n- J f, mm
-A^(./;.'-y:')+Tk|(./;--y;"')+...
Subtracting the latter from m times the former we deduce
/o+/i +J I + .-• +J>m{f,+J\+... +/.) - '"-^ (/;+/.)
12 m^'^' '^'■>
m* - 1 tv' ^ .,„ ^.,„^
80 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
Tliis is Woolhouse's Formula. OV)viou,sly by using the right-
hand expression instead of the left-hand one, the labour involved in
determining the sum ^ 4-/2 +./" 1 + ■■■ + A- is enormously reduced.
A somewhat similar formula, in which differences are employed
instead of differential coefficients, is due to Lubbock. Expressing
the derived functions in Woolhouse's formula in terms of the
successive differences of /q, J\, ■■■/,■ we obtain the formula in
question,
m — \
12m '"5 2 24 ?/t
(m-^-l)(19m^-l)
I 20 7H' ' ■' ^
(m--l){9m-- 1)
480 m"
{8V>-2 + ^V2)
3j0 1
E.rample 1. — To determine 2 — •
'/! = 300 "
Let w = 10, ?o=10.
Then from Woolhouse's formvila we have
350 J_ M +J_ +J_ +J_ +J_ + J_\
rMw ^ "^^^300 310 320 330 340 350/
2 VSOO 350/ 12 V 350- 300V
= 018512761-0-02785714
2432
= 0-15724615
Losing Lubbock's formula we get
2 IT = 10 (/o +/i + • • ■ +^5) - s (./'o +./;) - -/-vo (5 n - 5/a ) - i?^ roV\ + s% >
n = 300
720000 ^ -'-V •'4' 480000 ^"-'-i ^ "-'2^
^0-18512761 -0027S5714
193S
487
3
3
= 0 15724616
29, 30] NUMERICAL IXTEGRATION 81
To test the accuracy of these answers, the reciprocals of the numbers
300, 301, 302, ..., 350 were taken from Barlow's tables, and summed by
means of a comptometer. The result was found to be 0'15724616. The
closeness of the approximation is remarkable.
175 1
^.rrt??ip/e ^.—Obtain the value of 2 — •
n = ioo "
Examplf 3. — By giving n increasing positive values, show that
T 1
30. Formulae of Approximate Quadrature
The expansions which have been obtained in the preceding
articles furnish the value of an integral to any required degree of
accuracy, provided a sufficient number of terms is taken. We shall
now consider certain formulae which are in theii' nature only
aj>proxiinate, and which, though frequently useful, cannot be used
where great accuracy is required.
Suppose that it is required to integrate a f unction /(.x) between
limits, which may, without loss of generality, be taken to be - 1
and + 1. We shall endeavour to obtain expressions of the type
^o/(^'u) + H,f(K) + ...+ iiJiK)
which represent the integral as closely as possible, where Aq, h^, ... h„
are (n+l) values of x within the range of integration and h^, h^, ... h,„
Bq, //j, ... Zr„ do not depend on the function /(.r).
Let us first choose B„, B^, ... B,„ so as to make the formula
strictly accurate, so long as f{x) is a polynomial of degree less
than n + \.
Let F{x) = (.X - /g {x-h,) ... (x - h„).
F(x)
Then -— is a polynomial of degree less than n+l, so that the
X - h,. •' °
F (x)
formula must be strictly accurate when for /(xj we put — ^ .
" X - h,.
This gives
F(x)
B.L
BrF'{K).
82 INTERPOLATION AND NUMERICAL INTEGRATION [cii. iv
Hence, whatever hg, //j, ... h„ may be, in order tliat the above con-
dition may be satisfied, the values of Hq^ IIj, ... H^ must be
given by
In practice it is convenient to make the intervals between
successive values of the argument equal to each other. Suppose,
then, that A,i, ^i, ... A„ are chosen so as to divide the whole range
of integi'ation into n equal parts.
Then A^ = - 1
1
1
K =
-i + —
n
/i., =
-1 + 1
n
H,
'■ {h, -K) ... (h, - /i,,„, ) {h, - A,+, ) . . . (A, - A J
X (x- Ao) • • ■ (•>' - K-i)(^' - h,.+i) ■ . . (x - h„) i
r 2(r-l) 2- -2 -4 _ 2 (n - r)
?i 7J n n ti n
( - 1)"-^
("— Y\-! (ri-r)!
Now^ let t = -^{x-\-l).
Then
//.= ^~^^" '^^ , f%(^-l)... (<-r+l)(<-r-l) ... (t-n)dt.
n . r\ [n - r) I Jo
Hence, finally, if we put
1*0 = a, ... h^ = a + r w
30] NUMERICAL INTEGRATION
wo obtain
\y{x) dx^j: II, f {a + r w)
J a )-=o
whei'i
//,= ^-^— ^ -^ rV<-l) ... (t-r+l){t-~r-l) ... {t-n)dL
r I {71 - r) ! Jo
This is known as Cotes Formula of Quadrature. By means of
it we may calculate the value of the integral without first of all
forming a table of diifei'ences.
Special Cases.
(i.) Let?i=l. Tlien
f 0 4 w
f{x) d X = nj{a) + II, /{a + n-)
— iv r^
where H, = — — - {t-l)dt = h rr
0 ! I ! J 0
Hence
[^'\f{x)dx = ^ {f(^a)+f{a + iv)].
This is known as the Trapezoidal Rule.
(ii.) Letw = 2. Then
where
r^' "'./•(.r) d X = HJ{a) + HJ (a + w) + H.J(a + 2 w)
Hence
1+2 19 ?0 ,, 4 7t'
./ {x) dx =^ — y («) + ^y (a + «;) + y y (a + 2 w).
This is called Simpson'' s * i??i/e or the Parabolic Rule.
* It was, however, first given by James Gregory.
84 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
Corollary. — If we divide the whole range of integration, say
{a, b), into 2 n equal parts, and apply Simpson's Rule to each of
them, the value of the integral may be expressed as
/ = -g^ {yo + 2/2n+2(?/o + y4+ ...+?/.2„-:)+4(t/,+2/:;+..-+?/2„_l)}
where i/q, y^, ... y..,, are the values oi /(x) "when x has the values
h - a
2 n
form. The rule may be stated thus : —
Divide the area into an even number of strips by equidistant
ordinates y^„ t/j, ... v/.j„ : to the sum of the extreme orJinates add
ttvice the sum of the other ordinates with even suffixes, and four
times the sum of the ordinates with odd suffixes, and multiply this
total by one-third of the common distance betiveen the ordinates.*
(iii) Let n = 3. Then
fa+3w)
f{x) dx = HJ{a) + II J {a + ic) + HJ{a + 2iv) + HJ{a + 3 to)
where H, = -^^^^^{t-\){t -•2)(t -?,) dt = %w
Hence
T"^ VH dx ^%w {f(a) + 3/(« + tw) + 3 /'(a + 2 w) + /{a + 3 w)],
which is often termed the Three- Eighths Ride.
(iv) Let n = 6. Tlien it may be shown, as in the above
special cases, that
f(x)dx = lU {Uf(a) + 2\Gf{a + w) + 27 f{a -{- 2 tv)
+ 272/(a + 3 m;) + 27 f {a + 4 w) + 216 f {a + 5 w)
+ 41/(a + 6M7)}.
* The reader is warned that the notation varies in the different text-
books, so that in all cases aa examination of the notation used should first
of all be made.
i:
30,31] NUMERICAL INTEGRATION - 85
Adding
ih S!/*(« + 3 w) = -^Ijj {/{a) - 6/(a + w) + 15/(ffl + 2 w)
- 20,/"(« + 3 w) + 1 bj\a + 4 w) - 6/(a + 5 w;)
+/(a + 6M;)}
we obtain
Ttt + Ow
/(,«) ci X- + ^i^ S7 (a + 3 w;) = ^io {^V\<^) + ^lOy^a + t^)
^^ + 42/(a + 2 w) + 252/(« + 3 m;)
+ 42/(a + 4 w;) + 210/(a + 5 w)
+ 42/(a + 6to)}
= TO {/■(«) + 5/(« + ^<^) +/(« + - «*^")
+ 6/(a + 3 Mj) +/(« + 4 ztf)
+ 5/(a + 5 w) +/(« + 6 i«) } .
Neglecting the term -^\-^h^ f {a-^?> lo), which will be fairly
small, we get Weddl&s Rule of Quadrature for 6 intervals.
J{x)dx = j% { /(a) + 5/ {a + lo) +/(a + 2 w) + 6/ (a + ow)
"^ +/{a + 4 w) + of {a + 5 iv) +f{a + 'oiv)].
31. Comparison of the Accuracy of these Special
Formulae.
We may compare the accuracy of these formulae by evaluating, by each
n fix
of these methods, the definite integral I ~^^, , the correct value of which is
log^7 = 1-94591.
(i) Using the Trapezoidal Rule six times, i.e., dividing the region into
six equal portions and applying the rule to each of them, we get
f7 dx
— = \f{<i) +J\a + w) +f(a + 2iv)+ ... +f{a + 5w) + y'(a + 6w)
= 2-02143.
(ii) Using the Parabolic Rule repeated three times we get
n dx
= ^1+1 + 1 + 1 + 1 + 1 + 1)
= 1-95873.
86 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
(iii) Applying the Three-Eighths Rule twice we have
n d.c
1 T " a '-^("^ + ^J^'^ + "') + -Vi" + 2w) + 2/{a + 3 ic) + 3./'(a + 4 ir)
•' +3/(a-f 5t(j)+/(a + 6j«)}
= t{i + § + l + ! + § + if + i}
= 1-96607.
(iv) Weddle's Rule gives
\[^ = A{i + f + Hf + i + f + i}
= 1 •95286.
None of these gives the result with any high degree of accuracy, the
errors being
(i.) U -07552
(ii.) 0-01282
(iii.) 0-02016
(iv.) 0-00695
Example I. — E
evaluate j^^ by each of the methods of § 30.
Example ^.— Plot the curve y- = 2 .r -f 25 aud find the area in square inches
between the curve, the ^-axis, and the double ordinate at ,r= 12, the unit
being one inch.
Examjde 3.— The ordinates of the boundary of the deck of a ship are
8, 30, 36, 40, 42, 42, 40, 38, 34, 28, and 8 feet respectively, and the common
distance between them is 21 feet. Find the area of the deck.
Example J,. — The velocity of a train which starts from rest is given by the
following table, the time being reckoned in minutes from the etart and the
speed in miles per hour.
t 2 4 6 8 10 12 14 16 18 20
mlh 10 18 25 29 32 20 11 5 2 at rest.
Estimate the total distance run iu "20 minutes.
32. Gauss' Method of Quadrature.
The accuracy of the method given in § 30 may be increased if
we no longer suppose the intervals between the successive values
/t,„ ^1, ... ^„ of the argument to be equal.
Let us now endeavour to choose //„, //j, ... //„, //„, h^, ... h„ in
such a way that the formula
j J{^dx= Hj\h,) + HJ{h,)^ ... +IIJ{K)
31, 32] NUMERICAL INTEGRATION 87
shall be accurate foi- all functions /(x) of dcgi-ee equal to or less
than (2 n+ 1).
Taking
fix) = F{x)^{.i')
where </>(«) denotes any polynomial of degree <(n+ 1), and where
F{x) = {x - h,} {x - h,) (x -h.^ ... {x - h„)
we see that the condition to be satisfied may be written
f F{x)<fjij:)dx = 0.
Now it is well known that if F,,^^ (x) denotes the Legendre
polynomial * of degree (?t+ 1), we ha-se
r P„+A^^}4> {■'>') dx = 0
and, conversely, this condition is sufficient to determine the
Legendre polynomials. Hence we have (save for a multiplicative
constant)
F{x) = P„_,,(x)
and therefore the quantities Z/^, /tj, ... h„ must be the roots of the
equation
^,,-1-: (^■) - 0.
The foUuwing are the values of the constants h^, h^, ... h , //„, H^, ... H
in the simplest cases.
(i) Whenn^l.
K =
1
"v^3
K-
1
(ii) When7i = 2.
/in =
- ^'k
li, =
0
/l, =
•Ji
H,
= 1
H,
- 1
H,
_ 5
— y
H,
= 1!
H.,
_ 5
— y
* Cf. Whittaker and Watson's Modern Analysis or Byerly'a Fourier's
Series and Spherical Harmonics.
88 INTERPOLATION AND NUiVlERICAL INTEGRATION [ch.
(iii) When?t = 3.
K= - s^( 'i + is v^30) H, = h- rk ^/30
'>i= - V ( ? - i's v'30) Hi = i + sV v/30
''i = >/{■?- Tft V30) ^2 = i + gV x'30
'^.i = s'( f + y^ V^30) i^:, - i - 3^6 V30
(iv)
When w = 4.
K= -
- s'({; + ^n'70)
K= -
• N'(^-i^V70)
h.,=
0
h,=
xAa-Av^70)
^4 =
N^(H6lTV^70)
^^ ^ 322 + 13 v/70
900
900
128
225
322+13\/70
900
322-13V70
^^= 900
^' 900
Example I. — Determine the val
(i) Using 3 urdinates
r' ^ = 5 _J__
J-] 3+^ 9 -3- ^1-4
in
J-l 3+«
9 ■ 3 9 • 3+ x'f
= 0-69312165.
(ii) Using 5 (jrdinates
3-22- 13v/70
1
3-22+ 13 V
900
1
3+ V(|-
70
6-3
1
900
1-28
+ 225
322-
0-693147157.
3- N^(^
. 1 3
3 +
-13V70
900
+ 6^V70) '
22+13v70
900
L
3 - s!{^ - h sHO]
v/70)
3+ v'(| + ff-ov70)
The correct value is 0-69314718, so that using 3 ordinates the error is in the
5th decimal place : using 5 ordinates, in the 8th decimal place.
^xa»ip/e i^.— Determine the values of (i) I -r4> '
MISC. Exs.] NUMERICAL INTECJRATION
MISCELLANEOUS EXAMPLES.
L The values of a function /(.r) for the values 1-050, 1060, 1-070, 1-080,
1-090, 1-100 of the argument are 1-25386, 1-26996, 1-28619, 1-30254, 131903,
1-33565. Use Newton's formula of quadrature to show that the value of
fi-inu
.fix)dx
J 1-050
is 0-06 472.
2. The ordinates of a plane curve are of the following lengths, 3-5, 4-7,
58, 68, 7-6, 8-1, 80, 7-2 feet, and they are 4 feet apart. What is the area
included between the two end ordinates? and what is the distance of the
centre of gravity of that area (i) from the extreme left-hand ordinate, and
(ii) from the base line ?
3. Approximate to the values of the following definite integrals : —
(i) r dx.^{S-.v + x^)
(ii) c/a
i:
IT
f? dx
4. Use (i) Lubbock's formula, (ii) Woolhouse's formula, to determine the
value of an annuity-certain for 36 years, interest being at the rate of 3%.
5. Verify that the area of the curve y~A + Bx + Cx- + Dx^ between the
limits x = h and x= -h is equal to the product of h into the sum of the
ordinates at ,r = A/^'3 and x= -h/^/S.
In the case of the curve i/=f(x) = A + Bx+ Cx'^ + Dx^ + Ex^ ■\- Fx-' verify in
like manner that the area between x = h and x= -h is equal to
{ 5f(h Jl ) + 8/ (0) + 5/( - h J J ) } /t/9.
6. Corresponding values of .r and y are given in the following table, the
unit being the inch. Find the volume of the solid of revolution formed when
this curve revolves about the .r-axis. Find also the centre of gravity of this
uniform solid.
0
4
8
12
16
20
0-50
2-58
4-41
5-40
5-10
0-50
90 INTERPOLATION AND NUMERICAL INTEGRATION [ch. iv
7. A body of 200 lbs. is acted on by a variable force F. No other force
acts on the body. When the body has passed through the distance x feet
the force in pounds is as follows : —
0
0-1
0-2
0 3
0-4
Oo
0-6
07
0-8
0-9
1 0
20
21
21
20
19
18 5
ISO
1.3-5
9
4-5
0
Determine the work done upon the body from .r = 0 to x = 0'4. If the
velocity is 0 when x = 0, what is its value when a; = 0'4 ?
8. The semi-ordinates in feet of the deck plan of a ship are respectively
3, 16-6, 25-5, 28-6, 298, 30, 29-8, 29-5, 28-5, 24-2, and 6-8, the common
interval between them being 28 feet. Find the area of the deck.
9. The "half" ordinates in feet of the mid ship section of a vessel are
12-5, 12-8, 12-9, 13, 13, 128, 124, US, 104, 68, Oo respectively, and the
common interval is 2 feet. Find the area of the whole section and the
position of its centre of gravity.
10. The area of a ship's load-water plane is 8000 sq. ft. ; the body
below the load-water plane is divided into six portions by equidistant
horizontal sections, 3 feet apart, whose areas in sq. ft. are respectively
7600, 7000, 6000, 4500, 2800, and 250. Find (a) the displacement in tons,
and (b) the number of tons which must be taken out of the ship to lighten
her 4 inches.
Printed by Lindsay & Co., 17 Blackfriars Street, Edinburgh
14 DAY USE
RETURN TO DESK FROM WHICH BORROWED
LOAN DEPT.
This book is due on the last date stamped below, or
on the date to which renewed.
Renewed books are subject to immediate recall.
"aTFe^saWJ
REC'D LD
T-^m-
HECEIVED
MAKixms-
'VMdyWpW
DEC 2 7 '66 -4 PM
LOH DEPT.
RECrP \ TX
TfWHM96e-
^0^^ar'6St0
^REOU-LO
wty-'S^^i
LD 21A-50m-9,'£
(6889sl0)476B
General Library
University of Calif.
Berkeley
Live Music Archive
Librivox Free Audio
Metropolitan Museum
Cleveland Museum of Art
Internet Arcade
Console Living Room
Books to Borrow
Open Library
TV News
Understanding 9/11