Riemann Quadrature

This namespace includes functions for calculating the Riemann integral of a single-variable function. These are probably not methods that you'll want to use; see the documentation and defaults in emmy.numerical.quadrature for good recommendations. But they're clear and understandable. The goal of this namespace is to lay the groundwork for visualizable routines that you can use to step toward understanding of the tougher methods.

"Quadrature", in this context, means "numerical integration". The word is a historical term for calculating the area inside of some geometry shape. Riemann sums are a group of methods for numerical integration that use this strategy:

• partition the area under the curve of some function Loading... into Loading... "slices"
• generate some area estimate for each slice
• add up all of the slices to form an estimate of the integral
• increase the number of slices, and stop when the estimate stops changing.

The Riemann integral of a function Loading... is the limit of this process as Loading....

How do you estimate the area of a slice? All of these methods estimate the area by forming a rectangle. For the base, use Loading.... For the height, you might use:

• the function value at the left point, Loading... (Left Riemann sum)
• the right point, Loading... (Right Riemann sum)
• the max of either Loading... ("upper" Riemann sum)
• the minimum, Loading..., called the "lower" Riemann sums
• the function value at the midpoint: Loading...

This namespace builds up to implementations for left-integral, right-integral, upper-integral and lower-integral. midpoint.cljc holds an implementation of the Midpoint method.

A closely related method involves forming a trapezoid for each slice. This is equivalent to averaging the left and right Riemann sums. The trapezoid method lives in trapezoid.cljc.

Riemann Sum Implementation

We'll start with an inefficient-but-easily-understandable version of these methods. To form a Riemann sum we need to:

• partition some range Loading... into n slices
• call some area-generating function on each slice
• add all of the resulting area estimates together

windowed-sum implements this pattern:

Test this out with a function that returns 2 for every slice, and we get back an estimate (from the function returned by windowed-sum) of 2x the number of slices:

Now, let's implement the four classic "Riemann Integral" methods.

Let's say we want to integrate a function Loading.... The left and right Riemann sums estimate a slice's area as a rectangle with:

• width == Loading..., and
• height == Loading... or Loading..., respectively.

left-sum is simple to implement, given windowed-sum:

Every internal slice has the same width, so we can make the sum slightly more efficient by pulling out the constant and multiplying by it a single time.

Internally, we also generate all of the internal "left" points directly from the slice index, instead of pre-partitioning the range. This is fine since we don't need Loading....

right-sum is almost identical, except that it uses Loading... as the estimate of each rectangle's height:

Same trick here to get a more efficient version. This implementation also generates an internal function fx of the window index. The only difference from the left-sum implementation is an initial offset of h, pushing every point to the right side of the window.

The upper Riemann sum generates a slice estimate by taking the maximum of Loading... and Loading...:

Similarly, the lower Riemann sum uses the minimum of Loading... and Loading...:

Estimating Integrals with Riemann Sums

Given the tools above, let's attempt to estimate the integral of Loading... using the left and right Riemann sum methods. (The actual equation for the integral is Loading...).

The functions above return functions of n, the number of slices. We can use (us/powers 2) to return a sequence of (1, 2, 4, 8, ...) and map the function of n across this sequence to obtain successively better estimates for Loading.... The true value is Loading...:

Both estimates are bad at 32 slices and don't seem to be getting better. Even up to Loading... slices we haven't converged, and are still far from the true estimate:

This bad convergence behavior is why common wisdom states that you should never use left and right Riemann sums for real work.

But maybe we can do better.

Sequence Acceleration

One answer to this problem is to use "sequence acceleration" via Richardson extrapolation, as described in richardson.cljc.

pr/richardson-sequence takes a sequence of estimates of some function and "accelerates" the sequence by combining successive estimates.

The estimates have to be functions of some parameter Loading... that decreases by a factor of Loading... for each new element. In the example above, Loading... doubles each time; this is equivalent to thinking about the window width Loading... halving each time, so Loading....

This library's functional style lets us accelerate a sequence of estimates xs by simply wrapping it in a call to (pr/richardson-sequence xs 2). Amazing!

Does Richardson extrapolation help?

We now converge to the actual, true value of the integral in 4 terms!

This is going to be useful for each of our Riemann sums, so let's make a function that can accelerate a generic sequence of estimates. The following function takes:

• the sequence of estimates, estimate-seq
• a dictionary of "options"

This library is going to adopt an interface that allows the user to configure a potentially very complex integration function by sending a single dictionary of options down to each of its layers. Adopting that style now is going to allow this function to grow to accomodate other methods of sequence acceleration, like polynomial or rational function extrapolation.

For now, {:accelerate? true} configures Richardson extrapolation iff the user hasn't specified a custom sequence of integration slices using the :n option.

Check that this works:

Excellent!

Incremental Computation

