The American Term Life Insurance History Project

You are reading a page from An Elementary Treatise on Actuarial Mathematics by Harry Freeman (1932)
Part of the American Term Life Insurance History Project
Term Life Insurance

CHAPTER XVII
APPROXIMATE INTEGRATION

1. In order to obtain the area of a curve by the methods of the integral calculus as given in the preceding chapters two conmust necessarily hold. These conditions are

(i)the equation of the curve must be known; and in that event

the function

= u

representing the equation of the curve must be integrable.

In the theory of life contingencies these conditions are rarely satisfied. Rates of mortality, marriage, etc. are generally obtained from actual observations and the functions derived from these rates are seldom capable of expression in the form of a mathematical expansion. For example, if l_xbe the number of persons attaining exact age x in any year of time, the number living between ages x and x + 1 is given by

L_s= J 1_,s+tdt, .owhere x and t are independent.

Unless l~__t follows some definite mathematical law, we cannot evaluate the integral by the methods hitherto employed. By making certain assumptions, however, we can obtain approximations to the value of the integral which are accurate enough for practical purposes.Thus, if we expand l_x+tin terms of l_xand differences of l_x, we have

^lx+t = (^I+ ~1)^tl_x

tOl

to first differences.

But

Ol,.

x+i —

and

l„ — i

_f+l

d, ,

the number out of the l

persons who die before reaching age x

286 INTEGRAL CALCULUS
.. Lx = J ¹lx+t dt
0 r1
= J _o(l_x — td_x) dt as far as first differences
= Ctl_x — 2t^zd_x
L o
=lx—2dx
=lx^—(^Ix^—1x+₁) = 2 (^Ix + ^lx+i)•
This is an approximate value of L_x on the assumption that ^lx+e is a first difference function of l_x.

It is evident that a simple approximation to first differences will not be justifiable except in a limited number of instances : the conof formulae for approximate integration will necessarily depend on the data available and on the degree of accuracy required.

2. Simpson's rule.

Suppose that we are given two values only of u_s and that we represent these two values by means of the points H and K whose co-ordinates are (a, u_a) and (b, u_b) respectively. Suppose also that we have sufficient information to justify our representing the funcy = u_s by a curve of small slope and small changes of slope,
not affected by periodic changes. An infinite number of such mathecurves can be drawn to pass through H and K. The simplest curve passing through the points is a straight line, and for the purof interpolation there is usually no justification for adopting any formula other than a first difference
formula, as was done for example

in obtaining the approximation to O N b M Y
-----------
L_x in the preceding paragraph. If Fig. 36.

the straight line be drawn, the area cut off by the curve, the x-axis and the ordinates x = a and x = b will be that of the trapezium HNIl7K. The area of the trapezium is (ua + u_b) (b — a).

