finding roots of an expression y = f(x) for different values of y using scipy

Question

I have the following function.

import numpy as np
from scipy import optimize

def func(x):
    theta = np.arccos(1 - 2*x)
    return 2*(theta - np.sin(theta))

Question: I want to find the roots of the expression y = func(x) for different values of y. For example if y =1 I can modify the above function as

import numpy as np
from scipy import optimize

def func(x):
    theta = np.arccos(1 - 2*x)
    return 2*(theta - np.sin(theta)) - 1

initial_guess = 0
sol = optimize.root(func, initial_guess)

To obtain the desired root. Let us assume that such roots exist and we know appropriate initial guesses for all the different y values.

What modifications do I need to achieve this?

Answer 1

To vary the left hand side you could add an extra argument to your func and then pass this via the args keyword of scipy.optimize.root (keeping rths answer in mind) :

import numpy as np
from scipy import optimize

def func(x, lhs):
    theta = np.arccos(1 - 2*x)
    return 2*(theta - np.sin(theta)) -  lhs

initial_guess = 0.001
lhs = 1
sol = optimize.root(func, initial_guess, args=(lhs,))

Answer 2

The problem with your initial_guess, which is just the point of the singularity of your function. Therefore if you shift initial_guess a bit you can obtain a nice solution....

>>> initial_guess = 0.001                                                                                                                                                                                                                                
>>> sol = optimize.root(func, initial_guess)
>>> sol                                                                                                                                                                                                                                                                         
    fjac: array([[-1.]])                                                                                                                                                                                                                                                        
     fun: array([ -2.22044605e-16])                                                                                                                                                                                                                                             
 message: 'The solution converged.'                                                                                                                                                                                                                                             
    nfev: 11                                                                                                                                                                                                                                                                    
     qtf: array([  2.99071878e-12])                                                                                                                                                                                                                                             
       r: array([-3.71631332])                                                                                                                                                                                                                                                  
  status: 1                                                                                                                                                                                                                                                                     
 success: True                                                                                                                                                                                                                                                                  
       x: array([ 0.46328511])                                                                                                                                                                                                                                                  
>>>   

NOTE. Due to problem singularity, you cannot obtain numerical solution if your initial guess is negative or the zero.