Python Forum

Full Version: how can I solve fsolve function error?
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
import matplotlib.pyplot as plt  
import numpy as np
from scipy.optimize import fsolve

R = 0.08206
Tc = 150.687
Pc = 47.994
ω = -0.00219
T = 100

a = 0.42748*R**2*Tc**2/Pc
b = 0.08664*R*Tc/Pc
α = (1 + (0.48 + 1.574*ω - 0.176*ω**2)*(1 - (T/Tc)**(0.5)))**2


X = np.linspace(1, 731, 100) 
k = [] 
for dx in X: 
    P=dx
    def f(V):
        return P - R*T/(V - b)  + a*α/ V/(V + b)
    k.append(fsolve(f,0.05)*dx/R/T)
plt.xlabel("Pressure (atm)")
plt.ylabel("Compressibility Z")
plt.title("Z-p plot at 100K, Ar")
plt.plot(X, k, label = 'SRK')
plt.show()
This is my code for SRK Equation of state.
For some reason, when the temperature is high enough, (ex:T = 300) it won't have any problem.
[Image: P4zwO.png]
But when the temperature is low, like the code that I mentioned above, it will start to print out at error saying [Image: 5CKTN.png]
I've changed the number around a bit for my initial guess, but it won't help either

What changes should I make to solve this problem?
Giving a negative reputation for announcing that you cross-posted without disclosing it is childish.
In the stackoverflow thread you state that the 100k curve should look a lot like the 300k curve. I don't think that is true. Lowering the temperature not only lowers the compressibility, it also changes the shape of the curve, making it much sharper. It might be that fsolve is giving up because the equation is too non-linear and you are bouncing around but not getting any closer to the solution, or it could be that there is no solution. I noticed I can make the warning go away by starting with a really small starting value for V (0.000001). The resulting curve looks like a diagonal line to me. No jaggies and no curve.