Search code examples
pythonmathderivative

Derivatives in python

I am trying to find the coefficients of a finite series, $f(x) = \sum_n a_nx^n$. To get the $m$th coefficient, we can take the $m$th derivative evaluated at zero. Therefore, the $m$th coefficient is

$$
a_n = \frac{1}{2\pi i } \oint_C \frac{f(z)}{z^{n+1}} dz
$$

I believe this code takes the derivative of a function using the above contour integral.

import math
import numpy

import matplotlib.pyplot as plt

def F(x):
    mean=10
    return math.exp(mean*(x.real-1))

def p(n):
    mean=10
    return (math.pow(mean, n) * math.exp(-mean)) / math.factorial(n)

def integration(func, a, n, r, n_steps):
    z = r * numpy.exp(2j * numpy.pi * numpy.arange(0, 1, 1. / n_steps))
    return math.factorial(n) * numpy.mean(func(a + z) / z**n)

ns = list(range(20))
f2 = numpy.vectorize(F)

plt.plot(ns,[p(n) for n in ns], label='Actual')
plt.plot(ns,[integration(f2, a=0., n=n, r=1., n_steps=100).real/math.factorial(n) for n in ns], label='Numerical derivative')

plt.legend()

enter image description here

However, it is clear that the numerical derivative is completely off the actual values of the coefficients of the series. What am I doing wrong?

Solution

  • The formulas in the Mathematics Stack Exchange answer that you're using to derive the coefficients of the power series expansion of F are based on complex analysis - coming for example from Cauchy's residue theorem (though other derivations are possible). One of the assumptions necessary to make those formulas work is that you have a holomorphic (i.e., complex differentiable) function.

    Your definition of F gives a function that's not holomorphic. (For one thing, it always gives a real result for any complex input, which isn't possible for a non-constant holomorphic function.) But it's easily fixed to be holomorphic, while continuing to return the same result for real inputs.

    Here's a fixed version of F, which replaces x.real with x. Since the input to exp is now complex, it's also necessary to use cmath.exp instead of math.exp to avoid a TypeError:

    import cmath

def F(x):
    mean=10
    return cmath.exp(mean*(x-1))

    After that fix for F, if I run your code I get rather surprisingly accurate results. Here's the graph that I get. (I had to print out the values to double check that that graph really did show two lines on top of one another.)

    enter image description here