Using fsolve in Python

Question

I want to use fsolve to find the 2 values of a function z(q) in which I know the value of q (q_1000)

My code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
k_b = 1.38064852e-23 
n_0 = 2e+25 
m = 4.8e-26 
g = 9.81 
T = 300 
I_inf = 400 
sigma = 4e-22 
kappa = 0.2 
H = (k_b*T)/(m*g) 
val = 100
z_min = 1e+5
z_max = 4e+5
z = np.linspace(z_min, z_max, val)
# ______________________________ Functions ___________________________________________
n = lambda z : n_0*np.exp(-z/H) 
I = lambda z : I_inf*(np.exp(-(sigma*n(z)*H))) 
q = lambda z : kappa*sigma*n(z)*I(z) 
q_max = q(z).argmax() # q_max
q_1000 = q(z)[q_max]/1000 # q_max / 1000
result1 = fsolve(q,q_1000) # start at z = q_1000
plt.xscale('log')
plt.yscale('log')
plt.title('z(q)')
plt.axis([1e-11, 1e-1, 1e+5, 4e+5])
plt.plot(q(z)[q_max], z[q_max],'x', label="z*=q_max", color="red")
plt.plot (q(z),z,'k-')
plt.plot(result1,q_1000,'ro')
plt.legend()
plt.xlabel('q')
plt.ylabel('z')
plt.show()
print(q(z)[q_max])
print(q_1000)
print(result1)

That's my result

What I have

That's what I want, with the 2 points (let's say Z1 and Z2) What I want

• Peter O.: Include imports in your code. I have added some of them.
• Michel Chance: Ok, thank you, that's also "scipy.optimize"

Answer 1

The first parameter to fsolve should be a function for which the roots q(z) = 0 are sought. Therefore, q(z) = q_1000 should be written as q(z) - q_1000. The second parameter to fsolve is an approximation to the desired root. I'm not sure how to get a good initial value in general, but in this case the plot suggests 1.5e5 and 2.5-e5.

In code, including the plotting:

z1 = fsolve(lambda z: q(z) - q_1000,1.5e5)[0]
z2 = fsolve(lambda z: q(z) - q_1000,2.5e5)[0]
plt.plot((q_1000, q_1000), (z1, z2), linestyle='dashed', color='dodgerblue', label='z*=q_1000')
plt.plot(q_1000, z1,'ro')
plt.plot(q_1000, z2,'ro')

the resulting plot