CODEWITHSUNDEEP
pythonfunctionnumpyscipy
Abdallah
Abdallah

Reputation: 81

Find root of a function in a given interval

I am trying to find the root of a function between by [0, pi/2], all algorithms in scipy have this condition : f(a) and f(b) must have opposite signs. In my case f(0)*f(pi/2) > 0 is there any solution, I precise I don't need solution outside [0, pi/2].

The function:

def dG(thetaf,psi,gamma) :
   return 0.35*((cos(psi))**2)*(2*sin(3*thetaf/2+2*gamma)+(1+4*sin(gamma)**2)*sin(thetaf/2)-‌​sin(3*thetaf/2))+(sin(psi)**2)*sin(thetaf/2)

Upvotes: 2

Views: 1866

Answers (2)

Saullo G. P. Castro
Saullo G. P. Castro

Reputation: 58965

Based on the comments and on @Mike Graham's answer, you can do something that will check where the change of signs are. Given y = dG(x, psi, gamma):

x[y[:-1]*y[1:] < 0]

will return the positions where you had a change of sign. You can an iterative process to find the roots numerically up to the error tolerance that you need:

import numpy as np
from numpy import sin, cos

def find_roots(f, a, b, args=[], errTOL=1e-6):
    err = 1.e6
    x = np.linspace(a, b, 100)
    while True:
        y = f(x, *args)
        pos = y[:-1]*y[1:] < 0
        if not np.any(pos):
            print('No roots in this interval')
            return roots
        err = np.abs(y[pos]).max()
        if err <= errTOL:
            roots = 0.5*x[:-1][pos] + 0.5*x[1:][pos]
            return roots
        inf_sup = zip(x[:-1][pos], x[1:][pos])
        x = np.hstack([np.linspace(inf, sup, 10) for inf, sup in inf_sup])

Upvotes: 1

Mike Graham
Mike Graham

Reputation: 76753

There is a root only if, between a and b, there are values with different signs. If this happens there are almost certainly going to be multiple roots. Which one of those do you want to find?

You're going to have to take what you know about f to figure out how to deal with this. If you know there is exactly one root, you can just find the local minimumn. If you know there are two, you can find the minimum and use that's coordinate c to find one of the two roots (one between a and c, the other between c and what used to be called b).

You need to know what you're looking for to be able to find it.

Upvotes: 1

Related Questions