Anyone who has to use analysis as an instrument for investigating physical or technical phenomena is faced by the question whether and how the theory can be adapted to yield useful practical methods for actual numerical calculations. Yet, even from the point of view of the theoretical person, who only desires to recognize the link between natural phenomena and not to conquer them, these questions are of no trifling interest. For a systematic treatment of numerical methods, we refer the reader to special textbooks on the subject.* Here we can only discuss some particularly important points which are more or less closely linked to the preceding ideas. We wish to direct special attention to the fundamental fact that the meaning of an approximate calculation is not precise unless it is accompanied by a definite knowledge of the degree of accuracy attained.
* For example, Whittaker and Robinson, "The Calculus of Observations" (Blakie & Son, Ltd., 1929).
We have seen that even relatively simple functions cannot be integrated in terms of elementary functions and that it is futile to make this unattainable goal the aim of the integral calculus. On the other hand, the definite integral of a continuous function does exist and this fact raises the problem of finding methods for calculating it numerically. Here we shall discuss the simplest and most obvious of these methods with the aid of geometrical intuition and thenconsider error estimation.
Our objective is to calculate the integral
where a < b. We imagine the interval of integration to have been subdivided into n equal parts of length h = (b - a)/n and denote the points of subdivision by x1=a, x2=a+h, ··· , xn=b, the values of the function at the points of subdivision by f0, f1, ··· , fn and, similarly, the values of the function at the midpoints of the intervals by f1/2, f3/2, ··· , f(2n-1)/2. We interpret our integral as an area and cut up the region under the curve into strips of width h in the usual manner. We must now obtain an approximation for each such strip of surface, that is, for the integral
7.1.1 The Rectangle Rule: The crudest and most obvious method of approximating the integral I is directly linked to the definition of the integral; we replace the area of the strip In by the rectangle of area fnh and obtain for the integral the approximate expression
From here on, the symbol » means is approximately equal to.
7.1.2 The Trapezoid and Tangent Formulae: We obtain a closer approximation with no greater trouble if we do not replace the area of the strip In , as above, by a rectangular area, but by the trapezoid of area (fn+fn+1)h/2 (Fig. 1.) For the whole integral, this process yields the approximate expression
(trapezoid formula), since, when the areas of the trapezoids are added, each value of the function except the first and the last occurs twice.
As a rule, the approximation becomes even better if, instead of choosing the trapezoid under the chord AB as an approximation to the area of I, we select the trapezoid under the tangent to the curve at the point with the abscissa x=xn+h/2. The area of this trapezoid is simply hfn+½, so that the approximation for the entire integral is
which is called the tangent formula.
7.1.3
Simpson's Rule: By means of
Simpson's rule, we arrive with very little more trouble at a
numerical result which is generally much more exact. This rule
depends on estimating the area In + In+1
of the double strip between the abscissae x = xn and x
= xn+ 2 = xn+2 by
considering the upper boundary to be no longer a straight line
but a parabola - in order to be specific, that parabola which
passes through the three points of the curve with the abscissae 
(Fig. 2).
The equation of this parabola is
The student should verify by direct substitution that for the three values of x in question this eqnation gives the proper values of y, i.e., fn , fn+1 and fn+2, respectively.
If we integrate this second degree polynomial between the xn and xn + 2h, we obtain, after a brief calculation, for the area under the parabola
This represents the required approximation to the area of our strip In + In+1.
If we now assume that n = 2m, i.e., that n is an even number, we obtain, by addition of the areas of such strips, Simpson's rule
7.1.4
Examples: We now apply these
methods to the calculation of
loge2 = 
. If we divide this integral from 1 to 2
into ten equal parts, h will be equal 1/10 and we obtain
by the trapezoidal rule
The value, as was to be expected, is too large, since the curve has its convex side turned towards the x-axis.
By the tangent rule, we find
Owing to the convexity of the curve, this value is too small.
For the same set of sub-divisions, we obtain the most accurate result by means of Simpson's rule:
As a matter of fact,
loge2 =.693147···.
7.1.5 Estimation of the Error:. It is easy to obtain an estimate of the error for each of our methods of integration, if the derivatives of the function f(x) are known throughout the interval of the integration. We take M1, M2, ··· as the upper bounds of the absolute values of the first, second, ··· derivative, respectively, i.e., we assume that throughout the interval |f (n)(x)| < Mn . Then the estimates are:
For the rectangle rule:
For the tangent rule:
For the trapezoid rule:
For Simpson's rule
The last two estimates also lead to estimates for the entire integral 1. We see that Simpson's rule has an error of much higher order in the small quantity h than the other rules, so that when M4 is not too large, it is very advantageous for practical calculations. In order to avoid tiring the reader with the details of the proofs of these estimates, which are fundamentally quite simple, we shall restrict ourselves to the proof for the tangent formula. For this purpose, we expand the function f(x) in the (n + l)-th strip by Taylor's theorem:
where x is a certain intermediate value in the strip. If we integrate the right hand side over the interval xn£ x £ xn + h, the integral of the middle term is zero. Since, as is easily verified,
it follows immediately that
which proves the above assertion.
1. From the formula
compute p with h = 0.1 by (a) the trapezoid rule, (b) Simpson's rule.
2. Compute 
numerically to within 1/100 (10.6.5).
3. Compute 
—- numerically with an error of less
than 0.1.
7.2 Applications of the Mean Value theorem and of Taylor's Theorem
7.2.1 The Calculus of Errors: We now come to quite a different type of numerical calculations. These are applications of the mean value theorem or, more generally, of Taylor's theorem with remainder or, finally, of the infinite Taylor series. As an application which, although simple, is quite important in practice, we shall consider the calculus of errors. This rests on the idea - which lies at the root of the entire differential calculus - that a function f(x), which is differentiable a sufficient number of times, can be represented in the neighbourhood of a point by a linear function with an error of order less than the first, by a quadratic function with an error of order less than the second, etc. Consider the linear approximation to a function y=f(x). If y+Dy=f(x+Dx)=f(x+h), we find by Taylor's theorem
where x = x+ q h (0 < q < 1) is an intermediate value, which need not be known more precisely. If h = Dx is small, we obtain as a practical approximation
In other words, we replace the difference quotient by the derivative to which it is approximately equal and the increment of y by the approximately equal linear expression in h.
We use this fairly obvious fact for practical purposes in the following way. Let two physical quantities x and y be linked by the relation y = f(x). There arises then the question how an inaccuracy in the measurement of x affects the determination of y. If we happen to use instead of the true value of x the inaccurate value x + h, then the corresponding value of y differs from the true value y = f(x) by Dy = f(x + h) - f(x), whence the error is given approximately by the above relation.
We shall understand the use of these relations better if we consider a few examples.
Example 1. The tangent galvanometer: Whren we determine an electric current by a tangent galvanometer, we use the formula y = c tan a, where a is the angle of deflection of the magnetic needle, c is the constant of the apparatus and y = I the intensity of the current. Then
whence Dy » cDa/a². The percentage error in the measurement is given by
We see from this that the accuracy reaches the largest possible value, i.e., that there corresponds to a given error in the measurement of the angle the least possible error in the determination of the current, when the angle a is equal to p/4 or 45º.
In particular, let it be possible
to read the tangent galvanometer to within half a degree; then |Da| is
radians < ½´.01745··· and the percentage error will be 1.745/sin
2a. If the galvanometer reads 30°,
and the percentage
error is less than 2´1.745/1.732, i.e., about 2%.
Example 2. In a triangle ABC (Fig. 3), let the sides b
and c have been measured accurately, while the angle a = x
can only be measured with an error of |Dx|<d. Within
what error range, does the value 
We have
hence the percentage error is
If we consider the special case where b = 400 m., c = 500. and a = 60°, then, by the cosine formula, y = a = 458.2576 m and
If Dx can be measured to within ten seconds of arc, i.e., if
we find that at the worst
i.e., tbe error is at most about 0.004 %.
Example 3. The following illustrates a type of application of the above methods by which we can often save ourselves considerable trouble in physical problems.
It is known experimentally that, if an iron rod has the length l0 at temperature 0, then at temperature t its length will be l = l0(1 + at), where a depends only on the material of the rod. Now, if a pendulum clock keeps correct time at temperature t1, how many seconds will be lost per day if the temperature rises to t2?
We have for the period of oscillation of a pendulum the formula
Hence, if the change of length is Dl, the corresponding change in the period of oscillation is
where l1 = l0(1 + at1) and Dl =a l0(t2 - t1). This is the time lost during each oscillation. The time lost per a second is DT/T » Dl/2l1, whence in one day the clock loses 43,200 Dl/l1 seconds.
The application of our methods has saved here a number of multiplications and two extractions of square roots. Moreover, in the longer direct process, we should finally have to subtract T(l1) from the almost equal value T(l2) and a very small error in calculation would cause a relatively large percentage error in the result.
It is a point of this nature which makes the calculations of applied optics so extremely laborious.
In this case as well as in most cases where the function under consideration has several factors or fractional indices, we can reduce the calculation even more by taking before differentiation logarithms on both sides. In the present example, we have
differentiating, we find
Replacement of dT/dl by DT/Dl yields
in agreement with the preceding result.
7.2.2 Evaluation of p : Gregory's series, also called Leibnitz' series,
which we have obtained in 6.1.2 using the series for artan, is not suitable for the calculation of p due to the slowness of its convergence. However, we may calculate p with comparative ease by the artifice: From the addition theorem for the tangent
changing to the inverse functions a = artan u, b = artan v, we obtain
If we now choose u and v in such a manner that (u + v)/(1 - uv) = 1, we obtain the value p/4 on the right hand side and, if u and v are small numbers, we can easily compute the left hand side by means of known series.
For example, if we set u = 1/2, v = 1/3, as Euler did, we obtain
Moreover, if we note that
we have
so that
Using this formula, Georg Vega (1756 - 1802) calculated the number p to 140 places.
Moreover, using the equation
we obtain
or
This expansion is extremely useful for the computation of p by means of the series
in fact, if we substitute for x the values 1/5, 1/7 or 1/8, we obtain with just a few terms a high degree of accuracy, since the terms diminish rapidly. However, we can perform the calculation even more conveniently by basing it on the formula
obtained by similar arguments as those above.
7.2.3 Calculation of Logarithms: For the computation of logarithms, we transform the series
where 0 < x < 1, by the substitution
into the series
where 2p² - 1 > 1, that is, p² > 1. If p is an integer and p + 1 can be resolved into smaller integral factors, this last series expresses the logarithm of p in terms of the logarithms of smaller integers plus a series, the terms of which diminish very rapidly and the sum of which can therefore be calculated accurately enough from only a few terms. Hence this series allows us to compute the logarithms of any prime number and hence that of any number, provided we have already calculated the value of log 2.
The accuracy of this determination of log p
can be estimated more easily by means of the geometric series
than from the general remainder formula. In fact, we have for the
remainder Rn of the series, i.e.,
the sum of all the terms following the term 
and this formula yields immediately the required error estimate.
For example, let us compute loge 7, using the first four terms of the series. We have
whence
loge7 » 1.945,910,15.
Estimation of the error yields
However, we must note that each of the four numbers which we have added is only given to within an error of 5 ´ 10-9, so that the last place in the above value of loge 7 might be wrong by 2. However, as a matter of fact, the last place is also correct.
1. In order to measure the height of a hill, a tower 100 m. high on top of it is observed from the plain. The angle of elevation of the base of the tower is 42º and the tower itself subtends an angle of 6°. What are the limits of error in the determination of the height if the angle 42° is subject to an error of 1º?
2. Compute loge 2 to three decimal places by means of an expansion in series.
3. Compute loge 5 to six decimal places, using the values of loge 2 and loge 3 given in the text.
4. Compute p to five decimal places, using any one of the formulae in 7.2.2.
7.3 Numerical Solution of Equations
In conclusion, we shall add some remarks about the numerical solution of the equation f(x) = 0, where f(x) need not necessarily be a polynomial and we are, of course, only concerned with the determination of real roots. Every such numerical method is based on the scheme of starting with some known approximation x0 of one of the roots and then improving this approximation, whence this first approximation for the desired root of the equation is found; the quality of the approximation is not a concern. We may perhaps take a rough guess as a first approximation or, better still, obtain it from the graph of the function y =f(x), the intersection of which with the x-axis yields the required root (of course with an error depending on the scale and the accuracy of the drawing).
7.3.1 Newton's Method: The following procedure, due to Newton, is based on the fundamental principle of the differential calculus - the replacement of a curve by a straight line, the tangent - in the immediate neighbourhood of the point of contact. If we have an approximate value x0 for a root of the equation f(x) = 0, we consider the point on the graph of the function y = f(x) with the co-ordinates x = x0, y = f (x0). We wish to find the intersection of the curve with the x-axis; we find an approximation at the point x0, y0=f(x0), where the tangent intersects the x-axis. The abscissa x1 of this intersection of the tangent with the x-axis then represents a new, possibly better approximation to the required root of the equation than x0 .
By virtue of the
geometrical meaning of the derivative, Fig. 4 yields at once
whence we obtain the formula for the calculation of the new value x1:
If this procedure yields a better approximation than x0, we repeat the process and find x2, etc.; if the curve has the form shown in Fig. 4, these approximations will approach more and more closely to the required solution.
The usefulness of
this process depends essentially on the nature of the curve y=f(x).
In Fig. 4, we see that the successive estimates converge with
greater and greater accuracy to the required root. This is due to
the fact that the curve has its convex side turned towards the x-axis.
However, in Fig. 6, we see that if we choose the original value x0
badly, our construction does not lead at all to the required
root. We see from this that while using Newton's method we must
examine each individual case in order to determine the degree of
accuracy with which we have really solved the equation. We shall
return to this topic later
on.
7.3.2 The Rule of False
Position: Newton's method, in
which the tangent to the curve has a decisive role, is only the
limiting case of an older method,
known
as the rule of false
position, in which the secant takes the
place of the tangent. Let us assume that we know two points (x0,y0)
and (x1,y1) in the
neighbourhood of the required intersection with the x-axis.
If we replace the curve by the secant joining these two points,
the intersection of this secant with the x-axis will
under certain circumstances be an improved approximation to the
required root of the equation. If the abscissa of this point is
denoted by x, we have (Fig. 6) the equation
from which we compute x :
or
This formula, which determines the next approximation x from x0 and x1, is called the rule of false position. We can employ it with advantage if one value of the function is positive and the other negative, say, as in Fig. 6, where y0 > 0 and y1 < 0. Repetition of this process will always lead us to the required result, if we use at each step a positive and a negative value of the function, between which the required root must necessarily lie.
The above formula of Newton results from the
rule of false position as a limiting case, if we let x1
tend to x0. In fact, the denominator of the
second term on the right hand side of the statement of the rule
of false position tends to
f '(x0) as x1 tends
to x0.
7.3.3 The Method of Iteration: Another way of approximating roots of an equation f(x) = 0 is the method of iteration. i.e., set f (x) = f (x) + x and write it in the form x = f (x). We then assume that x is the true value of a solution of our equation and x0 a first approximation. We obtain a second approximation x1 by setting x1 = f (x0), a third approximation x2 by setting x2 = f (x1), etc. In order to investigate the convergence of these approximations, we apply the mean value theorem; recalling that x = f (x), we have
where 
lies between x and x0.
This shows that, if for
the absolute value of the derivative f '(x) is less than k < 1, then the successive approximations converge, because
whence the errors tend to zero, The smaller is the absolute value of the derivative f '(x) in the neighbourhood of x, the more rapid is the convergence.
If f '(x) > 1 in the neighbourhood of x , the approximations do not any longer tend to x. We can then use the inverse function or else the following device: We choose a first approximation x0, calculate A =f '(x0) and write
Then the equation f (x)=0 can be written in the form x=f (x), and now f'(x)=-f'(x)/A +1, which has the value 0 at x = x0, whence its absolute value will usually be less than a constant k < 1, provided |x - x| < |x - x0|.
Returning to Newton's method, we can now investigate its suitability for application at any given point. The equation f(x) = 0 is equivalent to
provided that f '(x) ¹ 0. Applying the method of iteration to this last equation, we obtain from a first approximation x0 a second approximation
x1 = x0 - f (x0)/f '(x0);
in other words, the same second approximation as Newton's method gives when it is applied to the equation f(x) = 0. We thus see that the smaller is the value of
the more rapidly converge the successive approximations. In words, Newton's formula converges rapidly for large values of f '(x0) and small values of f (x0) and curvature, as intuition would lead us to suspect.
We can also obtain an estimate of the accuracy of Newton's method, if we recall that, since f(x) = 0, the derivative f '(x) = 0. Applying Taylor's theorem, we have
where lies between x and x0.
Thus, if the error of the original estimate is small, the method
converges much more rapidly than the method of iteration applied
directly to f(x) = 0.
is everywhere less than 10, then a first approximation, which has an error by less than 0.001, will yield a second approximation with an error of less than (0·001)²´10¸2=0.000,005.
Consider as an example the equation
For x0 = 2, we have f(x0) = -1, while, for x1 = 2.l, we have f(x1) = 0.061. By Newton's method,
In order to estimate the error, we find from the expression (a) that the value of f"(x) is about 1 and certainly less than 2 near x = 2. Moreover, the error of our first approximation is certainly less than 1/160, because the secant joining the points x = 2, y = -1 and x = 2.1, y = -0.61 cuts the x-axis at a distance of less than 1/160 from x = 2.1 and the curve, lying under the secant, cuts it even closer to 2.1. Thus, the error of our second approximation is less than
Another way of estimating the error, without reference to the secant, is as follows: If we estimate the error to be less than 1/20, the error of our second approximation is less than 1/20² = 0.0025. Hence the root differs from 2.1 by less than (2.1 - 2.0945) + -0.0025 = 0.008. Hence the error was not merely less than 1/20, but less than 0.008, so that the error in x2 is less than (0.008)²=0.000,064.
If accuracy is insufficient, we can repeat the process, calculating f (x2) and f '(x2) for x2 = 2.094,569 and obtain a third approximation x3 with an error less than 1/(25,600)² < 0.000,000,002.
As a second example, let us solve the equation
We have f(3) = -0.6 and f(4) = +0.4, whence we use x0 = 3.5 as a first approximation. Then, using ten-figure logarithmic tables, we obtain the successive approximations
1. Using Newton's method, find the positive root of x6 + 6x - 8 = 0 to four decimal places.
2. Find to four decimal places the root of x = tan x between p and 2p. Prove that the result is accurate to four places.
3. Using Newton's method, find the value of x for which
4. Find the roots of the equation x = 2 sin x to two decimal places.
5. Determine the positive roots of the equation x5 - x - 0-2 = 0 by the method of iteration.
6. Determine the least positive root of x4 - 3x3 - 10x - 10 = 0 by the method of iteration.
7. Find the roots of x3 - 7x2 + 6x + 20 = 0 to four decimal places.