CODEWITHSUNDEEP
python-3.xnumpymathdifferentiation
ambitphysx
ambitphysx

Reputation: 51

Numerical differentiation using Cauchy (CIF)

I am trying to create a module with a mathematical class for Taylor series, to have it easily accessible for other projects. Hence I wish to optimize it as far as I can.

For those who are not too familiar with Taylor series, it will be a necessity to be able to differentiate a function in a point many times. Given that the normal definition of the mathematical derivative of a function will require immense precision for higher order derivatives, I've decided to use Cauchy's integral formula instead. With a little bit of work, I've managed to rearrange the formula a little bit, as you can see on this picture: Rearranged formula. This provided me with much more accurate results on higher order derivatives than the traditional definition of the derivative. Here is the function i am currently using to differentiate a function in a point:

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

I've tested this function in my main function with a carefully thought out mathematical function which I know the derivatives of, namely f(x) = sin(x).

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

These derivatives seems to freak out at around the derivative of degree 16. I've also tried to play around with dTheta, but with no luck. I would like to have higher orders as well, but I fear I've run into some kind of machine precission.

My question is in it's simplest form: What can I do to improve this function in order to get higher order of my derivatives?

Upvotes: 3

Views: 822

Answers (2)

ambitphysx
ambitphysx

Reputation: 51

I seem to have come up with a solution to the problem. I did this by rearranging Cauchy's integral formula in a different way, by exploiting that the initial contour integral can be an arbitrarily large circle around the point of differentiation. Be aware that it is very important that the function is analytic in the complex plane for this to be valid.

New formula

Also this gives a new function for differentiation:

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)

This gives me a very accurate differentiation of high orders. For instance I am able to differentiate f(x)=e^x 50 times without a problem.

Upvotes: 2

Robert Dodier
Robert Dodier

Reputation: 17595

Well, since you are working with a discrete approximation of the derivative (via dTheta), sooner or later you must run into trouble. I'm surprised you were able to get at least 15 accurate derivatives -- good work! But to get derivatives of all orders, either you have to put a limit on what you're willing to accept and say it's good enough, or else compute the derivatives symbolically. Take a look at Sympy for that. Sympy probably has some functions for computing Taylor series too.

Upvotes: 0

Related Questions