The results look quite nice; but notice how much redundant computation we're doing.

Consider the evaluation points of a left Riemann sum with 4 slices, next to a left sum with 8 slices:

Every time we double our number of number of evaluations, half of the windows share a left endpoint. The same is true for a right sum:

In both cases, the new points are simply the /midpoints/ of the existing slices.

This suggests a strategy for incrementally updating a left or right Riemann sum when doubling the number of points:

• Generate a new midpoint estimate of each n slices
• Add this estimate to the previous estimate
• Divide the sum by 2 to scale each NEW slice width down by 2 (since we're doubling the number of slices)

First, implement midpoint-sum. This is very close to the implementation for left-sum; internally the function adds an offset of Loading... to each slice before sampling its function value.

The next function returns a function that can perform the incremental update to a left or right Riemann sum (and to a midpoint method estimate, as we'll see in midpoint.cljc):

After using left-sum to generate an initial estimate, we can use Sn->S2n to generate all successive estimates, as long as we always double our slices. This suggests a function that takes an initial number of slices, n0, and then uses reductions to scan across (us/powers 2 n0) with the function returned by Sn->S2n:

Verify that this function returns an equivalent sequence of estimates to the non-incremental left-sum, when mapped across powers of 2:

Generalizing the Incremental Approach

We need to use the same style for right-sum, so let's try and extract the pattern above, of:

• generating an initial estimate of n0 slices using some function S-fn
• refining an estimate of n0 slices => n0 / 2 slices using some incremental updater, next-S-fn

In fact, because methods like the Midpoint method from midpoint.cljc can only incrementally update from n => n/3, let's make the factor general too.

geometric-estimate-seq captures the pattern above:

And another version of left-sequence, implemented using the new function:

Incremental Updates with Any Sequence

What if we want to combine the ability to reuse old results with the ability to take successively refined estimates that don't look like geometric series? The series 1, 2, 3... of natural numbers is an obvious choice of windows... but only the even powers are able to reuse estimates.

Integration methods like the Bulirsch-Stoer approach depend on sequences like 2, 3, 4, 6...

We absolutely want to be able to save potentially-expensive function evaluations.

One way to do this is to memoize the function f that you pass in to any of the methods above.

Alternatively, we could implement a version of geometric-estimate-seq that takes any sequence of estimate,s and maintains a sort of internal memoization cache.

For every n, check the cache for prev == n/factor. If it exists in the cache, use next-S-fn; else, use S-fn, just like we did in geometric-estimate-seq for the initial value.

general-estimate-seq does this:

We can combine general-estimate-seq and geometric-estimate-seq into a final method that decides which implementation to call, based on the type of the n0 argument.

If it's a number, use it as the n0 seed for a geometrically increasing series of estimates. Else, assume it's a sequence and pass it to general-estimate-seq.

Final Incremental Implementations

We can use incrementalize to write our final version of left-sequence, along with a matching version for right-sequence.

Notice that we're using accelerate from above. The interface should make more sense now:

lower-sequence and upper-sequence are similar. They can't take advantage of any incremental speedup, so we generate a sequence of ns internally and map lower-sum and upper-sum directly across these.

Integral API

Finally, we expose four API methods for each of the {left, right, lower, upper}-Riemann sums.

Each of these makes use a special qc/defintegrator "macro"; This style allows us to adopt one final improvement. If the interval Loading... is below some threshold, the integral API will take a single slice using the supplied :area-fn below and not attempt to converge. See common.cljc for more details.

These API interfaces are necessarily limiting. They force the assumptions that you:

• only want to use geometrical sequences that start with n0 = 1
• only want to (optionally) accelerate using Richardson extrapolation

I can imagine a better API, where it's much easier to configure generic sequence acceleration! This will almost certainly show up in the library at some point. For now, here are some notes:

• Richardson extrapolation requires a geometric series of estimates. If you want to use some other geometry series with left-sequence or right-sequence, you can still accelerate with Richardson. Just pass your new factor as t.

• For each of {left, right, lower, upper}-Riemann sums, the order of the error terms is 1, 2, 3, 4..., so always provide p=1 and q=1 to richardson-sequence. accelerate does this above.

• If you want to use some NON-geometric seq, you'll need to use the methods in polynomial.cljc and rational.cljc, which are more general forms of sequence acceleration that use polynomial or rational function extrapolation. Your sequence of xs for each of those methods should be n-seq.

upper and lower Riemann sums have the same interface; internally, they're not able to take advantage of incremental summation, since it's not possible to know in advance whether or not the left or right side of the interval should get reused.

Next Steps

For a discussion and implementation of the more advanced methods (the workhorse methods that you should actually use!), see midpoint.cljc and trapezoid.cljc. The midpoint method is the standard choice for open intervals, where you can't evaluate the function at its endpoints. The trapezoid method is standard for closed intervals.