Cubic Equation solver in a loop: Python

286 views Asked by At

I'm working on a van der waals equation of state calculator that that imports critical properties data from a spreadsheet, and uses those properties to calculate the volume and compressibility factor of a substance specified by the user. So far the code can read the sheet, retrieve the values and perform basic calculations with them. I'm running into trouble with solving the cubic equation for the compressibility factor and using the solution to solve for volume. It is important to note that I'm using a series of loops for the pressure and temperature so that I have a range of compressibility values, and a range of volumes across various conditions.

Issues: The cubic equation solver will return a list of imaginary roots, and the I'm having a hard time setting my volume equation equal to my list of compressibility values to solve for volume.

I'm totally new to python, any help is appreciated!

from sympy import var,solve,Eq
import sympy 
y = var('y', real=True)
x = var('x', real=True)
VANDERWALS
#Finds critical values based on substance
value = df[df[0] == name][1].values[0]
value2 = df[df[0] == name][2].values[0]
value3 = df[df[0] == name][3].values[0]
value4 = df[df[0] == name][4].values[0]
#Solves for compressability factor
pc = value3*10 #conversion to bar
R = 8.314e-5
a = (27*(pow(R,2))*(pow(value2,2)))/(64*pc)
b = (R*value2)/(8*pc)
print(a)
print(b)
#Solves of volume across require temp and pressure range
for tloop in range(7):
    tempval = (tloop-2)*10
    temp = tempval+273
    for ploop in range(30):
        presval = (ploop+1)
        R = 8.314E-5
        A = (a*presval)/pow(R*temp,2)
        B = (b*presval)/(R*temp)
        eq1 = Eq((pow(y,3))-(pow(y,2)*(1+B))+(A*y)-(A*B),0)
        sol = solve(eq1)
        print(sol[0])
        eq2= Eq(((presval*x)/(R*temp)),sol)
        sol2 = solve(eq2)
        
1

There are 1 answers

0
smichr On

real_roots(eq) will give a list of real roots; it works best if all numbers are Rationals or nsolve(eq, guess-for-solution). In your case, in a loop, you will have the last known solution so maybe you could use that as your initial guess for the next set of parameters.