HW08: Solve Equations#
Systems of equations are solved by finding variable values that satisfy all of the equations. With linear systems of equations, there is typically a unique solution that is easy to solve, even for large-scale systems. Nonlinear equations can be more challenging to solve and may have multiple solutions. Variable bounds and a good initial guess are sometimes required for nonlinear systems of equations. This section is to solve systems of:
Linear systems of equations
Nonlinear systems of equations
Problem #1#
Solve the following systems of linear algebraic equations \(\left(A\,x = b\right)\) using matrix calculations such as:
\(\left[ {\begin{array}{cc} 1 & 2 \\ 1.5 & -1 \\ \end{array} } \right] \left[ {\begin{array}{c} x_0 \\ x_1 \end{array} } \right] = \left[ {\begin{array}{c} 0 \\ 0.5 \end{array} } \right]\)
An array is defined with square brackets as a list or numpy array.
Vector: A vector (1D array) is:
\(b = \left[ {\begin{array}{c} 0 \\ 0.5 \end{array} } \right]\)
and is defined as a list:
b = [0,0.5]
or as a numpy array:
b = np.array([0,0.5])
Matrix: A matrix (2D array) is:
\(A = \left[ {\begin{array}{cc} 1 & 2 \\ 1.5 & -1 \\ \end{array} } \right]\)
and is defined as a list:
A = [[1.0, 2.0],[1.5,-1.0]]
or as a numpy array:
A = np.array([[1.0, 2.0],[1.5,-1.0]])
Solution Method #1: The solution to this set of equations can be found by inverting the matrix \(A\) and multiplying by \(b\) as \(x=A^{-1}\,b\).
\(\left[ {\begin{array}{c} x_0 \\ x_1 \end{array} } \right] = \left[ {\begin{array}{cc} 1 & 2 \\ 1.5 & -1 \\ \end{array} } \right]^{-1} \left[ {\begin{array}{c} 0 \\ 0.5 \end{array} } \right]\)
The inverse of a matrix can be found with the function np.linalg.inv. Two matrices are multiplied together with the dot product with function np.dot.
# practice problem
import numpy as np
A = [[1.0, 2.0],[1.5,-1.0]]
b = [0,0.5]
invA = np.linalg.inv(A)
x = np.dot(invA,b)
print(x)
[ 0.25 -0.125]
Solution Method #2: The function np.linalg.solve can also be used to find x for linear systems of equations.
x = np.linalg.solve(A,b)
print(x)
[ 0.25 -0.125]
Parse Solution: Individual values can be accessed by referencing the specific location in the array. An array index starts with 0.
print(x[0])
0.25
Last Value in Array: The last value of this array can be accessed with x[1] or with x[-1].
print(x[1])
print(x[-1])
-0.125
-0.125
Array References: The second to last value of any array can be accessed with x[-2].
print(x[-2])
0.25
Action: Find the values of x that satisfy the following system of linear equations (\(A \, x = b\)):
\(\left[ {\begin{array}{ccccc} 11 & 7 & 0 & 1 & 2 \\ 0 & 5 & 2 & 0 & 1 \\ 3 & 2 & 7 & 1 & 0 \\ 4 & 0 & 4 & 10 & 1 \\ 2 & 5& 1 & 3 & 14 \end{array} } \right] \left[ {\begin{array}{c} x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \end{array} } \right] = \left[ {\begin{array}{c} 45 \\ 30 \\ 15 \\ 20 \\ 92 \end{array} } \right]\)
A = [[11.0, 7.0, 0.0, 1.0, 2.0],\
[ 0.0, 5.0, 2.0, 0.0, 1.0],\
[ 3.0, 2.0, 7.0, 1.0, 0.0],\
[ 4.0, 0.0, 4.0, 10.0, 1.0],\
[ 2.0, 5.0, 1.0, 3.0, 14.0]]
b = [45.0, 30.0, 15.0, 20.0, 92.0]
Action: Solve the equations with a different method and compare the solution to the one above.
Action: Print the last value of x.
Problem #2#
A nonlinear equation solver such as scipy.optimize.fsolve can also solve linear equations. Returning to the introductory example:
\(\left[ {\begin{array}{cc} 1 & 2 \\ 1.5 & -1 \\ \end{array} } \right] \left[ {\begin{array}{c} x_0 \\ x_1 \end{array} } \right] = \left[ {\begin{array}{c} 0 \\ 0.5 \end{array} } \right]\)
This is solved with fsolve by placing all of the equations to one side as:
\(\left[ {\begin{array}{cc} 1 & 2 \\ 1.5 & -1 \\ \end{array} } \right] \left[ {\begin{array}{c} x_0 \\ x_1 \end{array} } \right] - \left[ {\begin{array}{c} 0 \\ 0.5 \end{array} } \right] = \left[ {\begin{array}{c} f_0 \\ f_1 \end{array} } \right]\)
The equation is solved when \(f_0\) and \(f_1\) are equal to zero. The fsolve function requires a function that returns the residuals (\(f_1\) and \(f_2\)).
#notice the use of the array or list to specify x0 and x1 (other variables may be specified like x and y in the equations as
#variables and you'd still use the array of x[:] or you could use a tuple)
def f(x):
f1 = x[0]+2*x[1]
f2 = 1.5*x[0]-x[1]-0.5
return [f1,f2]
The function solver fsolve also needs an initial guess. A few guess values are evaluated below.
# values of the function at a few points (guesses)
print(f([0.0,0.0]))
print(f([0.5,0.5]))
print(f([1.0,0.0]))
[0.0, -0.5]
[1.5, -0.25]
[1.0, 1.0]
Guess values for \(x_0\) and \(x_1\) are selected as [1,1] and fsolve finds the solution.
from scipy.optimize import fsolve
x = fsolve(f,[1,1])
print(x)
[ 0.25 -0.125]
#solving the same thing above using a tuple instead of an array
def ftuple(p):
x,y = p #unpack the tuple
f1 = x+2*y
f2 = 1.5*x-y-0.5
return (f1,f2)
xt = fsolve(ftuple,(1,1)) #pack the tuple for the guess value
print(xt)
[ 0.25 -0.125]
Action: Determine one solution for \(x\) and \(y\) to the following equations (there are 4 possible solutions).
\(2 x^2 + y^2 = 1\)
\((0.5x-0.5)^2+2(y-0.25)^2=1\)
The initial guess can be adjusted to find a different solution.
Problem #3#
The Redlich/Kwong Equation of State is
\(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.
A plot of the pressure as a function of the molar volume at a temperature of -87 C is shown below. Notice there are three roots at that temperature (three molar volumes where the curve yields 1 bar). The left most is the liquid phase and the right most is the vapor phase.
Find the molar volume of ethane for each phase (liquid and gas) that is present at \(T = -87 \, °C\) and \(P = 1 \, bar\). For ethane, \(a = 2.877×10^8 \, \frac{cm^6 bar K^{0.5}}{mol^2}\) and \(b = 60.211 \, \frac{cm^3}{mol}\).
Pay attention to units.
Use the ideal gas law to obtain the guess of the vapor volume.
Use a guess value of 1.01*b for the liquid volume.
You should obtain something around 16000 \(cm^3/mol\) for the vapor volume and 55 \(cm^3/mol\) for the liquid volume. The middle root (not required but for your information) is around 60 \(cm^3/mol\) (as shown in the above plot).
Compare your answer for the vapor volume (actual) to the ideal gas law (estimate) by calculating a percent error: (actual-estiate)/actual*100).
Also compare your answer for the liquid volume you solved for to the liquid volume of ethane at -87 °C (look up on the web or dippr). Also calculate a percent error.
import numpy as np
from scipy.optimize import fsolve
# constants
TC = -87 # degC
P = 1e5 # Pa
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
# define function for fsolve
def fRKP(V):
return # fill-in
Problem #4#
Imagine mixing liquid benzene (species 1) and toluene (species 2) together in an initially empty container. At equilibrium, some of the liquid from both species evaporates into the vapor phase and some is left liquid phase for certain temperatures and pressures. In thermodynamics, you learn that Raoult’s law may be used to describe the distribution of species in each phase. For the specific example mentioned above, Raoult’s law gives the following two expressions describing the equilibrium state
where \(y_i\) is the mole fraction of species \(i\) (either 1 or 2) in the vapor phase, \(x_i\) is the mole fraction of species \(i\) in the liquid phase, \(P_i^{sat}(T)\) is the vapor pressure of species \(i\) at the system temperature \(T\), and \(P\) is the system pressure.
The vapor pressures of benzene and toluene can be found using the Antione Equation given below:
Note that \(P_i^{sat}(T)\) is not \(P_i^{sat}\) times \(T\), but is simply showing that \(P_i^{sat}\) is a function of \(T\) in kPa. The constants, \(A\), \(B\), and \(C\), are compound specific and are found in the table below. \(T\) is in \(°C\) in the above equation.
Antoine Constants for Benzene and Toluene:
Compound |
A |
B |
C |
|---|---|---|---|
Benzene |
13.7819 |
2726.81 |
217.572 |
Toluene |
13.9320 |
3056.96 |
217.625 |
Setup Equations
Since the liquid (x’s) and vapor (y’s) mole fractions must sum to 1, we can write the following equations:
If a mixture of benzene (species 1) and toluene (species 2) has \(y_1 = 0.33\) and \(P = 120 \, kPa\), find \(x_1\) and \(T\).
For full credit, you must use an object-oriented approach to solve this problem; meaning you need to define a class and have a Psat function as part of that class. That Psat function should be written such that all the parameters except for the temperature are defined by self. That is, def Psat(self,T): and then the function should return the vapor pressure referring to the self.antione coefficients in the function.
Use a guess value of 0.5 for \(x_1\) and 100 °C for \(T\).
You should get about 0.15 for \(x_1\) and 110 °C for \(T\).
Report the solved for \(x_1\) and \(T\) values.
Credit to Dr. John Hedengren for some of these problems.