NORDISK MATEMATISK TIDSKRIFT
Volume 20, pages 138-142 (Oslo 1972)

GENERALIZED SIMPSON'S RULE

SUMMARY
       By rotation of the coordinate system in the customary Simpson integration scheme, a generalization of Simpson's rule is deduced to estimate the integral _C  y dx when the curve C is approximated by an ordered set C* of points (x_i, y_i) not necessarily with equidistant abscissas or unique ordinates. If T₁ and T₂ are trapezoidal estimates of _C  y dx using respectively all points and every second point of C*, the generalized Simpson's rule is shown to be _C  y dx  T₁ + (T₁-T₂) /3.

1. Introduction.
       Let a curve C in the xy-plane be approximated by a set C* of points (x_i , y_i) along the curve. This paper studies the problem of estimating the integral  _C  y dx by using the points of C*.
       We consider first the case when C is described by a unique function of x and the points of C* are given at equidistant intervals in x. The wellknown Simpson's rule, in which the function is approximated by a parabola with axis in the y-direction, gives for three successive points (x_i-1, y_i-1), (x_i, y_i) and (x_i+1, y_i+1) of C*:

(1)              _i-1ⁱ⁺¹  y dx    h (y_i-1 + 4y_i + y_i+1) /3,

where h = x_i+1 - x_i = x_i - x_i-1.
       For the case when the ponts of C* are not equidistant in x, various modifications of Simpson's rule have been suggested, such as Brun's rule [V Brun: A generalization of the formula of Simpson for non- equidistant ordinates. Nordisk Mat. Tidskr. 1 (1953), pp. 10-15] derived by an Archimedean quadrature:

(2)              _i-1ⁱ⁺¹  y dx    (a + b) (y_i-1 + 4y_i + y_i+1) /6 + (a - b) (y_i-1 - y_i+1) /3,

where a = x_i - x_i-1, b= x_i+1 - x_i.
       An analytical derivation of Brun's rule as well as generalizations to higher order formulae has been made by Selmer [E S Selmer: Numerical Integration by non-equidistant ordinates. Nordisk Mat. Tidskr. 6 (1958), pp. 97-108].

 2. A rotation-covariant formulation.
       The new rule presented below, which is applicable also in the general case when C is described by a non-unique function of x, is built on a reduction of the general quadrature problem to the well studied equidistant, unique case. For any three successive points A, B, C of C* in the general case, it is always possible to rotate the coordinate system in such a way that a new orthogonal frame   is formed in which these three points have equidistant  and unique , see fig 1 (the -axis is to be constructed parallel to the median from B in the triangle ABC). In the   -frame we can then apply Simpson's rule (1) and get

(3)              _i-1ⁱ⁺¹   d    h (_i-1 + 4 _i + _i+1) /3,

which we will now give a formulation which is independent of the coordinate system used.
       Let T₁ be the trapezoidal estimate of the left hand side of (3) using the points A, B, C, and similarly T₂ be the trapezoidal estimate using the points A, C, i.e.

                   T₁ = h (_i-1 + 2 _i + _i+1) /3,

                   T₂ = h (_i-1 + _i+1).

From this and (3) then follows that the Simpson estimate can be expressed in the trapezoidal estimates, and we obtain

(4)               _i-1ⁱ⁺¹   d    T₁ + (T₁ - T₂) /3.

       This equation is covariant under rotation of the coordinate system and expresses the geometrical fact (cf Archimedes: Die Quadratur der Parabel. Ostwald's Klassiker der exakten Wissenschaften, Nr. 203, Leipzig 1923, § 17) that if A, B, and C are points on a parabola with its axis parallel to the median from B in the triangle ABC, then the segments S₁ and S₂ in fig. 1 are of equal size and, taken together, equal to 1/3 of the area of the triangle ABC. (This is easily proved by considering an isosceles triangle inscribed in a parabola, fig. 2, for which this is true. Fig. 1 is obtained from fig. 2 by an affine transformation, which will leave the parabolic form invariant and not change the relative proportions between areas.)
       For three points in succession the rule (4) gives, in the equidistant case, Simpson's rule (1), and in the non-equidistant case, Brun's rule. (Selmer has pointed out [op. cit., p. 102] that Brun's rule can be obtained by replacing the given curve by a certain parabola with sloping axis. It is not difficult to see that his parabola is the same as the one in fig. 1.)
       The formulation of (4), however, is valid for an arbitrary number of points of C*, (when calculating T₂ using every second point of C*, possible single excess points in the beginning and/or end of C* are to be included in T₂ as well), and valid ragardless of the shape of C which may have an arbitrary number of loops. We thus obtain the general quadrature rule

(5)              _C  y dx    T₁ + (T₁ - T₂) /3.

       Noting that there are two possible trapezoidal estimates T₂' and T₂" when every second point of C* is used, it is sometimes C* advantageous to use for T₂ in (5) their arithmetic mean, T₂ = (T₂' + T₂")/2.
       An estimate of the error committed in the approximation (5) may be derived, cf. Selmer op. cit. However, such error terms tend to become rather complicated.

 3. Two examples.
       The main value of a non-equidistant integration formula lies in the fact that any set of predetermined points on the integration contour can be used directly without interpolation in the integration scheme, so that for instance a predetermined interval division which is expected to follow the functional behaviour of the integrand can be used. To illustrate this we study the numerical integration of  I = ₀ ¹x^-½ dx. In order to follow the asymptotic behaviour of the integrand at x=0, we use an Achilles-and-the-tortoise type of subdivision in which the intervals are successively halved towards x=0. The x-values are thus 1, 1/2, 1/4, 1/8, 1/16, etc, with the corresponding values of the integrand 1, 2, 2, 22, 4, etc.
       The trapezoidal rule using every point then gives

        T₁ = ½ {(1 + 2 )/2 + (2 + 2)/4 + (2 + 22 )/8 + (22 + 4)/16 + ...}

             = ½ { ½ (1 + 2 )/(1 - ½ 2 )} = 1 + 3 2 / 4  2.061,

whereas using every second point we obtain the trapezoidal estimate

        T₂ = ½ {3 (1 + 2) / 4 + 3 (2 + 4) / 16 + ...}

             = ½ { ³/₄ (1 + 2 ) / (1 - ¹/₄ 2)} =  ⁹/₄  = 2.250.

The two trapezoidal estimates T₁ and T₂ thus deviate by respectively 3% and 12.5% from the exact value I = 2. Using the generalized Simpson's rule (5) presented above, we get

          I    T₁ + (T₁ - T₂) / 3 = ⁹/₄ + 2  1.998,

i.e. 0.1% from the exact value.

       To illustrate the case when the integration contour C is not a unique function, we also consider the numerical quadrature of a unit circle by using eight equally spaced points on the circumference of the circle. The area of the rectangular octagon which is obtained when the eight points are connected corresponds to the trapezoidal estimate T₁ above, and the inscribed square obtained by connecting every second of the eight points correponds to the trapezoidal estimate T₂. Here T₁ = 2 2   2.83 and T₂ = 2.00, values which deviate by about 10% and 36% from the exact value . Using the generalized Simpson's rule we get T₁ + (T₁ - T₂)/3  3.11, i.e. a 1% deviation from the exact value.

Back to Scientor's home page

Copyright © 1972-97 Scientor Innovation AB, Stockholm, Sweden.