Pythonic way to manage arbitrary amount of variables, used for equation solving.

Question

This is a bit difficult to explain without a direct example. So let's put the very simplistic ideal-gas law as example. For an ideal gas under normal circumstances the following equation holds:

PV = RT

This means that if we know 3 of the 4 variables (pressure, volume, specific gas constant and temperature) we can solve for the other one.

How would I put this inside an object? I want to have an object where I can just insert 3 of the variables, and then it calculates the 4th. I wonder if this can be achieved through properties?

My current best guess is to insert it like:

class gasProperties(object):
    __init__(self, P=None, V=None, R=None, T=None)
        self.setFlowParams(P, V, R, T)
    def setFlowParams(self, P=None, V=None, R=None, T=None)
        if P is None:
            self._P = R*T/V
            self._V = V
            self._R = R
            self._T = T
        elif V is None:
            self._V = R*T/P
            self._P = P
            self._R = R
            self._T = T
        #etc

Though this is quite cumbersome, and error prone (I have to add checks to see that exactly one of the parameters is set to "None").

Is there a better, cleaner way?

I see this "problem" happening quite often, in all kinds of various ways, and especially once the number of variables grows (adding density, reynolds number, viscosity to the mix) the number of different if-statements grows quickly. (IE if I have 8 variables and any 5 make the system unique I would need 8 nCr 5 = 56 if statements).

Answer 1

Using sympy, you can create a class for each of your equations. Create the symbols of the equation with ω, π = sp.symbols('ω π') etc., the equation itself and then use function f() to do the rest:

import sympy as sp    

# Create all symbols.
P, V, n, R, T = sp.symbols('P V n R T')
# Create all equations
IDEAL_GAS_EQUATION = P*V - n*R*T   

def f(x, values_dct, eq_lst):
    """
    Solves equations in eq_lst for x, substitutes values from values_dct, 
    and returns value of x.

    :param x: Sympy symbol
    :param values_dct: Dict with sympy symbols as keys, and numbers as values.
    """

    lst = []
    lst += eq_lst

    for i, j in values_dct.items():
        lst.append(sp.Eq(i, j))

    try:
        return sp.solve(lst)[0][x]
    except IndexError:
        print('This equation has no solutions.')

To try this out... :

vals = {P: 2, n: 3, R: 1, T:4}

r = f(V, values_dct=vals, eq_lst=[IDEAL_GAS_EQUATION, ])
print(r)   # Prints 6

If you do not provide enough parameters through values_dct you ll get a result like 3*T/2, checking its type() you get <class 'sympy.core.mul.Mul'>.

If you do provide all parameters you get as a result 6 and its type is <class 'sympy.core.numbers.Integer'>, so you can raise exceptions, or whatever you need. You could also, convert it to an int with int() (it would raise an error if instead of 6 you had 3*T/2 so you can test it that way too).

Alternatively, you can simply check if None values in values_dct are more than 1.

---

To combine multiple equations, for example PV=nRT and P=2m, you can create the extra symbol m like the previous symbols and assign 2m to the new equation name MY_EQ_2, then insert it in the eq_lst of the function:

m = sp.symbols('m')
MY_EQ_2 = P - 2 * m

vals = {n: 3, R: 1, T:4}

r = f(V, values_dct=vals, eq_lst=[IDEAL_GAS_EQUATION, MY_EQ_2])
print(r)   # Prints 6/m

Answer 2

A basic solution using sympy, and kwargs to check what information the user has provided:

from sympy.solvers import solve
from sympy import Symbol
def solve_gas_properties(**kwargs):
    properties = []
    missing = None
    for letter in 'PVRT':
        if letter in kwargs:
            properties.append(kwargs[letter])
        elif missing is not None:
            raise ValueError("Expected 3 out of 4 arguments.")
        else:
            missing = Symbol(letter)
            properties.append(missing)
    if missing is None:
        raise ValueError("Expected 3 out of 4 arguments.")
    P, V, R, T  = properties
    return solve(P * V - R * T, missing)

print solve_gas_properties(P=3, V=2, R=1) # returns [6], the solution for T

