A Newton-Raphson root finder computes the eccentric anomaly from Kepler's equation…

PHYSICS · Multi-part · Fermi University Physics · Mathematical Methods for Physics · Computational Physics · Research software practice, verification and validation

Scientific Computing for Physicists, Ch Is the equation solved right? Testing, debugging and code verification

Problem

A Newton-Raphson root finder computes the eccentric anomaly

E

from Kepler's equation

M = E - e\sin E

for a highly eccentric orbit with

e = 0.9

, at the test case

M = 0.01\,\mathrm{rad}

, whose root is

E = 0.0986\,\mathrm{rad}

. At the converged root, a rounding error in evaluating the residual

E - e\sin E - M

, of order

\epsilon_{\rm mach}|E|

, is amplified into

E

by the local condition number

\kappa = 1/(1 - e\cos E)

. Chapter 2's backward-error bound for an iterative root finder that has converged to machine precision gives a relative error bound of

c\,\kappa\,\epsilon_{\rm mach}

, with

c = 6

absorbing the arithmetic operations of the final Newton step, and

\epsilon_{\rm mach} = 2^{-52}

. ((a)) Evaluate

\kappa

at the test case, derive the tolerance

c\,\kappa\,\epsilon_{\rm mach}/2

for a unit test comparing this solver's output with an independent high-precision solve, and give its numerical value. ((b)) A previous bug in this solver (a wrong sign in the derivative

1 - e\cos E

used for the Newton update) was fixed in a commit that added exactly the

e = 0.9

test case described here. Using the repository workflow, state why this test case must be added to the continuous-integration job that runs on every pull request, rather than simply run once by hand after the fix, for it to serve as a genuine regression test.

Browse the problem library