Calculating 1/r*d/dr(r*f) numerically in python when r=0. f is a function of r

Question

Usually when you do this by hand there's no problem as the 1/r usually gets cancelled with another r. But doing this numerically with scipy.misc.derivative works like a charm for r different from zero. But of course, as soon as I ask for r = 0, I get division by zero, which I expected. So how else could you calculate this numerically. I insist on the fact that everything has to be done numerically as my function are now so complicated that I won't be able to find a derivative manually. Thank you!

My code:

rAtheta = lambda _r: _r*Atheta(_r,theta,z,t)
if r != 0:
    return derivative(rAtheta,r,dx=1e-10,order=3)/r
else:
    #What should go here so that it doesn't blow up when calculating the gradient?

Answer 1

tl;dr: use symbolic differentiation, or complex step differentiation if that fails

If you insist on using numerical methods, you really have to approximate the limit of the derivative as r->0 one way or the other.

I suggest trying complex step differentiation. The idea is to use complex arguments inside the function you're trying to differentiate, but it usually gets rid of the numerical instability that is imposed by standard finite difference schemes. The result is a procedure that needs complex arithmetic (hooray numpy, and python in general!) but in turn can be much more stable at small dx values.

Here's another point: complex step differentiation uses

F′(x0) = Im(F(x0+ih))/h + O(h^2)

Let's apply this to your r=0 case:

F′(0) = Im(F(ih))/h + O(h^2)

There are no singularities even for r=0! Choose a small h, possibly the same dx you're passing to your function, and use that:

def rAtheta(_r):
    # note that named lambdas are usually frowned upon
    return _r*Atheta(_r,theta,z,t)

tol = 1e-10
dr = 1e-12

if np.abs(r) > tol: # or math.abs or your favourite other abs
    return derivative(rAtheta,r,dx=dr,order=3)/r
else:
    return rAtheta(r + 1j*dr).imag/dr/r

Here is the above in action for f = r*ln(r):

The result is straightforwardly smooth, even though the points below r=1e-10 were computed with complex step differentiation.

Very important note: notice the separation between tol and dr in the code. The former is used to determine when to switch between methods, and the latter is used as a step in complex step differentiation. Look what happens when tol=dr=1e-10:

the result is a smoothly wrong function below r=1e-10! That's why you always have to be careful with numerical differentiation. And I wouldn't advise going too much below that in dr, as machine precision will bite you sooner or later.

But why stop here? I'm fairly certain that your functions could be written in a vectorized way, i.e. they could accept an array of radial points. Using complex step differentiation you don't have to loop over the radial points (which you would have to resort to using scipy.misc.derivative). Example:

import numpy as np
import matplotlib.pyplot as plt

def Atheta(r,*args):
    return r*np.log(r) # <-- vectorized expression

def rAtheta(r):
    return r*Atheta(r) #,theta,z,t) # <-- vectorized as much as Atheta is

def vectorized_difffun(rlist):
    r = np.asarray(rlist)
    dr = 1e-12

    return (rAtheta(r + 1j*dr)).imag/dr/r

rarr = np.logspace(-12,-2,20)
darr = vectorized_difffun(rarr)
plt.figure()
plt.loglog(rarr,np.abs(darr),'.-')
plt.xlabel(r'$r$')
plt.ylabel(r'$|\frac{1}{r} \frac{d}{dr}(r^2 \ln r)|$')
plt.tight_layout()

plt.show()

The result should be familiar:

---

Having cleared the fun weirdness that is complex step differentiation, I should note that you should strongly consider using symbolic math. In cases like this when 1/r factors disappear exactly, it wouldn't hurt if you reached this conclusion exactly. After all double precision is still just double precision.

For this you'll need the sympy module, define your function symbolically once, differentiate it symbolically once, turn your simplified result into a numpy function using sympy.lambdify, and use this numerical function as much as you need (assuming that this whole process runs in finite time and the resulting function is not too slow to use). Example:

import sympy as sym

# only needed for the example:
import numpy as np
import matplotlib.pyplot as plt

r = sym.symbols('r')
f = r*sym.ln(r)
df = sym.diff(r*f,r)
res_sym = sym.simplify(df/r)
res_num = sym.lambdify(r,res_sym,'numpy')

rarr = np.logspace(-12,-2,20)
darr = res_num(rarr)

plt.figure()
plt.loglog(rarr,np.abs(darr),'.-')
plt.xlabel(r'$r$')
plt.ylabel(r'$|\frac{1}{r} \frac{d}{dr}(r^2 \ln r)|$')    
plt.tight_layout()
plt.show()

resulting in

As you see, the resulting function was vectorized thanks to lambdify using numpy during the conversion from symbolic to numeric function. Obviously, the best solution is the symbolic one as long as the resulting function is not so complicated to make its practical use impossible. I urge you to first try the symbolic version, and if for some reason it's not applicable, switch to complex step differentiation, with due caution.