I wonder how I could manage to implement an iterative root finder using fsolve over an interval, up until it found N roots ? (Assuming I know how small the steps should be to get every roots during the procedure) Is there a way to do so with a simple double for loop ?
Here is what it would look like for a "simple" function (cos(x)*x):
import numpy as np
from scipy.optimize import fsolve
import numpy as np
def f(x):
return np.cos(x)*x
for n in range(1,10) :
a = 0
k = 0
while k < 1000 :
k = fsolve(f,a)
if k == a :
a = a 0.01
k = fsolve(f,a)
else :
print(k)
But I can't make it works this way. I can't use chebpy because my real function is more complexe (involving bessel function) and chepby doesnt seem to accept such function as an argument.
Edit : Corrected indentation , this program yields 0 (first solution) an infinite number of time without stopping.
CodePudding user response:
can you share your error? may be something related to the indentation of your f(x) function, try changing your code to this:
def f(x):
return np.cos(x)*x
CodePudding user response:
I found a solution that does the job as of now. It consist of passing an array corresponding to the search interval, and then sorting out the solutions (in my case, only looking at positive solutions, removing duplicates etc etc)
It might not be the best way but it works for me :
In this example i'm looking for the 10 first positive solutions of cos(x)*x=0 assuming they would be in [0,100].
import numpy as np
from scipy.optimize import fsolve
def f(x):
return np.cos(x)*x
int = np.arange(0,100,1)
arr = np.array([int])
roots=fsolve(f,arr)
# print(roots)
roots=np.around(roots, decimals=5, out=None)
a = roots[roots >= 0]
b = np.unique(a)
b=(b[:10])
print(b)
Result :
[ 0. 1.5708 4.71239 7.85398 10.99557 14.13717 17.27876 20.42035
23.56194 26.70354]
I had to use np.around otherwise np.unique would not works. Thanks again.