Search code examples
pythoninterpolationroot

Find all roots of an arbitrary interpolated function in a given interval

I'm trying to find all roots of an arbitrary function within a defined interval. I already search for a solution, but didn't find anything that could solve my problem easily.

So for example I have got points and I use interp1d to create an function of these. Let's say a cosinus function:

x_p = np.linspace(0,2*pi,100)
y_p = np.cos(x_p)

Next, I use interp1d to generate the desired function

f = interp1d(x_p,y_p,kind="cubic")

I already try to use fsolve, but it only finds one root depending on the start point

fsolve(lambda x: f(x),0.1)

When I try

brentq(f,0,6)

I get the error

f(a) and f(b) must have different signs

I guess brentq is not the correct answer, since within the interval I got more than 1 roots. Is there an elegant solution for this problem?

Thank you very much in advance.

Solution

  • In general, finding all roots of arbitrary functions is impossible (See 1 & 2).

    For interpolation functions, they are usually constructed by concatenating polynomials, which one can determine the number of roots and find all roots. But you have to know the polynomial coefficients under interp1d in order to solve the problem.

    In reality, an unknown function is sampled to obtain a series of data points. In your case, the samples of the unknown function is given by cos(). The interpolant function f is used to simulate the behavior. The minimum spacing of discoverable roots is limited by sampling interval. This is because the existence theorem of root-finding tells that if f(a) * f(b) <= 0, then there is a root between [a, b]. In other words, if the sign does not change at f(a) and f(b), we have no conclusion about roots between [a,b]. This explains the error f(a) and f(b) must have different signs.

    For fsolve(lambda x: f(x),0.1), I got

    ValueError: A value in x_new is above the interpolation range.

    This mean the initial guess 0.1 is bad and the algorithm behind probably run x outside the interval [0, 2*pi].

    To make a good educated guess, it should be as close to the root as possible. Therefore, the left / right points where the sign change should be chosen.

    f_p = f(x_p)
indices, = np.nonzero(f_p[:-1] * f_p[1:] <= 0)
for i in range(indices.shape[0]):
      guess = x_p[indices[i]]
      print('root', i + 1, 'x=', fsolve(f, guess))

    Output:

    root 1 x= [1.57079633]
root 2 x= [4.71238898]