Find the root of a function in a given interval
I am trying to find the root of a function between [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
there is any solution, I definitely don't need a solution outside [ 0, pi/2]
.
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)
source to share
Based on the comments and answer by @Mike Graham, you can do something that will check where the sign change occurs. Given y = dG(x, psi, gamma)
:
x[y[:-1]*y[1:] < 0]
will return the positions where you had a sign change. You can iterate through the process to find the roots numerically, up to an acceptable deviation:
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])
source to share
There is a root only if between a
and b
there are meanings with different signs. If this happens, there will almost certainly be multiple roots. Which one do you want to find?
You need to take what you know about f
in order to figure out how to deal with it. If you know that there is exactly one root, you can simply find the local minimum of n. If you know there are two of them, you can find the minimum and use the coordinate c
to find one of the two roots (one between a
and c
, the other between c
and what was previously called b
).
You need to know what you are looking for in order to find it.
source to share