This module builds up to an implementation of numerical derivatives. The final function, D-numeric, uses Richardson extrapolation to speed up convergence of successively tighter estimates of Loading... using a few different methods.
The inspiration for this style was Halfant and Sussman's "Abstraction in Numerical Methods", starting on page 10.
We'll proceed by deriving the methods symbolically, and then implement them numerically.
First, a function that will print nicely rendered infix versions of (simplified) symbolic expressions:
And a function to play with:
The key to all of these methods involves the taylor series expansion of an arbitrary function Loading... around a point Loading...; we know the taylor series will include a term for Loading..., so the goal is to see if we can isolate it.
Here's the taylor series expansions of Loading...:
Use show to print out its Loading... representation:
We can solve this for Loading... by subtracting Loading... and dividing out Loading...:
Voila! The remaining terms include Loading... along with a series of progressively-smaller "error terms" (since Loading...). The first of these terms is Loading.... It will come to dominate the error as Loading..., so we say that the approximation we've just derived has error of Loading....
This particular formula, in the limit as Loading..., is called the "forward difference approximation" to Loading.... Here's the Clojure implementation:
We could also expand Loading...:
and solve for Loading...:
To get a similar method, called the "backward difference" formula. Here's the implementation:
Notice that the two expansions, of Loading... and Loading..., share every term paired with an even power of Loading.... The terms associated with odd powers of Loading... alternate in sign (because of the Loading... in the expansion of Loading...).
We can find yet another method for approximating Loading... if we subtract these two series. We're trying to solve for Loading..., and Loading... appears paired with Loading..., an odd-powered term... so subtracting Loading... should double that term, not erase it. Let's see:
Amazing! Now solve for Loading...:
We're left with Loading..., a quadratic error term in Loading.... (Of course if we'd expanded to more than initial terms in the taylor series we'd see a long error series with only even powers.)
This formula is called the "central difference" approximation to the first derivative. Here's the implementation:
There's one more approximation we can extract from these two expansions. We noted earlier that the terms associated with odd powers of Loading... alternate in sign. If we add the two series, these odd terms should all cancel out. Let's see:
Interesting. The Loading... term is gone. Remember that we have Loading... available; the first unknown term in the series is now Loading.... Solve for that term:
This is the "central difference" approximation to the second derivative of Loading.... Note that the error term here is quadratic in Loading.... Here it is in code:
Let's attempt to use these estimates and see how accurate they are. (This section follows Sussman starting on page 10.)
The following function returns a new function that approximates Loading... using the central difference method, with a fixed value of Loading...:
The error here is not great, even for a simple function:
Let's experiment instead with letting Loading.... This next function takes a function Loading..., a value of Loading... and an initial Loading..., and generates a stream of central difference approximations to Loading... using successively halved values of Loading..., i.e., Loading...
Let's print 20 of the first 60 terms (taking every 3 so we see any pattern):
At first, the series converges toward the proper value. But as Loading... gets smaller, Loading... and Loading... get so close together that their difference is less than the minimum epsilon allowed by the system's floating point representation.
As Sussman states:
"Hence we are in a race between truncation error, which starts out large and gets smaller, and roundoff error, which starts small and gets larger." ~Sussman, p12
We can actually analyze and quantify how many halvings we can apply to Loading... before roundoff error swamps our calculation.
Why does roundoff error occur? From Sussman: "Any real number Loading..., represented in the machine, is rounded to a value Loading..., where Loading... is effectively a random variable whose absolute value is on the order of the machine epsilon Loading...: that smallest positive number for which 1.0 and Loading... can be distinguished."
In the current library, u/machine-epsilon holds this value.
Our goal, then, is to see if we can figure out when the error due to roundoff grows so large that it exceeds the tolerance we want to apply to our calculation.
For the central difference formula:
Loading...without any roundoff error, the numerator should be equal to Loading.... In reality, for small values of Loading..., Loading... and Loading... both have machine representations in error by about Loading.... Their difference doesn't change the order, so we can say that their difference also has error of Loading....
Dividing these two together, the relative error is:
Loading...The relative error doubles each time Loading... is halved. This is technically just the relative error of the numerator of the central difference method, but we know the denominator Loading... to full precision, so we can ignore it here.
If we actually calculate this ratio, we'll find the INITIAL relative error due to roundoff for a given h. Also, because we want to make sure that we're working in integer multiples of machine epsilon, let's actually take the next-highest-integer of the ratio above. The following method takes the ratio above as an argument, and returns:
Loading...That calculation, as the documentation states, returns the number of "roundoff units". Let's call it Loading....
Each iteration doubles the relative error contributed by roundoff. Given some tolerance, how many roundoff error doublings (or, equivalently, halvings of Loading...) can we tolerate before roundoff error swamps our calculation?
Here's the solution:
Let's combine these two ideas into a final function, terms-before-roundoff, that calculates how items we can pull from a sequence like central-diff-stream above before roundoff swamps us. (This is 1 + max iterations, because we need to include the first point.)
How many terms are we allowed to examine for an estimate of the derivative of Loading..., with an initial Loading...?
6 terms, or 5 halvings, down to Loading.... How many terms does the sequence take to converge?
15 is far beyond the level where roundoff error has rendered our results untrustworthy.
We need a way to converge more quickly. richardson.cljc lays out a general method of "sequence acceleration" that we can use here, since we know the arithmetic progression of the terms in the error series for each of our methods above.
For the central difference method, our highest-order error term has an exponent of Loading..., and successive terms are all even. r/richardson-sequence takes p and q for an arithmetic sequence of error exponents Loading...
It also needs the initial size Loading... of our sequence progression.
Given that information, we can transform our existing sequence of estimates into an accelerated sequence of estimates using Richardson extrapolation. Does it converge in fewer terms?
Happily, it does, in only 6 terms instead of 15! This brings convergence in to our limit of 6 total terms.
If you're interested in more details of Richardson extrapolation, please see richardson.cljc! For now we'll proceed.
We're ready to write our final numeric differentiation routine, D-numeric. First, some supporting structure. We support four methods, so let's describe them using keywords in a set:
To apply one of the methods, we need to be able to:
generate the method's estimate as a function of Loading...
calculate the "relative error ratio" that we used above to calculate a maximum number of terms to analyze
know the order Loading... of the highest order error term, and
the increment Loading... of successive error terms
Once again, richardson.cljc for a discussion of Loading... and Loading....
This configs function bundles all of this together. I don't know that this is the best abstraction, but I don't know yet of any other methods for numeric differentiation, so it'll do for now.
Note here that Loading... for both central difference methods, just like we determined above. the forward and backward difference methods both have all of the remaining terms from the taylor expansion in their error series, so they only get Loading....
More resources about numerical differentiation: