使用 Cauchy (CIF) 进行数值微分

Question

我正在尝试为泰勒级数创建一个带有数学类的模块，以便其他项目可以轻松访问它。因此，我希望尽可能优化它。

对于不太熟悉泰勒级数的人来说，能够多次在一个点上对一个函数进行微分是很有必要的。鉴于函数的数学导数的正常定义对于高阶导数需要极高的精度，我决定改用Cauchy's integral formula。通过一些工作，我设法重新排列了公式，正如您在这张图片中看到的那样：Rearranged formula。这为我提供了比传统导数定义更准确的高阶导数结果。这是我目前用来区分一个点的函数：

def myDerivative(f, x, dTheta, degree):
  riemannSum = 0 
  theta = 0

  while theta < 2*np.pi:
      functionArgument = np.complex128(x + np.exp(1j*theta))
      secondFactor = np.complex128(np.exp(-1j * degree * theta))
      riemannSum += f(functionArgument) * secondFactor * dTheta
      theta += dTheta
  return factorial(degree)/(2*np.pi) * riemannSum.real

我已经在我的主函数中测试了这个函数，并使用了一个经过深思熟虑的数学函数，我知道它的导数，即 f(x) = sin(x)。

def main():
  print(myDerivative(f, 0, 2*np.pi/(4*4096), 16))
  pass

这些导数似乎在 16 次导数附近吓坏了。我也尝试过使用 dTheta，但没有运气。我也想要更高的订单，但我担心我遇到了某种机器精度。

我的问题是最简单的形式：我可以做些什么来改进这个函数以获得更高阶的导数？

Answer 1

我似乎想出了解决问题的办法。我通过以不同的方式重新排列柯西积分公式来做到这一点，通过利用初始轮廓积分可以是围绕微分点的任意大圆。请注意，该函数必须在复平面中解析才能使其有效。

New formula

这也提供了一个新的微分功能：

def myDerivative(f, x, dTheta, degree, contourRadius):
  riemannSum = 0 
  theta = 0

  while theta < 2*np.pi:
      functionArgument = np.complex128(x + contourRadius*np.exp(1j*theta))
      secondFactor = (1/contourRadius)**degree*np.complex128(np.exp(-1j * degree * theta))
      riemannSum += f(functionArgument) * secondFactor * dTheta
      theta += dTheta

  return factorial(degree) * riemannSum.real / (2*np.pi)

这让我可以非常准确地区分高阶。例如，我可以毫无问题地区分 f(x)=e^x 50 次。

Answer 2

好吧，由于您使用的是导数的离散近似值（通过 dTheta），因此迟早您一定会遇到麻烦。我很惊讶你能够得到至少 15 个准确的导数——干得好！但是要获得所有订单的导数，要么你必须限制你愿意接受的东西并说它足够好，要么象征性地计算导数。看看Sympy。 Sympy 可能也有一些计算泰勒级数的函数。