Taking N to the limit

A simple numerical algorithm to compute lim_{x → a} f(x) in the case when f(a) is undefined is to evaluate f(a+ε) for some fixed small value ε (or f(1/ε) in case a = ∞). The expression f(a+ε) or f(1/ε) is typically extremely poorly conditioned and hence a challenge for a numerical evaluation routine.

As a stress test for N, I have tried numerically evaluating all the limits in test_demidovich.py, a set of tests for SymPy's symbolic limit function.

I chose to evaluate each limit to 10 accurate digits, using ε = 10^-50. Some simplifications were necessary:

• Since the numerical limit algorithm described above cannot generally detect convergence to 0 or ∞ (giving pseudorandom tiny or huge values instead), I chose to interpret any magnitude outside the range 10^-10 to 10¹⁰ as 0 or ∞.



• SymPy's limit function supports limits containing parameters, such as lim_x→0 (cos(mx)-cos(nx))/x² = (n²-m²)/2. In all such cases, I replaced the parameters with arbitrary values.

The nlimit function with tests is available in the file limtest.py (requiring evalf.py, mpmath and the hg version of SymPy). A straightforward limit evaluation looks like this:

>>> nlimit(sin(x)/(3*x), x, 0, 10)
'0.3333333333'

The results of the tests? After fixing two minor bugs in N that manifested themselves, nlimit passes 50 out of 53 tests. It only fails three tests involving the functions log and asin which are not yet implemented. SymPy's limit function fails 8 out of the 53 tests; in each of these cases, nlimit gives the correct value to 10 digits.

nlimit is 100 times faster than limit, processing all the test cases in 0.18 (versus 18) seconds.

Despite only 10 digits being requested, running N with verbose=True shows that upwards of 800 bits of working precision are required for some of the limits, indicating very clearly the need for adaptive numerical evaluation.

The heuristic of using a fixed, finite ε will not work in case a limit converges extremely slowly. And of course, limit gives a nice, symbolic expression instead of an approximation (nlimit could give an exact answer in simple cases by passing its output through number recognition functions in mpmath). Due to these limitations, a numerical limit algorithm is at best a complement to a symbolic algorithm. The point, at this moment, is just to test N, although providing a numerical limit function in SymPy would also be a good idea.

I should note that there are much more sophisticated algorithms for numerical limits than the brute force method described here. Such algorithms are necessary to use especially when evaluating limits of indexed sequences where each element is expensive to compute (e.g. for summation of infinite series). A few acceleration methods for sequences and series are available in mpmath.