Calculus Knight
Calculus Knight

Reputation: 113

Python fsolve does not take array of floats. How to implement it?

I used fsolve to find the zeros of an example sinus function, and worked great. However, I wanted to do the same with a dataset. Two lists of floats, later converted to arrays with numpy.asarray(), containing the (x,y) values, namely 't' and 'ys'.

Although I found some related questions, I failed to implement the code provided in them, as I try to show here. Our arrays of interest are stored in a 2D list (data[i][j], where 'i' corresponds to a variable (e.g. data[0]==t==time==x values) and 'j' are the values of said variable along the x axis (e.g. data[1]==Force). Keep in mind that each data[i] is an array of floats.

Could you offer an example code that takes two inputs (the two mentioned arrays) and returns its intersecting points with a defined function (e.g. 'y=0').

I include some testing I made regarding the other related question. ( @HYRY 's answer)

I do not think it is relevant, but I'm using Spyder through Anaconda.

Thanks in advance!

"""
Following the answer provided by @HYRY in the 'related questions' (see link above).
At this point of the code, the variable 'data' has already been defined as stated before.
"""
from scipy.optimize import fsolve

def tfun(x):
    return data[0][x]

def yfun(x):
    return data[14][x]

def findIntersection(fun1, fun2, x0):
    return [fsolve(lambda x:fun1(x)-fun2(x, y), x0) for y in range(1, 10)]

print findIntersection(tfun, yfun, 0)

Which returns the next error

  File "E:/Data/Anaconda/[...]/00-Latest/fsolvestacktest001.py", line 36, in tfun
    return data[0][x]

IndexError: arrays used as indices must be of integer (or boolean) type

The full output is as it follows:

Traceback (most recent call last):

  File "<ipython-input-16-105803b235a9>", line 1, in <module>
    runfile('E:/Data/Anaconda/[...]/00-Latest/fsolvestacktest001.py', wdir='E:/Data/Anaconda/[...]/00-Latest')

  File "C:\Anaconda\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 580, in runfile
    execfile(filename, namespace)

  File "E:/Data/Anaconda/[...]/00-Latest/fsolvestacktest001.py", line 44, in <module>
    print findIntersection(tfun, yfun, 0)

  File "E:/Data/Anaconda/[...]/00-Latest/fsolvestacktest001.py", line 42, in findIntersection
    return [fsolve(lambda x:fun1(x)-fun2(x, y), x0) for y in range(1, 10)]

  File "C:\Anaconda\lib\site-packages\scipy\optimize\minpack.py", line 140, in fsolve
    res = _root_hybr(func, x0, args, jac=fprime, **options)

  File "C:\Anaconda\lib\site-packages\scipy\optimize\minpack.py", line 209, in _root_hybr
    ml, mu, epsfcn, factor, diag)

  File "E:/Data/Anaconda/[...]/00-Latest/fsolvestacktest001.py", line 42, in <lambda>
    return [fsolve(lambda x:fun1(x)-fun2(x, y), x0) for y in range(1, 10)]

  File "E:/Data/Anaconda/[...]/00-Latest/fsolvestacktest001.py", line 36, in tfun
    return data[0][x]

IndexError: arrays used as indices must be of integer (or boolean) type

Upvotes: 0

Views: 2483

Answers (1)

Lanting
Lanting

Reputation: 3068

You can 'convert' a datasets (arrays) to continuous functions by means of interpolation. scipy.interpolate.interp1d is a factory that provides you with the resulting function, which you could then use with your root finding algorithm. --edit-- an example for computing an intersection of sin and cos from 20 samples (I've used cubic spline interpolation, as piecewise linear gives warnings about the smoothness):

>>> import numpy, scipy.optimize, scipy.interpolate
>>> x = numpy.linspace(0,2*numpy.pi, 20)
>>> x
array([ 0.        ,  0.33069396,  0.66138793,  0.99208189,  1.32277585,
    1.65346982,  1.98416378,  2.31485774,  2.64555171,  2.97624567,
    3.30693964,  3.6376336 ,  3.96832756,  4.29902153,  4.62971549,
    4.96040945,  5.29110342,  5.62179738,  5.95249134,  6.28318531])
>>> y1sampled = numpy.sin(x)
>>> y2sampled = numpy.cos(x)
>>> y1int = scipy.interpolate.interp1d(x,y1sampled,kind='cubic')
>>> y2int = scipy.interpolate.interp1d(x,y2sampled,kind='cubic')
>>> scipy.optimize.fsolve(lambda x: y1int(x) - y2int(x), numpy.pi)
array([ 3.9269884])
>>> scipy.optimize.fsolve(lambda x: numpy.sin(x) - numpy.cos(x), numpy.pi)
array([ 3.92699082])

Note that interpolation will give you 'guesses' about what data should be between the sampling points. No way to tell how good these guesses are. (but for my example, you can see it's a pretty good estimation)

Upvotes: 1

Related Questions