This namespace builds on the ideas introduced in riemann.cljc
.
riemann.cljc
described four different integration schemes ({left, right, upper, lower} Riemann sums) that were each conceptually simple, but aren't often used in practice, even in their "accelerated" forms.
One reason for this is that their error terms fall off as Loading..., where Loading... is the width of an integration slice. Each order of sequence acceleration can cancel out one of these terms at a time; but still, the performance is not great.
It turns out that by taking the midpoint of each interval, instead of either side, you can reduce the order of the error series to Loading.... This is too good to pass up.
Additionally, because the error terms fall off as Loading..., each order of acceleration is worth quite a bit more than in the Riemann sum case.
This namespace follows the same development as riemann.cljc
:
Here's an implementation of a function that can take the midpoint of a single slice:
And a full (though inefficient) integrator using windowed-sum
:
Let's integrate a triangle!
=> true
It turns out that we already had to implement an efficient version of midpoint-sum
in riemann.cljc
; the incremental version of left and right Riemann sums added the midpoints of each interval when doubling the number of slices.
We can check our implementation against qr/midpoint-sum
:
We'll use qr/midpoint-sum
in the upcoming functions.
Unlike the left and right Riemann sums, the Midpoint method can't reuse function evaluations when the number of slices doubles. This is because each evaluation point, on a doubling, becomes the new border between slices:
If you triple the number of slices from Loading... to Loading..., you can in fact reuse the previous Loading... evaluations:
By scaling Loading... down by a factor of 3, and adding it to a new sum that only includes the new points (using the new slice width).
BTW: The only place I found this idea mentioned is in Section 4.4 of Press's "Numerical Recipes". I haven't found other references to this trick, or implementations. I'd love to hear about them (via a Github issue) if you find any!
We'll follow the interface we used for qr/Sn->S2n
and write Sn->S3n
. This function of Loading... will return a function that performs the incremental update.
The returned function generates Loading... across Loading... with Loading... intervals, and picking out two new points at Loading... and Loading... of the way across the old interval. These are the midpoints of the two new slices with width Loading....
Sum them all up and add them to Loading... to generate Loading...:
Now we can write midpoint-sequence
, analogous to qr/left-sequence
. This implementation reuses all the tricks from qr/incrementalize
; this means it will be smart about using the new incremental logic any time it sees any Loading... multiple of 3, just as the docstring describes.
The following example shows that for the sequence Loading... (used in the Bulirsch-Stoer method!), the incrementally-augmented midpoint-sequence
only performs 253 function evaluations, vs the 315 of the non-incremental (midpoint-sum f2 0 1)
mapped across the points.
The final version is analogous the qr/left-integral
and friends, including an option to :accelerate?
the final sequence with Richardson extrapolation.
I'm not sure what to call this accelerated method. Accelerating the trapezoid method in this way is called "Romberg integration". Using an Loading... sequence of powers of 2 and accelerating the midpoint method by a single step - taking the second column (index 1) of the Richardson tableau - produces "Milne's method".
The ability to combine these methods makes it easy to produce powerful methods without known names. Beware, and enjoy!
We noted above that the the terms of the error series for the midpoint method increase as Loading...... Because of this, we pass Loading... into pr/richardson-sequence
below. Additionally, integral
hardcodes the factor of 3
and doesn't currently allow for a custom sequence of Loading.... This requires 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 midpoint evaluations, investigate polynomial.cljc
or rational.cljc
. The "Bulirsch-Stoer" method uses either of these to extrapolate the midpoint method using a non-geometric sequence.
If you start with the midpoint method, one single step of Richardson extrapolation (taking the second column of the Richardson tableau) is equivalent to "Milne's rule" (see milne.cljc
).
The full Richardson-accelerated Midpoint method is an open-interval variant of "Romberg integration" (see romberg.cljc
).
See the wikipedia entry on Open Newton-Cotes Formulas for more details.