The Trapezoid Method

ToC

TOC

Simple Implementation

Efficient Trapezoid Method

Incremental Trapezoid Rule

Final Trapezoid API:

Note on Richardson Extrapolation

Next Steps

Back to index/Generated with Clerk from src/emmy/numerical/quadrature/trapezoid.cljc@f807468

The Trapezoid Method

This namespace builds on the ideas introduced in riemann.cljc and midpoint.cljc, and follows the pattern of those namespaces:

implement a simple, easy-to-understand version of the Trapezoid method
make the computation more efficient
write an incremental version that can reuse prior results
wrap everything up behind a nice, exposed API

Let's begin.

Simple Implementation

A nice integration scheme related to the Midpoint method is the "Trapezoid" method. The idea here is to estimate the area of each slice by fitting a trapezoid between the function values at the left and right sides of the slice.

Alternatively, you can think of drawing a line between Loading... and Loading... and taking the area under the line.

What's the area of a trapezoid? The two slice endpoints are

Loading... and
Loading...

The trapezoid consists of a lower rectangle and a capping triangle. The lower rectangle's area is:

Just like in the left Riemann sum. The upper triangle's area is one half base times height:

The sum of these simplifies to:

Or, in Clojure:

(defn single-trapezoid [f xl xr]

  (g// (g/* (g/- xr xl)

            (g/+ (f xl) (f xr)))

2))

#object[emmy.numerical.quadrature.trapezoid$single_trapezoid 0x5a2ca371 "

emmy.numerical.quadrature.trapezoid$single_trapezoid@5a2ca371"

]

We can use the symbolic algebra facilities in the library to show that this simplification is valid:

