numerical root finding for positive definite function in Python

Question

I have a very complicated positive definite linear continuous function of a single variable k for which I am trying to find all roots in a given range of k; say -4 < k < 4.

Up until now I have first estimated the minima of the function, by searching for points k_j where both $k_{j+1}>k_j$ and $k_{j-1}>k_j$. Then using each of these points as a starting point I apply the optimisation function, scipy.optimize.newton. To some extent this method has worked. However, as my functions get more complicated the search for minima becomes more and more time consuming and possibly inaccurate.

Is there any built in function in numpy or scipy which searches in a given domain (eg -4 < k < 4) of a function and finds all the roots. I am willing to sacrifice some computational efficiency so that I do not have to specify exact points to search near.

Thanks

Answer 1

You can use range:

k_list = range(-4, 4)

But that only does integers, the problem here is specifying your step. Obviously there are infinite decimals between -4 and 4 so you need to specify how many decimals you want to go to.

You can use numpy.arange in order to make a list from a range and set the incrementing value

For example

k_list = numpy.arange(-4, 4, 0.5)

will increment by 0.5

>>> numpy.arange(-4, 4, 0.5)
>>> [-4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4]

if you want to increment by a smaller amount and therefore get a bigger range of values then make the 0.5 smaller.

You will have to specify an increment because there's infinite decimals in that range, as previously stated.

After specifying your list you can then iterate through your list with your function to find the roots.

for k in k_list:
    some_function(k)
    return root

EDIT:

For this to work you will of course need a function that finds the root of k, however if I've understood your question correctly this should just be your linear equation, to use a simple example : root = 2k (the mathematical way of writing this is of course y=2x.

For the sake of simplicity let's just say your function is y=2x at which point your script would become

k_list = numpy.arange(-4, 4, 0.5)

for k in k_list:
    root = 2*k
    return root

And then you just specify your own value for 0.5 to decide to what decimal your k values are

Unless of course you are looking at a form of quadratic. In which case we may have something like

y = x^2 - 2x +2

This makes your question slightly more confusing. You'd obviously find the root of x by setting y=0 however, now you have one variable which I imagine is what you mean by k, which you specify, leaving a sum rather than a formula.

In this case I'd let y=k and then specify your k value and solve to find your root.

eg:

y = 32x^2 - 7.2x + 12
let y = k
k = 32x^2 - 7.2x +12
let k = -4 (we'd iterate through your range in reality)
4 = 32x^2 - 7.2x + 12
0 = 32x^2 - 7.2x + 8
and solve for x (aka. root)

I expect there is a numpy or scipy method for solving formulas with multiple instances of the same variable. I am not an expert in either lib though so I wouldn't be able to advise you in that direction.

See also: http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html#numpy.polyval