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.
Solve Linear Equations with Python
Source Code for Linear Solutions
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)
Solve Nonlinear Equations with Python
Source Code for Nonlinear Solution (fsolve)
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)
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$$
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.
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')
Additional Tutorials
Linear and nonlinear equations can also be solved with Excel and MATLAB. Click on the appropriate link for additional information and source code.
The Gekko Optimization Suite with a Python interface is optimization software for mixed-integer and differential algebraic equations. It is coupled with large-scale solvers for linear, quadratic, nonlinear, and mixed integer programming (LP, QP, NLP, MILP, MINLP). Modes of operation include data reconciliation, real-time optimization, dynamic simulation, and nonlinear predictive control.
Generative AI Learning
Use these prompts to test your understanding after completing the tutorial. Equation solving is the cleanest place to practice verification: every solution can be substituted back into the equations it claims to solve.
Tip: Make print(f(x_solution)) the last line of every solve script - a residual near machine precision is the receipt. When an AI writes the solver code for you later (projects, HW), the residual check and a guess-sensitivity test are how you audit it in ten seconds. For symbolic work (class 20), the same habit holds: differentiate with sympy, then check the derivative numerically at one point.
What to Turn In
Submit a short report (PDF, 1-2 pages) that curates your results into a demonstration of what you learned. You may use Generative AI to help write the report, but you must guide it to the correct solutions, justifications, and assumptions. Answer these questions:
- Solve a linear system with np.linalg.solve and report the solution AND the residual A x - b (all elements near zero).
- Solve one nonlinear equation or system with fsolve from at least two different initial guesses. Report guesses, roots, and residuals. Did the guesses find different roots?
- For the Redlich-Kwong (or your assigned) equation: report all roots you found and state the physical meaning of each - which molar volume would you report for the vapor phase and why?
- From the audit prompt: what did the reviewer questioning expose about your solution or your understanding?
- From the quiz prompt: one question you missed and the corrected answer.
