0

I am having a problem solving for the variable x_Prime using the scipy.optimize.fsolve() function. Below is what I have:

import math
import scipy.optimize as so

c_val = 100
f_NFW_val = math.log(1+c_val)-c_val/(1+c_val)

psi_prime = 3.0608421399604424

eqn = lambda x_Prime: psi_prime  -  1/f_NFW_val * ( math.log(1+c_val*x_Prime)/x_Prime - c_val/(1+c_val*x_Prime) )
sol = so.fsolve(eqn, 1)[0]

I get the error: Error code ouput

It seems like there is no solution for psi_prime = 3.0608421399604424. However, if I try plotting the curve using the following code:

c_val = 100
f_NFW_val = math.log(1+c_val)-c_val/(1+c_val)

psi_prime = 3.0608421399604424

xs = np.linspace(0.1633804055843348, 1, 100)
plt.plot(xs, 1/f_NFW_val * ( np.log(1+c_val*xs)/xs - c_val/(1+c_val*xs) ))
plt.plot(xs, psi_prime*np.ones(np.size(xs)), color="grey", linestyle="dashed")
plt.show()

Error graph

So, there should be a solution. Can someone tell me why fsolve() is not able to solve the equation?

ellipse314
  • 47
  • 3

1 Answers1

2

The domain error is caused by an invalid argument of the logarithm. You should use numpy.log to handle cases like this (it will return NaN instead of raising an exception).

The second point is that you use an approach that is likely to diverge and stuck somewhere around the initial point. You should try different methods and analyze the function graph to come up with a way to find the root numerically, there is just no universal method to do so. In your case a default method from scipy.root_scalar should work just fine.

eqn = lambda x_Prime: psi_prime  -  1/f_NFW_val * ( np.log(1+c_val*x_Prime)/x_Prime - c_val/(1+c_val*x_Prime) )
sol = so.root_scalar(eqn, bracket=(0, 2))
print(sol)

It converges at 0.17997695477705547

root_scalar documentation

Kuroneko
  • 146
  • 6
  • Thanks for your answer! Quick question. What is the tuple assigned to the `bracket` argument. Is it just that the required answer lies between the two values of the tuple? In general, I don't know where the required answer would lie. – ellipse314 Mar 01 '23 at 14:21
  • Sorry. Just read the documentation. I guess the tuple is the expected interval in which the answer lies. My code is automated. In general, my `x_Prime` lies any anywhere between `0.1633804055843348` and `infinity`. But if I give the second value of the tuple as `np.inf`, I get the answer as `nan`. – ellipse314 Mar 01 '23 at 14:27
  • I played around with it in Python. I set the second value of the tuple to an insanely high value like `1e20`. Python did not complain when I automated the code. – ellipse314 Mar 01 '23 at 14:55