I'm writing a library that involves some approximations of variational calculus problems. Whenever I implement routines to evaluate the derivative or Hessian of an action functional $A$, I write a test to check that these are working correctly. To do this, I evaluate the error in using the first and second derivatives of the action to compute a local quadratic approximation to the objective functional:
$E = A(u + \delta v) - \left(A(u) + \delta\cdot dA(u)\cdot v + \frac{1}{2}\delta^2(d^2A(u)v)\cdot v\right)$
and check that $E/\delta^3$ doesn't grow for $\delta = 2^{-1},\ldots,2^{-N}$ for some big enough $N$.
This works well for moderately sized $N$ (up to 16 or so), but for $N$ too large, truncation error begins to dominate and the error in the approximation gets worse.
What is a good way to test derivative approximations in the face of truncation error?
The way I see it, I have only bad options:
- Write an initial version of the test that prints the approximation errors, find the $N$ for which truncation error takes over, then use this $N$ as the cutoff point in the final version of the unit test.
- Have some way of automatically diagnosing where truncation error takes over and discard all the results thereafter.
- Don't make my unit tests return failure unless it's really obvious. Graph the results and let the user decide.
- Determine analytically where truncation error is likely to occur.
I don't like option 1 because it feels like cherry-picking my tests so that they will pass. I don't like option 2 because a genuinely failing implementation's inherent mathematical errors might be erroneously assessed as truncation error, so that automated regression testing would miss the bug. At that rate, I'd have to check the results manually, in which case I might as well go with option 3. But I don't like option 3 because I want automated testing. I think option 4 is practical for functions of a single variable -- evaluating $\sqrt{\epsilon|f(x)/f''(x)|}$ where $\epsilon$ is the machine epsilon gives you a pretty good upper bound -- but for big PDE problems where the unknown is a field with thousands of variables, not so much.