Solve Equations in Python

The following tutorials are an introduction to solving linear and nonlinear equations with Python. The solution to linear equations is through matrix operations while sets of nonlinear equations require a solver to numerically find a solution.

Source Code for Linear Solutions

import numpy as np

A = np.array([ [3,-9], [2,4] ])
b = np.array([-42,2])
z = np.linalg.solve(A,b)
print(z)

M = np.array([ [1,-2,-1], [2,2,-1], [-1,-1,2] ])
c = np.array([6,1,1])
y = np.linalg.solve(M,c)
print(y)

Source Code for Nonlinear Solution (fsolve)

import numpy as np
from scipy.optimize import fsolve

def myFunction(z):
x = z[0]
y = z[1]
w = z[2]

F = np.empty((3))
F[0] = x**2+y**2-20
F[1] = y - x**2
F[2] = w + 5 - x*y
return F

zGuess = np.array([1,1,1])
z = fsolve(myFunction,zGuess)
print(z)

Source Code for Nonlinear Solution (gekko)

from gekko import GEKKO
m = GEKKO()
x,y,w = [m.Var(1) for i in range(3)]
m.Equations([x**2+y**2==20,\
y-x**2==0,\
w+5-x*y==0])
m.solve(disp=False)
print(x.value,y.value,w.value)

Symbolic Solution with Sympy

Sympy is a package for symbolic solutions in Python that can be used to solve systems of equations.

$$2x^2+y+z=1$$ $$x+2y+z=c_1$$ $$-2x+y=-z$$

import sympy as sym
sym.init_printing()
x,y,z = sym.symbols('x,y,z')
c1 = sym.Symbol('c1')
f = sym.Eq(2*x**2+y+z,1)
g = sym.Eq(x+2*y+z,c1)
h = sym.Eq(-2*x+y,-z)

sym.solve([f,g,h],(x,y,z))

When solved in an IPython environment such as a Jupyter notebook, the result is displayed as:

$$x=-\frac{1}{2}+\frac{\sqrt{3}}{2}$$ $$y=c_1 - \frac{3 \sqrt{3}}{2} +\frac{3}{2}$$ $$z=-c_1 - \frac{5}{2} +\frac{5 \sqrt{3}}{2}$$

and a second solution:

$$x=-\frac{1}{2}-\frac{\sqrt{3}}{2}$$ $$y=c_1 + \frac{3 \sqrt{3}}{2} +\frac{3}{2}$$ $$z=-c_1 - \frac{5}{2} -\frac{5 \sqrt{3}}{2}$$

The same approach applies to linear or nonlinear equations.

Equations with Multiple Solutions (Roots)

Nonlinear equations may have multiple solutions. An example from engineering is the Redlich/Kwong Equation of State (EOS) that has 3 roots (solutions) that are the liquid volume, the vapor volume, and an extra root that isn't physically meaningful.

$$P = \frac{R_g \, T}{V-b} - \frac{a}{T^{1/2} \, V \, \left(V+b\right)}$$

where T is the temperature, V is the molar volume, R_g is the universal gas constant, and a and b are compound-specific constants. The solution for the molar volume of ethane for each phase at T = 77°C and P = 1 bar is shown below with Scipy fsolve and Gekko. For ethane, a = 2.877e8 cm^6 K^0.5 bar / mol^2 and b = 60.211 cm^3 / mol. The ideal gas law is a guess of the vapor volume and 1.1*b is a guess of the liquid volume.

import numpy as np
from scipy.optimize import fsolve

# constants
TC = 77 # degC
P = 1.0 # bar
a = 2.877e8 # cm^6 bar K^0.5 / mol^2
b = 60.211  # cm^3 / mol
Rg = 83.144598 # cm^3 bar / K-mol
# derived quantities
TK = TC+273.15 # K
Vg = Rg*TK/P # ideal gas guess

# Scipy fsolve solution
def f(V):
return P-Rg*TK/(V-b) + a/(np.sqrt(TK)*V*(V+b))
V = fsolve(f,Vg)
print(f'Scipy Vapor Solution: {V[0]:0.1f} cm^3/mol')
V = fsolve(f,b*1.1)
print(f'Scipy Liquid Solution: {V[0]:0.1f} cm^3/mol')

# Gekko solutions
from gekko import GEKKO
m = GEKKO(remote=False)
V = m.Var(value=Vg,lb=1e4)
m.Equation(P == Rg*TK/(V-b) - a/(np.sqrt(TK)*V*(V+b)))
m.options.SOLVER=1; m.solve(disp=False)
print(f'Gekko Vapor Solution: {V.value[0]:0.1f} cm^3/mol')
# update bounds and initial guess
V.lower=b*0.5; V.upper = b*2; V.value=b*1.1
m.options.SOLVER=1; m.solve(disp=False)
print(f'Gekko Liquid Solution: {V.value[0]:0.1f} cm^3/mol')