This could then be converted into a class method, drawing on class properties instead of keyword arguments, if you want to store and manipulate the different values in the system.

The above can also be rewritten as:

def gas_properties(**kwargs):
    missing = [Symbol(letter) for letter in 'PVRT' if letter not in kwargs]
    if len(missing) != 1:
        raise ValueError("Expected 3 out of 4 arguments.")
    missing = missing[0]
    P, V, R, T = [kwargs.get(letter, missing) for letter in 'PVRT']
    return solve(P * V - R * T, missing)

Answer 3

One solution could be the use of a dictionary to store variable names and their values. This allows you to easily add other variables at any time. Also, you can check that exactly one variable has value "None" by counting the number of "None" items in your dictionary.

Answer 4

My approach would be fairly simple:

class GasProperties(object):
    def __init__(self, P=None, V=None, R=None, T=None):
        self.setFlowParams(P, V, R, T)
    def setFlowParams(self, P=None, V=None, R=None, T=None):
        if sum(1 for arg in (P, V, R, T) if arg is None) != 1:
            raise ValueError("Expected 3 out of 4 arguments.")
        self._P = P
        self._V = V
        self._R = R
        self._T = T
    @property
    def P(self):
        return self._P is self._P is not None else self._R*self._T/self._V

You similarly define properties for V, R and T.

Answer 5

This approach allows you to set up object's attributes:

def setFlowParams(self, P=None, V=None, R=None, T=None):
    params = self.setFlowParams.func_code.co_varnames[1:5]
    if sum([locals()[param] is None for param in params]) > 1:
        raise ValueError("3 arguments required")
    for param in params:
        setattr(self, '_'+param, locals()[param])

In addition, you need to define getters for attributes with formulas. Like this:

@property
def P(self):
    if self._P is None:
        self._P = self._R*self._T/self._V
    return self._P

Or calculate all values in setFlowParams.

Answer 6

Numerical methods

You might want to do this without sympy, as and exercise for example, with numerical root finding. The beauty of this method is that it works for a extremely wide range of equations, even ones sympy would have trouble with. Everybody i know was taught this in uni on bachelor maths course*, unfortunately not many can apply this in practice.

So first we get the rootfinder you can find code examples on wikipedia and on the net at large this is fairly well known stuff. Many math packages have these built in see for example scipy.optimize for good root finders. I'm going to use the secant method for ease of implementation (in this case i don't really need iterations but ill use generic versions anyway if you happen to want to use some other formulas).

"""Equation solving with numeric root finding using vanilla python 2.7"""

def secant_rootfind(f, a, incr=0.1, accuracy=1e-15):
    """ secant root finding method """
    b=a+incr;
    while abs(f(b)) > accuracy :
        a, b = ( b, b - f(b) * (b - a)/(f(b) - f(a)) )

class gasProperties(object):
    def __init__(self, P=None,V=None,n=None,T=None):
        self.vars = [P, V, n, 8.314, T]
        unknowns = 0
        for i,v  in enumerate(self.vars):
            if v is None :
                self._unknown_=i
                unknowns += 1
        if unknowns > 1:
             raise ValueError("too many unknowns")

    def equation(self, a):
        self.vars[self._unknown_] = a
        P, V, n, R, T = self.vars
        return P*V - n*R*T # = 0

    def __str__(self):
        return str((
                   "P = %f\nV = %f\nn = %f\n"+
                   "R = %f\nT = %f ")%tuple(self.vars))

    def solve(self):
        secant_rootfind(self.equation, 0.2)
        print  str(self)

if __name__=="__main__": # run tests
    gasProperties(P=1013.25, V=1., T=273.15).solve()
    print "--- test2---"
    gasProperties( V=1,n = 0.446175, T=273.15).solve()

The benefit of root finding is that even if your formula wouldn't be so easy it would still work, so any number of formulas could be done with no more code than writing formulation. This is generally a very useful skill to have. SYMPY is good but symbolic math is not always easily solvable

The root solver is easily extendable to vector and multi equation cases, even matrix solving. The ready made scipy functions built for optimization allready do this by default.

Here is some more resources:

• Some numerical methods in python

* most were introduced at minimum to Newton–Raphson method