SIMPSON^'S RULE 287 b
Alternatively, this area is ( u_xdx.
-a
An approximation to the integral is therefore (ua + u_b) (b — a). If y = u_x= ^lx+s and the limits b and a are 1 and o respectively,
we have the approximation to J l_x+tdt given above, namely
0

i (^lx + ^lx+1)^.

Let us consider the effect on the approximation of introducing a third known value of u_x. In other words, suppose our data to be u_a, u_b, ue. Then there is still
an unlimited number of curves which can be drawn to pass through the three points H, K, L. The simplest of these curves will now no longer be a straight line, but will be a second difference curve. The form of this curve will be y = A + Bx + Cx², and since three points on the curve are given we
are in possession of sufficient data to

determine the coefficients A, B, C. ^Fi^g.37.

In practice the data will generally be available at equidistant intervals, and since we may choose our origin where we please, the problem is therefore reduced to one of finding the second difference curve which passes through the points (o, uo), (I, u₁),

(², ^u2)^.

Having found the equation of the curve we can then find the area enclosed by the curve, the extreme ordinates and the x-axis, by integrating the function thus obtained between the limits o and 2.

Assume therefore that y = ux = a + bx + cx² is the required equation.
Then, by substitution,

uo=a, u₁=a+b -!- c,

u₂=a+2b+4C,

z88 INTEGRAL CALCULUS
whence a = uo, ^b=(-^ut+4^u1-3^uo),

r ^c=2(^u2^—2161^-+0)^. 2

Also J _ou_x dx = t_o(a + bx + cx²) dx
= ax + 2bx² + icx3l 2 0
=za+2b+3c
= 3 (u₀ + 4u₁ + u2),
on substituting for a, b and c.

This is Simpson's rule for approximate integration. We can obtain this formula by alternative methods. For example, let

ux = (1 + 0)^x u₀ = u₀ + x,uo + ix (x — I) 0²u₀ as far as second differences.
    2    2
Then I_ouxdx=_Jo{uo+x,uo+2x(x-I)O²uo}dx
    ^J2
= J^._o {u₀ + x (Au₀ — Z,²u₀) + x²20²u₀} dx
[xuo + ix2(Au, — L^,2u0) + i^x3A2u0 0
2u0 + 2 (Xu₀ — 2A2160) + 3^L2u0 = 2u₀ + 2A11₀ + 1.1²U₀
,
and since    uo = uo,
u_l = uo + :duo ,
u2 = u₀ + 2:,u₀ + r,=uo ,
this reduces to (11.o + 4u₁ + u₂) as before.
A third method is to assume that the integral can be expressed in the form

muo + nu₁ + pu_t .

Then if y = a + bx + cx² we have eventually that J₀x = za + zb
;-Vic,
as in the first method.

CHANGE OF UNIT 289 Substituting in the assumed expression for the integral: ma-j-n(a+b+c)+p(a+2b+ 4c)=2a+zb+Rc.
Whence, by equating coefficients of a, b, c and solving the resulting equations,
^m=3; ⁿ⁼p⁼3^.3. Change of unit.
r2
The formula J u_xdx = 3 (uo + 4u₁ + u₂) is an approximate 0

formula obtained by considering unit intervals. If we wish to transform the formula to a form in which the interval is changed to, say, n, then the given values of u_x will be uo; u_n; u_2n. Our new variable is z, where z = nx.

dz = n dx.
The new limits are evidently o and 2n.
The formula is
j2n ₁(
- tl_Zd2 = (^uo + 4^un + ⁿ2n)) o ⁿ12n
i.e.    J u,dz = An (uo + 4ⁿn + ^u2n),
0
or, on changing the variable to x, r2n
J uxdx = An (1d₀ + 4nn + ^u2n)• 0 This principle is of universal application and any formula of
approximate integration can immediately be transferred from unit intervals to nthly intervals and vice versa.
For example, the formula
uxdx = 22 (u_l + u₄ + u₆ + u₉)
(20
becomes    ^fuxdx = 5 (u₂ + U₉+ 11₁₂+ u₁₉)
0
on doubling the interval;
iOn
and    Jo u
dx = 2n (u+ ^u4+ ^u-+ ^u9)
when the interval is n.
F    X9

r10 0

290    INTEGRAL CALCULUS 4. Change of origin.
Consider again Simpson^'s rule:
u_sdx = (uo + 4U₁ + 11₂). 0
We have obtained a formula which gives in terms of u₀, u_i and u₂ the area of the curve cut off by the Y ordinates x = o and x = 2. If we    ^Y'

change the origin so that the point (1, u_i) becomes the point (o, vo), as in Fig. 38, we shall have a formula of integration between the limits -1

and I in terms of v_₁, vo and v_i.    v    u
The formula will be otherwise un- 110    zr, °    2 v,
altered and we shall have found an approximation to the same area with reference to a new y-axis 0₁ Y_i .
In its new form the formula be- 0    0,    %
comes    Fig. 38,
Jvdx = (^v-1 + 4^vo + ^vi), -i or, changing back to u^'s,
~ u_xdx = (^u-i+4¹¹0+¹¹1)• __i

Now when it is desired to obtain an approximate integration formula it is evident that integration between — r and I (or between — n and n) involves much less labour than integration between o and 2 (or between o and 2n). Thus, if

	uy⁼a+bx+cx²,
then	¹¹0^=a,
	11_₁=a—b+c,
and	U₁=a+b+c,
and	(u__i + 4^uo + ¹¹1) =	{4^uo + (^u-1 + ^u1)}
	=	(4^a + za + zc)

= s (6a + 2C).

CHANGE OF ORIGIN 291
Also _J¹₁(a+bx+cx²)dx= [ax + 2bx2+ ^3cx3J¹1 = 2a + sc,
which proves the approximation.

Again, since in this form the coefficients of odd powers of x disappear in the definite integral and also in the paired terms u_t and u_, , we could equally well employ as our assumed expansion for u_x the third difference function a + bx + cx² + dx³.

We should have
4uo+(u_₁+III) =4a+a—b+c—d+ a+b+c+d=6a+zc as before, and

[ax + 2bx² + 3cx³ + 4dx⁴]¹ _ = 2a + 3c. 1

This leads to the important fact that Simpson^'s rule is true to one more order of differences than was originally assumed, i.e. it is true to third differences.

For the above reasons, namely

(i)for simplicity in working;

(ii)to enable us to find the true order of differences to which the approximation is correct;

it is advisable to integrate between —

and

in preference to any other limits.

This can always be done by a suitable change of limits; for example, integration between _zand I can be simplified to integrabetween say z and 4 by increasing the interval from to 2. Then, by changing the origin, a formula can be obtained for intebetween — i and 1. To express the formula in the required form the process can then be reversed.

Example

'2'5

Show that }

= 2¢

₅

+ 2211

₂

₃

)

approximately.

    -2₅    ₁₅v,. dx = - / u~ dx (if we double the interval).
    .1'5    ².3
This will produce a formula in u₂; u₄; us.
In the first place we will obtain a formula for j u_x dx. 1
19-2

292    INTEGRAL CALCULUS
This formula will have to be of the form mu_₂ + nuo + pu₂, since we have moved our origin to the point where x was originally 4.

Let	u_x = a + bx + cx² + dx³.
Then	u₀ = a, u_₂=a—2b+4c— 8d,
and so that	u₂=a+zb+4c+8d,

(u_₂ + 2211₆ + u₂) = _ⁱ} (22a + 2a + 8c) = a + lc. f1 u_xdx = [ax + 2bx² + 3cx³ + edx⁴J¹ ₁= za +ic. 24 (u_₂ + 22u₀ + ^u2) = 2 I_ 1 u_x dx.
Changing the origin back again to the original origin:
5
24 (11₂ + 2211₄+ u₆) = 2 /3 u_ydx.
2^.5
(711+22712+V3)= [ uxdx. 24 . 1^.5

Moreover, although only three values of v_x have been given (so that only second differences are known), the formula is true to at least third differences.

5. Some well-known approximate integration formulae. (1) Extension of Simpson's rule.
2 We have    f
o u_xdx = (uo + 4u₁ + u₂). J
Therefore by changing the origin

uxdx= ₃(u₂+411₃+11₄).
16
Similarly    u_xdx = (u₄ + 4u₅ + 11₆). 4
.................................... In general, by addition,
r2n
11xdx = (u₀ + 411₁ + 2112 + 4113 + 2114 + ... + ^u2n)^. .o
This is a formula for approximate integration used extensively in practical mathematics. In engineering problems, where the
j 4

THE ^"THREE-EIGHTHS^" RULE    ²93

curve to be integrated has actually been sketched as the result of experiments, a small unit may be chosen and a large number of ordinates drawn. In that event the extended formula can be apThis will in general give better results than the single formula

₁2n
n

o uxdx = 3 (11₀ + 4U,, + ^u2n)

It should be noted that Simpson^'s rule, applied many times (as above), does not assume that a smooth curve can be drawn between all the points u₀, u₁, ... u_>,,. The method of obtaining the formula has in fact been to draw a number of disjointed curves between u₀, u_1iu₂; u₂, u₃, u₄; u₂,,_i and the curve passing through three points such as u_o, u₁, u₂will not as a rule pass through the next series
u₃and u₄.
(2) The "three-eighths^" rule.
The following symmetrical formula can be derived by working on the above lines, when four consecutive points are given.

If the four points were u₀, u_1iu₂, u₃ we should obviously inbetween o and 3. Change the unit and origin so that we pass firstly to the limits o and 6 and secondly to the limits — 3 and 3
3
We require ^j'-u_xdx, given u_₃ _iu_₁, u₁, u₃.
3

Let    ux = a + bx + cx² + dx³,
so that    u_₁+21₁⁼2(a+c)
and    u_₃ + u₃ = 2 (a + 9c).
3
If    Jr ₃u_xdx = m (u_₁ + u₁) + n (u_₃ + u₃), we have
_C3

ax + zbx² + cx³ + idx⁴ = m (^u_1 + u₁) + n (u_₃ + u₃);

-3
i.e.    6a + 18c = m (u_₁ + u₁) + n (21_₃ + u₃)
= m (2a + 2c) + n (2a ± 18c)
= a (2m + 2u) + c (zm + 18n), from which, easily, m = 18/8 and n = 6/8.

²94    INTEGRAL CALCULUS
J - uxdx = {6 (u_1 + u1) + 2 (u-3 + i.e.    3
i.e.    JO
uxdx = $ {6 (u₂ + u₄) + 2 (up + u₆)},
0
or    u_xdx = s {3 (u₁ + u₂) + (uo + u₃)} = $ (u0+3¹11+3u₂+u3). (3) Weddle's rule.

Suppose that seven equidistant ordinates are given. In order to reduce the algebra to a minimum, we integrate between limits — 3 and + 3, and write the assumed formula for u_x thus :

_ux=a+bx+cx²₊dx(x^2—I)₊ex²(x^2—I)₊fx(x^2—I)(x²-4) 2    6    24    120
₊^gx2(^x2—I)(^x2—4)
720
This is Stirling^'s formula with the constants a, b, c, ... g re-placing the differences of u_x.
Then
Judx=        bx² cx³ d x⁴ x²e x⁵ x³l
Cax+ z + 6 +6(4—2)+24(5 3/
    f x^s5x⁴ 4x²₁    g      (x7 ₅x5    x3;1 3
+ 120 \ 6    q⁺2/⁺ 720 7    5    3 J 3 g = 6a + 54 (₆72 + 540) + ⁴⁸⁶(120 720) + 280^g ^=6a+9^c+ ^e+ 4^'6- g^.Replacing a, c, e, g by the differences of u_x in Stirling's formula, we have
3
(_ 3 u_xdx=6a+ 9c+ o^e+ i⁴hu^g
= 6u₀ + 9^{A2 11_}1 + 0⁴^u-2 + 14., ^{A6 71_}3

If now A2u_1, _{04u_2}^and_{A'u_3}are expressed in terms of ^u_3, ^u-2, ... U2, 11_3iwe shall obtain a formula correct to seventh differences. It will be found, however, that the coefficients of the terms are large. This is due to the awkward fraction 1₄¹₀multi-

WEDDLE^'S RULE ²95

plying „⁶11_₃. As sixth differences are usually small, the error involved in replacing _{r`'o-06u_3}^by_{4?;06u_3,}i.e. by -₁'',0⁶u__3i will, in general, be negligible. By this substitution the coefficients in the final formula will be much simplified. The modified formula will involve an error of ₁1,Q⁶u_₃ and will therefore be correct to fifth differences only.

The terms in the expression thus adjusted, namely,
6uo + 90²11_₁ + 'It”; ^A4u_2 + A^4u-3, may be collected as shown in the following table:

	^u_3	^u_2	^u_1	u0 ^u,		^u2	U₃
6u,				6•o
9^J21L-1			9'0	—i8•o	9^.0
3'3-^1'U—2		3'3	— 13^.2	19^.8	— ¹3'²	3'3
3~^eu~	•3	-1^.8	4^.5	— 6•o	4^.3	-1^.8	'3
Total ...	•3	1^.5	•3	1^.8	•3	1^.5	•3

3
Therefore _{J_3}u_xdx is approximately equal to
•311_₃ + I^.511_₂ + _/•311_1 + I 811₀ + •3111 + 1'5112 + •3113
11^3-[(u-3 + 113) + 5 (¹¹-2 + u2) + (^u-1 + u_l) + ^6u0~
.⁶or    I₀usdx =    [(ⁿo + ¹¹6) + 5 (u_l + 11₅) + (112 + u₄) + 6u₃].

This is WVeddle's rule.
(4) Hardy's formulae.

Certain approximate formulae due to G. F. Hardy have been used extensively in actuarial work.
Let    u_x= a + bx + cx² + dx³ + ex⁴ + fx⁵
3
so that    ( u_xdx = 6a + 18c + !Is. -3

11₀= a; 11_₂ + 11₂= 2a + 8c + 32e; 11_₃ + 11₃= 2a + 18c + 162e.

Solving for a, c and e, and substituting:
((³

J uxdx = 2^.211₀ + I^.62 (u-₂ + 11₂) + •28 (u_₃ + u₃) -3

296    INTEGRAL CALCULUS (6
or    I u_xdx = 2'2u₃ + I^.6z (u₁ + u₅) + •28 (u₀ + ZZ₆) o
= '28 (u₀ + u₆) + I^.62 (u₁ + u₅) + 2^.2u₃ ,

which is Hardy^'s ^"formula (37)."
If the interval of differencing be n, this becomes
r6n    l
J u_xdx =n{'28(u₀+n_6n)+I^.62(un+ u5n)+2.2%n). 0 Similarly
12n
uxdx = n {^.28 (u₆n ± ^u12n) + I^.62 (u7n + ^u11,,) + 2'2ugn},
on
and so on. Since
6n    12,    18n
J u_xdx = J u_xdx + J uxdx + J u_xdx + ... , 0    0    6n    12n Jusdx = n {28 (Ito +^2u6n+^2Z112n+ •••)
+ 1'62 (u_n + u₅n + ^u7n + ...) + 2'2 (^u3n + uyn + ...)}. This is Hardy's ^"formula (38)."

If we are dealing with functions derived from a mortality table, and we choose n so that 7n falls just within or just without the limits of the table—so that 7n is in fact about too—we can write the above formula thus :

_J⁰u_xdx = n (^.28% + I^.62u_n + z•2u₃n + I^.62u₅n + •56up_n + I^.62u₇„).
0
In this form the formula is known as Hardy^'s " formula (39a)." (See G. King, Life Contingencies, p. 488.)

Note. Since all these approximations give in effect the area of the curve bounded by u_x, the limiting ordinates and the x-axis, they are often termed ^"quadrature formulae," the word ^"quadrabeing defined as the exact or approximate calculation of the area of the square equal in area to that of the given figure.
6. Practical applications of the formulae.

Integration formulae of the above type may be used to advantage in obtaining approximations to the values of certain complicated forms of integral which occur in the theory of life contingencies.

PRACTICAL APPLICATIONS ²97

Where the functions involved are such that the summation extends to the end of the life table, it is customary to calculate the values of the functions by Hardy^'s formula (39a). If, on the other hand, the upper limit is at an age-point short of the limiting age of the life table, any of the simpler integration formulae can be employed to advantage. (See, for example, Spurgeon, Life Contingencies, pp. 262, 263.)
It is evident that we can neither integrate the function u_x nor interbetween given values of u_xif we are absolutely without informaregarding the value of the function between the ordinates. For example, in a mortality table we may integrate between, say, 1₃₃and ¹31, assuming first differences constant, for we know that in general the decrements between ages 30 and 31 are small and may be fairly conas being evenly distributed over the interval. We could not reasonably adopt the same assumption for interpolation between to and 1₁, since the deaths of infants in the first year of age are not evenly distributed over the year. Again, a reliable estimate of the population of a seaside resort in June of any calendar year would not be obtained by a first difference interpolation between the population figures for January and December of the same year. In the two latter illustrations further information would be necessary before we could proceed to interpolation or integration.
In applying the formulae of approximate integration to the solution of a problem it is therefore essential that we have sufficient knowledge of the function to justify our assumptions regarding the nature of the curve passing through the given points.
When the function to be integrated is one in which the slope of the curve is gradual between the limits of the integration, almost any of the above formulae will give satisfactory results. If the values of the function are changing rapidly between the limits there may be considerable differences between the approximate results. The following examples will show this clearly.
Example 2.

The formula for a continuous annuity-certain is _    (
a,-.1 = m
I    e-Ft dt,
.o
where 8 is the force of interest corresponding to the rate of interest i at which the annuity is valued.

298    INTEGRAL CALCULUS
By means of an approximate integration formula find the value of a_E at 4 per cent.
Since we have to evaluate a,, we shall require some or all of the values
of u_t from t = o to t = 6, where u_t = e^—81.
By a well-known formula in the theory of interest,
e—at = v^t, where v = (1 + 0^—1.
Values of v^t at 4 per cent. are available from any tables of interest. We have, therefore:
u° = ^e°    = ¹,
11₁= e_s = v = '9⁶¹54,

¹¹2 = e—21 = ^v2= ^'9²45⁶, u, = e^—38= v³ = '88900, u, = e—41 = v⁴ = '8548o, us = e—55 = v⁵ = '⁸²¹93, 116 = e—63 = vs = ^'79⁰3^1.

The following results will be obtained:

(i)Simpson's rule applied three times :

3 (

0+4

¹¹

2112+4113+

₄

+4115 +

= 5^.34630.

(ii)

The three-eighths rule (with interval of differencing

2):

8 {Z (Tt_° + 311₂ + 3U₄ + u₆)}= 5'34629.

(iii)

Weddle's rule:

(u0

511

₁

₂

611

₃

₄

+ 5U

₅

_°

)

= 5'34631.

The result obtained by integrating

^5t

will be L

^—

est]s or

v°-v6

₈

and since the value of 8 at 4 per cent. as given by an interest table is '039221, the integral

_ 1 —

'79031 _ '20969

'039221•039221 = 5'34

⁶

The three approximate formulae all give results differing from this in the fourth or fifth place of decimals. The excellence of the approximais due to the fact that the slope of the curvey = u_t = e -^8tis gradual between the limits t = o and t = 6.

THE EULER-MACLAURIN EXPANSION ²99
Example 3.
Evaluate '^s ^dx (i) exactly, (ii) by an approximate integration
o (1+4x)2
⁽¹⁾.10 (1+4x)^2—L (1+4x)]0 = 4(I—21) =24. (ii) By Simpson's rule applied three times, we obtain ²3745⁸;

„ the three-eighths rule „	²'43²49
„ Weddle's rule „	2^.36382.

It will be seen that the differences between these results and the accurate result are considerable. None of the values is correct to more than the first place of decimals, and, moreover, the errors are not all in the same direction. The function y = (1 + 4x)-² is changing rapidly over the given interval, and it is not to be expected that any of the above formulae will give a near approximation.

7. The Euler-Maclaurin expansion.

We have considered quadrature formulae in which the number of ordinates used is a simple multiple of the first two or three natural numbers, 2, 3, 4. We will now proceed to investigate the problem of approximate integration more generally.
The basic formula for expressing a definite integral in terms of given ordinates is the Euler-Maclaurin expansion. This expansion is in effect of similar form to the expressions already obtained, a greater degree of accuracy being ensured by the addition of funcnot of the ordinates themselves, but of derivatives of certain of the ordinates.

The formula may be derived by the expansion of operators, thus : x=n-1
f(x)=f(o)+f(1)+...+f(n—i)
x=0

= F (n) — F (o), where f (x) is AF (x).

Since f (x) = AF (x),
F (x) = A-¹1^-(x)
= (e^D — 1) ¹f (x) (since 1 + A = e^D. Chap. xi, para. 13)
D2 D³

= [(I + D + 2~+-+...) — I f(^x) 3 •
formula.

300    INTEGRAL CALCULUS
I    D    D2    ...1
^I+₂,-₃,+    f(^x) D²=D-¹ 1 ~ _–d
... f(^x)
720 D D³
    ^D-1— _+ 12 —    {
720 ... J (^x)
D³
= D-1 f (^x) – if (^x) + 12 f (^x) – 720 f (.x) ...
i = Jf(x) dx – f (^x) + i z df dx~⁾ 72ⁱ 0 d (x) ... . Between limits o and n, we have therefore
F (n) — F (^o) = Jf(x) dx — {f (ⁿ) — f (^o)} + { f' (ⁿ) — f' (^o)}
— h If"' ⁽ⁿ⁾-f'" (o)} ... .
For F (n) — F (o) we may write
x=n-1
f (x) ^orf (⁰) +f (^I) +f (²) + ... + f (n — 1).
Jf(x)f(0)+f(1)+f(2)f(n_1){f(n)_f(0)}
— ~'ⁱ {f^' (ⁿ) ^—f^'(0)} + -7-i u If „, (n) -f'” (0)1...
=2f(o)+f(I)+f(2)+...+f(n—I) +f(n)
— , If ' (n) -f' (0)} +    If"' ⁽ⁿ⁾-f "' (o)} ... .
This is a simple form of the Euler-Maclaurin expansion.
A more general form can be obtained by changing the origin to the point a and the unit of measurement to r, thus :
Ja+nr r a f(^x)^dx=f(a)+f(a+r)+f(a-; 2r)+... +f (a+n — ^ir)+2f(^a+^nr)^—i2{f'(a+nr)—f'(a)}
r³
    +    _o { _f,,, (a + nr) — f,,, (a)} ....
It will be noted that, since
s
    x²x⁴ = I — Zx+B₁₂~^—B₂ , +B3hj —... 4•
x
e^x—
(Chap. xi, p. 185),

    ILLUSTRATIVE EXAMPLES    30I

we can express the coefficients in (e^D — 1)-¹ in terms of Bernouilli's numbers. As, however, the resulting approximation formula is to be used for numerical computation, it is of advantage to give the coefficients their actual numerical values.

8. The following examples are illustrative of the use of the formula.
Example 4.

Evaluate I _oII dx x to five places of decimals.
Choose a convenient unit, say o•I. Then in the Euler-Maclaurin
expansion we have
a = o, n = Io, r = o•I, and u_x =
I + x^'
du_x d³ u 6
dx - - (I + x)^2' dx³ = — (1 + x)^4'
I 11 dx I I I I I I I I ^0'I-.^IO ^I+x 2.I+I—I+I•22+~ ₃+...+ I9⁺2 2

^—*^°O2 -22+j2J+0200I[— 6 ⁺6]^...

50000 — 120 ' 4 + 72000016

_'90909

⁸3333
.₇6₉₂₃71429

66667 •62500 •58824 5555⁶
52632 •25000 ⁶'93773

= 6

93773 — •00625 +

'00001

⁶

'93

rIdx = •69315 to five places of decimals. 01+x

This agrees with the true result (log

2) to the required degree of accuracy.

302    INTEGRAL CALCULUS
Example 5.
Find the sum of the fourth powers of the first n natural numbers by means of the Euler-Maclaurin formula.
In the formula
1 fa+nr
u_x dx = s^ua + ^ua+r + ^ua+2r + ^rr
+ ?ua+nr    ₁₂(^ua+nr    ^11'a) + 720 (^1la+nr — 11 a) • • • ,
put a =_r 0, r = I, u_x = x⁴; then
1 ~t x⁴dx = 1⁴ + 2⁴ + 3⁴ + ...
I
+ (n — 1)⁴ + ~n⁴ — ₁- L4x³]_x=n+ 720 [4.3.2x]_xan
since higher differential coefficients of u_xwill be zero.
x,5 n r=n    I    I
I.e.    [ ] = E r⁴ — ¹-n⁴ — 4n³ +    24n.
5 o r=1    12    720

r=n I I I I
.. E Y⁴=-n⁵+ - n⁴+-n³--n.
r=1 5 ²3 3⁰

By proceeding on the above lines it can easily be seen that the general formula for the sum of the pth powers of the first n natural numbers is
r=n
E Y^P= n²'⁺¹+2n^P+ pn^p-1— I ^p(^p—I)(^p—2)^nV-3r=1 p + 12 720
₊^p(P — 1)(^p—2)(^p -3)(P—4)np_5
9. Lubbock's formula.

The previous formulae have been developed for the purpose of relating a definite integral to the sum of a number of ordinates at finite distances apart. We have, in effect, obtained approximate formulae for the value of

Lt h (ua + ^ua+h + ^ua+2h + • • • + ub). h-~o

In addition to formulae of this type we can find expressions which enable us to find the value of the sum of a number of orat finite distances apart in terms of the sum of ordinates at greater or less intervals apart. Thus, if for the curve y = u_xthere are nt unit intervals, so that we have Eu_x from x = o to x = n1 — 1,

30,240

LUBBOCK^'S FORMULA 303
we may develop a relationship between this sum and the sum when the intervals are, say, h, where nh = 1.

In place of the sum h (u_a + ^ua+h + ^ua+2h + ... + us), we may consider, without loss of generality, the simpler form when the series commences with u₀.
Let there be m ordinates at unit distance apart, and let each of these unit distances be divided into n equal parts, so that the new ordinates are

U0, _Ul'n,U₂/n,    . ur/n, ... ^um—l/n^.
Replace 1/n by h for convenience in working.
Then h (uo + u_h+ ^Ugh + ... + ^Um—h) =h(I+E^h+E^2h+...+E^m-^h)u_o
E^m — I h
E_{h_}_Iu_o= h _Ehl i (um — ^u0)
1
^h₁(I + 0)^h — _I(^um — u₀)
= h (hA + h₂A² + 11₃:,³ + ...)^—1(^um — ^Z10)
h _        h' —
= h.h^—10^—1 (1 — 2 I 0 + I2 ¹02 h² -1 3 `
I    -2'4' Z ... I (um — ^uo)
2 — //
= 0^—1(um — uo) — 2 (um — u₀) + 12 I 0 (um — u0)
_Z2
^—I    — -    A²(um^—ZLo)...
²4
_ (uo + u_l + u₂ + ... + ^Um^—1) — ^h₂I (nm — ^uo)
+ ^h212 I    h
(,u_m — Auo) — ~4 (A² u_m — A² uo) ... ,
24
since A(uo+ul+u2+... +^um^—1) _ (11l+u2+... +^um^—1+^um) —(u0+u1+... +^um^—1) =Um^—U₀.

^I/n2— ^I—    (A² um — A²u₀)

₂₄3⁰4
But h = 1/n.

INTEGRAL CALCULUS

n(uo+ul'n+u2,n+... +um-i/n)=(uo+u1+u2+...+um-1)
1/n2    ^I(um — u₀) + 1/n
—        ₁₂(Sum — ^Duo)
or
uo + U1!n + u₂^,₁, ± ... + ^um^–1/n = n (uo + u₁+ ... + ^um^–1)
I —n (um^— u0)+I —n2 (Aum—0110)—I— n2 (02 _umA2 uo)... 2    I2n    24n
n _ I    2 — /
= n (uo + u_i+ ... + ^um^–1) + 2 (um — u₀) —     (
ⁿI2n I ^Um — Duo)
This is Lubbock^'s formula.
The coefficients of higher differences than the second are cumber-some, the terms in A³u and A⁴u being respectively
(n² — I) (1911² — I) (^L.3um — ^03uo)
720 n³and    (n2 48onⁿ² — I) (04 u_m— 0⁴1(₀).
If the interval of differencing be originally n instead of unity, Lubbock^'s formula becomes, on changing the unit, u₀+u₁+11₂+ ... +^umn–1
=n(u₀ +u +u +...+u(m–1)n)+-     ₂(umn — ^U0)
n    2n
^n2—I _(DumnDuo) + n^2—I Z2umn— L121(o).... 12n    _24n10. Woolhouse's formula.

Although Lubbock^'s formula has the advantage that it may be used when the function is not capable of expression as a matheexpansion—as for example when the data are based on a

n² — I
+ - (A²um — L1²uo) .... 2411

WOOLHOUSE^'S FORMULA 305
mortality table—there are disadvantages in adopting the formula. In the first place, if it is necessary to proceed further than second differences the calculations are heavy, and secondly it may happen that the differences converge slowly, so that if we stop at second or third differences we are likely to obtain a result differing considerably from the true value of the function. Where a mathematical function is available an alternative summation formula can be adopted in which differential coefficients replace the finite difin Lubbock's formula.
The formula involving differential coefficients is due to Wool-house, and may be developed directly from the Euler-Maclaurin expansion.

The Euler-Maclaurin expansion is
111
uxdx = 0
2u₀ + n₁ + 11₂ + ... + ^u,n—1 +Zuni    IZ (um^' — u₀^')
¹r„    ,,,
+720(unZ— uo) ....
If the interval of differencing be I/n, the formula becomes
    rm    _In J uxdx = (2110 + n_l n + u₂,_n+ ... + Zum)    (um' — u₀^')
    o    I2n 0
"'    rrr
+72on(11 —uo )....
If, however, we multiply both sides of the first expression by n we have
n
n    u_xdx=n( uo+u1+u2+••.—zu„r)(um'—110')
    j 0    ₁₂n
+    (^um"' —a^„') ....
₇₂₀0
Equating the two values of n Jndx:
yu0+ul,n+u2+...+ urn^—_IX11(ur^'—u₀^')+     2 (urn"'�,rr)...
720n

=n(?u0+u1+112+...+gum)— n (it '�')+'~20(um,,, ^—u₀^rrr)...
12 F 20

306    INTEGRAL CALCULUS
or
I uo+ul/n+u2/„+...+um— (uo+um) ^—12n^(um'—nog)+ ^I(u nr — uo^ur)
720113    ...
m
n n(uo+ul+u2+...+um)—2(uo+um) _ ₁₂(n_n,^'no^r)
+ ..n    u rrr — u_o''')
~₂O(m
Re-arranging : uo+uln+u2n+...+um=n(uo+n_l+u₂+...+um)
1Z —    (Ito + um) — 11²-

-(um' _ _uo,₎+ 71⁴ — I (_ulna' — norm) ... .

I2n

720/1

If the unit of measurement be changed to

we have

uo+u,+u

₂

+...+umn

⁼

n(uo+un+u2n±...+umn)

—

(uo

mn) —

—

(u, ,, — uo

)

+ - -

(um,nrrr —

u0ur)

... ,

720

the usual form of Woolhouse's formula.

It should be noted that, by replacing the derivatives of u by their values in finite differences, Lubbock^'s formula can be obdirectly from the above.

In applying these formulae to certain actuarial functions the values of u,

dx '

(1,0

at the end of the mortality table will disappear. Woolhouse's

formula may then be written as

n—n

^2—

^4—

(uo

+ ...) —

+ -

uo' —

^rrr

....

2n1271

72011

This is a convenient form for expressing the value of a benefit paid

nthly

intervals in terms of the values at intervals of a year.

11. Other formulae for approximate integration.

It will have been observed that in the formulae of the type of Simpson^'s, Weddle^'s, etc., the function of the u^'s is symmetrical about the central value. If, however, a fixed number of ordinates be given and it is desired to obtain an approximation to the area

APPROXIMATE INTEGRATION FORMULAE 307

of a curve in terms of these ordinates, the resulting form will not necessarily be symmetrical. Again, the formula for the area of the curve may be related not only to ordinates falling within the area to be measured, but to ordinates outside the area. It may be noted therefore that although standard formulae are avail-able, it is not difficult to devise approximations to fit the particular problems under investigation.

Example 6.
If u_x is of the form a + bx + cx², find a convenient formula for
t
f u_xdx in terms of It_o, u₁ and u₂.
The interpolation formula which involves the terms u_o, u_l and u₂ is
u_x = 110 + x,11₀ + 2x (x — I) L.² 11₀.
Jo^t
r
u_xdx = Lxuo + ix²Au₀ + (bx³ — x2) ^A2uol₀
=uo+?Atl_i— 11,'02uo J

=110+ (^u1^—u0)^—3 (^u2^—2¹¹1+110)
1

= i°-^lto + ~,_ _ltl— 1. ¹¹2^,
The required formula is therefore
t
1₀u_xdx=,x;(511₀+811₁u,).
It should be noted that (i) the expression is unsymmetrical in the u's, and (ii) the ordinate u₂falls without the area to be integrated.

It is obviously impossible to quote all the formulae that are in current use. Further examples and illustrations of various approximation integration formulae will be found in the following sources : Whittaker and Robinson^'s Calculus of Observations, chapter vli; C. H. Wickens, J.I.A. vol. LIv, pp. 209—13 ; A. E. King, T.F.A. vol. Ix, pp. 218—31.

12. Alternative methods of proof of the formulae.

It has been stated above that the Euler-Maclaurin expansion is the basic quadrature formula. It will be instructive to develop Simpson's formula from this expansion.

J 2mn n2
uxdx⁼n(zu0+un+ll2n + . ••+u(2m—1)n +2u2mn) ₁₂(^u'2mn — ^u0^') 0
approximately.
20-2

308 INTEGRAL CALCULUS
Writing zn for n, but preserving the same range o to 21nn
j 2nzn 0
4ⁿ², _I2(^u2mn — u₀).
Subtracting this from four times the first:
r 2)nn

3 J u,dx= 2n(zn0+tun+u2n+2113n+... +^u(2m—2)n 0

i.e.
f 2nen rz
J 0 uxdx = ₃(u₀ + 4u_n + 211_2n + 4¹¹3n + ... + 211₍₂,n_₂) n
+ ^2u(2m^—1)n + 2^u2mn),
which is Simpson^'s formula.

Again, we have shown in para. 2 (p. 288) that we may adopt the expression for u_xin terms of uo and differences of u₀ as the assumed form of function. Since Lagrange's formula is based on the same assumption it is evident that we could obtain any approximate formula by the use of the Lagrange formula.

For example, given uo, u, , u₂,
(x — I) (x — 2)    x (x — 2)    x (x — I)
^Zl.c = u₀(— I) (— 2) + Zl_l     1.(- I) + u₂     2. I
= 2u₀ (x² — 3X + 2) + u₁(2x — x2) + 2₁₆₂(x2 — x),
so that
J 2    (x3    ^~2_x3x3 x212
o^uxdx=`
[^u06 -34      +x)+u1(x2— 3)+u2(6 41 0
= 3 (uo + 4u₁ + u₂), which is Simpson^'s rule. If we integrate between o and I we shall obtain
)I
u_xdx=(5uo+811₁—u..),
the formula given in Example 6.
u_xdx = 211 (2u₀ + ^u2n + Zl_ln + •.. + ¹¹(2m—2)n + l^u2mn)
+ 4^u(2,n^—I)n + ^lmn),

EXAMPLES    309
EXAMPLES 17 1. Prove that
i
e ¹⁰dx.
1

3.Show that [ u_x dx =(fiu1 + 8u_° — u__₁) approximately. Find the approximate mileage travelled between 12.0 and 12.30 by use of the above formula, from the following:

TimeSpeed (m.p.h.)

11.50	24'2
12.0	35^.0
12.10	41'3
12.20	42^.8
12.30	39^.2

4. Prove that, if u_x is a rational integral function of x, then ~e^a u_x dx = ae^a (I — aD + a²D² — a³ D³+ ...) u_x,
where D = d dx^'

Show that the area of a curve, divided into n parts by n + 1 equiordinates u_°, u_t, .., u,,, is given approximately by the series

ⁿ²

Au,

+ {

ⁿ³

—

ⁿ

^{02 u0}

+ {

ⁿ

—

ⁱ

^A3

°'

1.24

1.2.3

to n

1 terms.

Between the limits x = o and x = n the functions u_x and dux/dx are continuously increasing.

Show that /

dx is

less than

_°

ⁿ

+gun.

Obtain the approximate formula

(

-1)

(u3

₇

showing up to what order of differences it holds.

₍

dx = I° (u

2a—x)

^dx

and illustrate the result geometrically.

= a + bx + cx

prove that

₂

_°

—

₂

₄

and hence find an approximate value for

310 INTEGRAL CALCULUS
5
8. Assuming u_x to be of the fourth degree in x, express [ u_xdx in .o
terms of u₀, u_1iu₂,u₃andu₄.

A plane area is bounded by a curve, the axis of x, and two ordinates. The area is divided into five figures by equidistant ordinates 2 inches apart, the lengths of the ordinates being 21^.65, 21^.04, 20^.35, 19^.61, 18^.75 and 17^.8o inches respectively. Apply the method of integration to obtain an approximate value of the area.

Io.

Prove the approximate formula

(ii) Simpson

s rule :

_Io

dx =

(11

₀

411

₁

₂

applied three times.

Iz.

Which of the two following formulae would you expect to give

the better approximation

for'

dx?

^y{5 (^u0 + ^u4) + 4 (u_I + 11₃) + 18u₂},

^(b){2 (u0 + u₂ + ¹¹4) + 12 (u_I + ^u3)}

13.Prove that

n {~u

₀

+ }

(

¹¹

,1 —

zn +

3rz)},

and hence find /

given the following table: .o

27,65031,25235,

5439,3

⁶⁸

14.Prove that

(x) dx{

f (— )) + f O} +

— Ala)}

approximately. Hence find

11f

(x) dx,

when

(o)

= 105,

(I)

212,

(2)

= 421, f (3) = 749 and

(4) = 1050.

2'5

(1l

₄

+ 1[g + 11

₉

and show that the formula involves a small second difference error.

-s

11.

Find the value of

_to

x)-

dx.

Obtain approximations to the value by applying (i) Weddle's rule:

dx=0'3(aio+5u

₂

+611

₃

+11

₄

+511

₅

+us),

EXAMPLES 31I

15.If u_s either increases continually or decreases continually as x

increases, prove

that'

^~L

differs from

^nEI

by less than the differ-

.I1

ence between u

₁

and u,,.

Prove that the difference between log

and

1 + + + ... + 1

(n -<

however great

may be.

16.If third differences are constant, prove that

dx =. ?} (u_

2311

₁

2311 +

u5).

Adapt this formula to find the approximate value of log

from the

integral 2a dx ' _x

17.

Prove that, approximately,

/ u

= 0

111_

₃

+ 4'411_

₂

5tt_

₁

+Ono

+ ¢•q.u

₂

+ 1

111

₃

.-3

18.f (x) is a rational integral function of the fifth degree in x. Prove that

f¹f(x)dx= f(o)+ {f('Vo•6)+f(--Vo•6)}.

19.

Use Simpson's rule to prove that log, 7 is approximately

1'95.

2o. Apply

the Euler-Maclaurin formula to find a formula for the sum of the fifth powers of the first

natural numbers.

21.Obtain Shovelton^'s integration formula:

-Io

dx =

i4o- {

⁸

(

io) + 35

(

₃

₂

₄

+ U

₉

)

+ 36u

₅}

22.By means of Hardy^'s formula

dx = •28uo +

62u

+ 2•zu

₃

+ 1

6zu

₅

+ •

z8u

, .o

calculate the value of

)

^-1

₀

correct to four places of decimals.

23.If f (x) be a function of the third degree in x, and if

ⁱ

J-3t

f (x) dx, tto =

ttt

(

)

^dx

^3t

(

)

^dx

show that

f (o) = - {u

312 INTEGRAL CALCULUS

24. AB is the base of a semicircle, centre 0 and radius unity. The points P and Q bisect OA and OB respectively. The area between the
semicircle, the base PQ and the ordinates at P and Q is 6 + -³. Use
4
Weddle's rule that
.f(x)dx=0'3{f(o)+5f(I)+f(2)+6.f(3)+f(4)+5f(5)+f(6)} to find an approximate value of 7T to three places of decimals.

25.The following values of zz . are given:

2345

⁶

146

•161•176

•190

•204

•217

•230

Use an approximate integration formula to find the value off

dx.

It is found that, for the values given, y = log₁₀ ('05X + P4) fits the data. Verify that this is so by integrating log₁₀ (•05.x + 1^.4) between the limits o and 6. (loglo ^e= .4343; log₁₀ P7 = •2304; Iog₁₀ 1^.4 = •1461.)

26.

is a function whose fifth differences are constant,

ⁱ¹

can be expressed in the form

pu_„ + qu₀ + pu_t.

Find the values of

and a.

Use this formula, after making the necessary changes in the origin and scale, to find the value of log_e 2 to four places of decimals from the equation

l I

+ x

dx =log

27.Prove that, if a= o,

ra+r

= Iu

₀

+ u

₂

+ ... +

_{r_I}

u,. —

(

^Du

r_I —

^Du

_r-2

+u0)(.

^3u

r_3 —

^Q3u

i u u (

^A'u

r_4

A4u0) ... .

If u

be the function

+ x

)

^—I

find an approximate value for

28.Obtain an approximate formula for' u_xdx in the form

a (u_₂ + u₂) + b (u_₃ + ^1[3)

and find the values of a and b.