(->clerk-only

 (let [f (f/literal-function 'f)

       square    (g/* (f 'x_l)

                      (g/- 'x_r 'x_l))

       triangle  (g/* (g// 1 2)

                      (g/- 'x_r 'x_l)

                      (g/- (f 'x_r) (f 'x_l)))]

   (g/simplify

    (g/- (single-trapezoid f 'x_l 'x_r)

         (g/+ square triangle)))))

We can use qr/windowed-sum to turn this function into an (inefficient) integrator:

(defn- trapezoid-sum* [f a b]

  (qr/windowed-sum (partial single-trapezoid f)

                   a b))

#object[emmy.numerical.quadrature.trapezoid$trapezoid_sum_STAR_ 0x786cd893 "

emmy.numerical.quadrature.trapezoid$trapezoid_sum_STAR_@786cd893"

]

Fitting triangles is easy:

(->clerk-only

 ((trapezoid-sum* identity 0.0 10.0) 10))

In fact, we can even use our estimator to estimate Loading...:

(def ^:private pi-estimator*

  (let [f (fn [x] (/ 4 (+ 1 (* x x))))]

    (trapezoid-sum* f 0.0 1.0)))

#object[emmy.numerical.quadrature.riemann$windowed_sum$fn__81025 0x746637c0 "

emmy.numerical.quadrature.riemann$windowed_sum$fn__81025@746637c0"

]

The accuracy is not bad, for 10 slices:

(->clerk-only

 (pi-estimator* 10))

3.1399259889071587

(->clerk-only

 (- Math/PI (pi-estimator* 10)))

0.0016666646826344333

10000 slices gets us closer:

(->clerk-only

 (- Math/PI (pi-estimator* 10000)))

1.6666663604780751e-9

Fun fact: the trapezoid method is equal to the average of the left and right Riemann sums. You can see that in the equation, but lets verify:

(->clerk

 (defn- basically-identical? [l-seq r-seq]

   (every? #(< % 1e-15)

           (map - l-seq r-seq))))

#'emmy.numerical.quadrature.trapezoid/basically-identical?

(->clerk-only

 (let [points  (take 5 (iterate inc 1))

       average (fn [l r]

                 (/ (+ l r) 2))

       f       (fn [x] (/ 4 (+ 1 (* x x))))

       [a b]   [0 1]

       left-estimates  (qr/left-sequence f a b {:n points})

       right-estimates (qr/right-sequence f a b {:n points})]

   (basically-identical?

    (map (trapezoid-sum* f a b) points)

    (map average

         left-estimates

         right-estimates))))

true

Efficient Trapezoid Method

Next let's attempt a more efficient implementation. Looking at single-trapezoid, it's clear that each slice evaluates both of its endpoints. This means that each point on a border between two slices earns a contribution of Loading... from each slice.

A more efficient implementation would evaluate both endpoints once and then sum (without halving) each interior point.

This interior sum is identical to a left Riemann sum (without the Loading... evaluation), or a right Riemann sum (without Loading...).

Here is this idea implemented in Clojure:

(defn trapezoid-sum

  "Returns a function of `n`, some number of slices of the total integration

  range, that returns an estimate for the definite integral of $f$ over the

  range $(a, b)$ using the trapezoid method."

  [f a b]

  (let [width (- b a)]

    (fn [n]

      (let [h  (/ width n)

            fx (fn [i] (f (+ a (* i h))))]

        (* h (+ (/ (+ (f a) (f b)) 2)

                (ua/sum fx 1 n)))))))

#object[emmy.numerical.quadrature.trapezoid$trapezoid_sum 0x3669ac45 "

emmy.numerical.quadrature.trapezoid$trapezoid_sum@3669ac45"

]

We can define a new pi-estimator and check it against our less efficient version:

(def ^:private pi-estimator

  (let [f (fn [x] (/ 4 (+ 1 (* x x))))]

    (trapezoid-sum* f 0.0 1.0)))

#object[emmy.numerical.quadrature.riemann$windowed_sum$fn__81025 0x961a6d2 "

emmy.numerical.quadrature.riemann$windowed_sum$fn__81025@961a6d2"

]

(->clerk-only

 (basically-identical?

  (map pi-estimator (range 1 100))

  (map pi-estimator* (range 1 100))))

true

Incremental Trapezoid Rule

Next let's develop an incremental updater for the Trapezoid rule that lets us reuse evaluation points as we increase the number of slices.

Because interior points of the Trapezoid method mirror the interior points of the left and right Riemann sums, we can piggyback on the incremental implementations for those two methods in developing an incremental Trapezoid implementation.

Consider the evaluation points of the trapezoid method with 2 slices, next to the points of a 4 slice pass:

x-------x-------x
x---x---x---x---x

The new points are simply the midpoints of the existing slices, just like we had for the left (and right) Riemann sums. This means that we can reuse qr/Sn->S2n in our definition of the incrementally-enabled trapezoid-sequence:

(defn trapezoid-sequence

  "Returns a (lazy) sequence of successively refined estimates of the integral of

  `f` over the open interval $(a, b)$ using the Trapezoid method.

  ### Optional arguments:

  `:n`: If `:n` is a number, returns estimates with $n, 2n, 4n, ...$ slices,

  geometrically increasing by a factor of 2 with each estimate.

  If `:n` is a sequence, the resulting sequence will hold an estimate for each

  integer number of slices in that sequence.

  `:accelerate?`: if supplied (and `n` is a number), attempts to accelerate

  convergence using Richardson extrapolation. If `n` is a sequence this option

  is ignored."

  ([f a b] (trapezoid-sequence f a b {:n 1}))

  ([f a b {:keys [n accelerate?] :or {n 1}}]

   (let [S      (trapezoid-sum f a b)

         next-S (qr/Sn->S2n f a b)

         xs     (qr/incrementalize S next-S 2 n)]

     (if (and accelerate? (number? n))

       (pr/richardson-sequence xs 2 2 2)

       xs))))

#object[emmy.numerical.quadrature.trapezoid$trapezoid_sequence 0x17575f2d "

emmy.numerical.quadrature.trapezoid$trapezoid_sequence@17575f2d"

]

The following example shows that for the sequence Loading..., the incrementally-augmented trapezoid-sequence only performs Loading... function evaluations; i.e., the same number of evaluations as the non-incremental (trapezoid-sum f2 0 1) would perform for Loading... slices. (why Loading...? each interior point is shared, so each trapezoid contributes one evaluation, plus a final evaluation for the right side.)

The example also shows that evaluating every Loading... in the sequence costs Loading... evaluations. As Loading... gets large, this is roughly twice what the incremental implementation costs.

When Loading..., the incremental implementation uses 2049 evaluations, while the non-incremental takes 4017.

(->clerk-only

 (let [n-elements 11

       f (fn [x] (/ 4 (+ 1 (* x x))))

       [counter1 f1] (u/counted f)

       [counter2 f2] (u/counted f)

       [counter3 f3] (u/counted f)

       n-seq (take (inc n-elements)

                   (iterate (fn [x] (* 2 x)) 1))]

   ;; Incremental version evaluating every `n` in the sequence $1, 2, 4, ...$:

   (dorun (trapezoid-sequence f1 0 1 {:n n-seq}))

   ;; Non-incremental version evaluating every `n` in the sequence $1, 2, 4, ...$:

   (run! (trapezoid-sum f2 0 1) n-seq)

   ;; A single evaluation of the final `n`

   ((trapezoid-sum f3 0 1) (last n-seq))

   (let [two**n+1 (inc (g/expt 2 n-elements))

         n+2**n (+ n-elements (g/expt 2 (inc n-elements)))]

     [[two**n+1 n+2**n two**n+1]

      [@counter1 @counter2 @counter3]])))

[[2049 4107 2049] [2049 4107 2049]]

Another short example that hints of work to come. The incremental implementation is useful in cases where the sequence includes doublings nested in among other values.

For the sequence Loading... (used in the Bulirsch-Stoer method!), the incrementally-augmented trapezoid-sequence only performs 162 function evaluations, vs the 327 of the non-incremental (trapezoid-sum f2 0 1) mapped across the points.

This is a good bit more efficient than the Midpoint method's incremental savings, since factors of 2 come up more often than factors of 3.

(->clerk-only

 (let [f (fn [x] (/ 4 (+ 1 (* x x))))

       [counter1 f1] (u/counted f)

       [counter2 f2] (u/counted f)

       n-seq (take 12 (interleave

                       (iterate (fn [x] (* 2 x)) 2)

                       (iterate (fn [x] (* 2 x)) 3)))]

   (dorun (trapezoid-sequence f1 0 1 {:n n-seq}))

   (run! (trapezoid-sum f2 0 1) n-seq)

   [@counter1 @counter2]))

[162 327]

Final Trapezoid API:

The final version is analogous the qr/left-integral and friends, including an option to :accelerate? the final sequence with Richardson extrapolation. (Accelerating the trapezoid method in this way is called "Romberg integration".)

Note on Richardson Extrapolation

The terms of the error series for the Trapezoid method increase as Loading...... (see https://en.wikipedia.org/wiki/Trapezoidal_rule#Error_analysis). Because of this, we pass Loading... into pr/richardson-sequence below. Additionally, integral hardcodes the factor of 2 and doesn't currently allow for a custom sequence of Loading.... This is configured by passing Loading... into pr/richardson-sequence.

If you want to accelerate some other geometric sequence, call pr/richardson-sequence with some other value of t.

To accelerate an arbitrary sequence of trapezoid evaluations, investigate polynomial.cljc or rational.cljc. The "Bulirsch-Stoer" method uses either of these to extrapolate the Trapezoid method using a non-geometric sequence.

(qc/defintegrator integral

  "Returns an estimate of the integral of `f` over the closed interval $[a, b]$

  using the Trapezoid method with $1, 2, 4 ... 2^n$ windows for each estimate.

  Optionally accepts `opts`, a dict of optional arguments. All of these get

  passed on to [[emmy.util.stream/seq-limit]] to configure convergence

  checking.

  See [[trapezoid-sequence]] for information on the optional args in `opts` that

  customize this function's behavior."

  :area-fn single-trapezoid

  :seq-fn trapezoid-sequence)

#object[emmy.numerical.quadrature.common$make_integrator_fn$call__53914 0x52cab230 "

emmy.numerical.quadrature.common$make_integrator_fn$call__53914@52cab230"

]

Next Steps

If you start with the trapezoid method, one single step of Richardson extrapolation (taking the second column of the Richardson tableau) is equivalent to "Simpson's rule". One step using t=3, i.e., when you triple the number of integration slices per step, gets you "Simpson's 3/8 Rule". Two steps of Richardson extrapolation gives you "Boole's rule".

The full Richardson-accelerated Trapezoid method is also known as "Romberg integration" (see romberg.cljc).

These methods will appear in their respective namespaces in the quadrature package.

See the wikipedia entry on Closed Newton-Cotes Formulas for more details.