当 r=0 时，在 python 中数值计算 1/r*d/dr(r*f)。 f 是 r 的函数
Calculating 1/r*d/dr(r*f) numerically in python when r=0. f is a function of r
通常当您手动执行此操作时没有问题，因为 1/r 通常会被另一个 r 取消。但是用 scipy.misc.derivative 以数字方式执行此操作就像 r 不同于零的魅力一样。但是，当然，只要我要求 r = 0，我就会得到除以零的结果，这是我所期望的。那么你怎么能用数字计算这个呢。我坚持认为所有事情都必须用数字来完成，因为我的函数现在太复杂了，我无法手动找到导数。谢谢！
我的代码：
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?
tl;dr：使用符号微分，如果失败则使用复步微分
如果你坚持使用数值方法，你真的必须以一种或另一种方式逼近导数的极限r->0。
我建议尝试 complex step differentiation。这个想法是在你试图微分的函数内部使用复杂的参数，但它通常会消除标准有限差分方案强加的数值不稳定性。结果是一个需要复杂算术的过程（hooray numpy，通常 python！）但反过来在小 dx 值时可以更加稳定。
这里还有一点：复步微分使用
F′(x0) = Im(F(x0+ih))/h + O(h^2)
让我们将其应用到您的 r=0 案例中：
F′(0) = Im(F(ih))/h + O(h^2)
即使 r=0 也没有奇点！选择一个小 h，可能与您传递给函数的 dx 相同，然后使用它：
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
这是 f = r*ln(r) 的上述操作：
即使 r=1e-10 以下的点是使用复杂的步长微分计算的，结果还是非常平滑的。
非常重要的注意事项：请注意代码中 tol 和 dr 之间的分隔。前者用于确定何时在方法之间切换，后者用作复杂步骤分化中的一个步骤。看看 tol=dr=1e-10:
时会发生什么
结果是r=1e-10下面的顺利错函数！这就是为什么您始终必须小心数值微分的原因。而且我不建议在 dr 中低于太多，因为机器精度迟早会咬你。
但为什么要到此为止呢？我相当确定您的函数可以用矢量化的方式编写，即它们可以接受一组径向点。使用复杂的步长微分，您不必遍历径向点（您必须求助于使用 scipy.misc.derivative）。示例：
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()
结果应该很熟悉：
清除了复杂步骤微分的有趣怪癖后，我应该注意到您应该认真考虑使用符号数学。在这种 1/r 因素完全消失的情况下，如果您准确地得出这个结论也不会受到伤害。毕竟双精度仍然只是双精度。
为此，您需要 sympy module, define your function symbolically once, differentiate it symbolically once, turn your simplified result into a numpy function using sympy.lambdify，并根据需要使用 this 数值函数（假设整个过程在有限时间内运行，并且生成的函数使用起来不会太慢）。示例：
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()
导致
如您所见，由于 lambdify 在从符号函数到数字函数的转换过程中使用 numpy，结果函数被矢量化。显然，最好的解决方案是符号解决方案，只要生成的函数不会复杂到无法实际使用的程度。我敦促您首先尝试符号版本，如果由于某种原因它不适用，请谨慎切换到复步微分。
通常当您手动执行此操作时没有问题，因为 1/r 通常会被另一个 r 取消。但是用 scipy.misc.derivative 以数字方式执行此操作就像 r 不同于零的魅力一样。但是，当然，只要我要求 r = 0，我就会得到除以零的结果，这是我所期望的。那么你怎么能用数字计算这个呢。我坚持认为所有事情都必须用数字来完成，因为我的函数现在太复杂了，我无法手动找到导数。谢谢！
我的代码：
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?
tl;dr：使用符号微分，如果失败则使用复步微分
如果你坚持使用数值方法，你真的必须以一种或另一种方式逼近导数的极限r->0。
我建议尝试 complex step differentiation。这个想法是在你试图微分的函数内部使用复杂的参数，但它通常会消除标准有限差分方案强加的数值不稳定性。结果是一个需要复杂算术的过程（hooray numpy，通常 python！）但反过来在小 dx 值时可以更加稳定。
这里还有一点：复步微分使用
F′(x0) = Im(F(x0+ih))/h + O(h^2)
让我们将其应用到您的 r=0 案例中：
F′(0) = Im(F(ih))/h + O(h^2)
即使 r=0 也没有奇点！选择一个小 h，可能与您传递给函数的 dx 相同，然后使用它：
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
这是 f = r*ln(r) 的上述操作：
即使 r=1e-10 以下的点是使用复杂的步长微分计算的，结果还是非常平滑的。
非常重要的注意事项：请注意代码中 tol 和 dr 之间的分隔。前者用于确定何时在方法之间切换，后者用作复杂步骤分化中的一个步骤。看看 tol=dr=1e-10:
结果是r=1e-10下面的顺利错函数！这就是为什么您始终必须小心数值微分的原因。而且我不建议在 dr 中低于太多，因为机器精度迟早会咬你。
但为什么要到此为止呢？我相当确定您的函数可以用矢量化的方式编写，即它们可以接受一组径向点。使用复杂的步长微分，您不必遍历径向点（您必须求助于使用 scipy.misc.derivative）。示例：
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()
结果应该很熟悉：
清除了复杂步骤微分的有趣怪癖后，我应该注意到您应该认真考虑使用符号数学。在这种 1/r 因素完全消失的情况下，如果您准确地得出这个结论也不会受到伤害。毕竟双精度仍然只是双精度。
为此，您需要 sympy module, define your function symbolically once, differentiate it symbolically once, turn your simplified result into a numpy function using sympy.lambdify，并根据需要使用 this 数值函数（假设整个过程在有限时间内运行，并且生成的函数使用起来不会太慢）。示例：
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()
导致
如您所见，由于 lambdify 在从符号函数到数字函数的转换过程中使用 numpy，结果函数被矢量化。显然，最好的解决方案是符号解决方案，只要生成的函数不会复杂到无法实际使用的程度。我敦促您首先尝试符号版本，如果由于某种原因它不适用，请谨慎切换到复步微分。