import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import fsolve
To solve a system $Ax=b$ for x, we do the following:
Example
Solve $$\left[ \begin{array}[ccccc] --2 & 1 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 \\ 0 & 1 & -2& 1 & 0 \\ 0 & 0 & 1 & -2& 1 \\ 0 & 0 & 0 & 1 & -2 \end{array} \right]\cdot \left[ \begin{array}[c] \\x_0\\x_1\\x_2\\x_2\\x_4\end{array}\right] = \left[\begin{array}[c] \\1\\1\\1\\1\\1\end{array}\right]$$
A = np.array([[-2,1,0,0,0],[1,-2,1,0,0],[0,1,-2,1,0],[0,0,1,-2,1],[0,0,0,1,-2]])
b = np.ones(5)
x = np.linalg.solve(A,b)
x
array([-2.5, -4. , -4.5, -4. , -2.5])
Here's a link to more linear algebra operations.
fsolve
, based on what you need to solve equations:from scipy.optimize import fsolve
def f(x):
return x**2 - 5
x_guess = 1
x = fsolve(f, x_guess)
x
array([2.23606798])
Solving a nonlinear equation involves first putting it in the form $f(x)=0$. We know the function $f$ and we want to find the value of $x$ that gives $f(x)=0$.
To solve this, do the following:
from scipy.optimize import fsolve
fsolve
function: x = fsolve(f, x0)
Solve for both roots of $$f(x) = x^2-5.$$
# the problem is solved above. Here's the solution using an anonymous function:
fsolve(lambda x: x**2 - 5, x_guess)
array([2.23606798])
fsolve
function can handle this?fsolve
how would you account for all these possibilities?$$g(y,z) = \sin(y)/z$$
which becomes: $$\left[\begin{matrix} f_0(x_0,x_1) \\ f_1(x_0,x_1)\end{matrix}\right] = \left[\begin{matrix}x_0+2x_1 \\ \sin(x_0)/x_1\end{matrix}\right] = \left[\begin{matrix}0 \\ 0\end{matrix}\right],$$ or $$\vec{f}(\vec{x}) = \left[\begin{matrix}x_0+2x_1 \\ \sin(x_0)/x_1\end{matrix}\right] = \vec{0},$$ that is, $$\vec{f}(\vec{x})=\vec{0}.$$
Note Here is a link to more information about the fsolve function with extra arguments that can give more control over the solution process.
$$g(y,z) = \sin(y)/z.$$
def hg(yz):
y = yz[0]
z = yz[1]
h = y + 2*z
g = np.sin(y)/z
return np.array([h,g])
#---------------
yz_guess = np.array([1,2])
yz = fsolve(hg, yz_guess)
print("yz =", yz, "hg =", hg(yz))
yz = [ 9.42477796 -4.71238898] hg = [0.00000000e+00 4.82244674e-15]
Solve the following system of three equations in three unknowns:
$$x^2 + y^2 = 1,$$$$xy + yz = -1.1,$$$$y^2 + z^2 = 2.$$A reasonable guess for all variables is $x=y=z=2.$
def f(xyz):
x = xyz[0]
y = xyz[1]
z = xyz[2]
f0 = x**2 + y**2 - 1
f1 = x*y + y*z + 1.1
f2 = y**2 + z**2 - 2
return np.array([f0, f1, f2])
#------------
xyz_guess = np.array([2.0, 2.0, 2.0])
xyz = fsolve(f, xyz_guess)
print(xyz)
[ 0.10056089 -0.99493091 1.00504353]
The Colbrook equation is used to compute pressure drops across pipes, and therefore to compute pumping requirements and costs.
The Colebrook equation is given by $$\frac{1}{\sqrt{f}} = -2\log_{10}\left(\frac{\epsilon/D}{3.7}+\frac{2.51}{Re\sqrt{f}}\right),$$ or, $$F(f, Re, \epsilon/D) = \frac{1}{\sqrt{f}} + 2\log_{10}\left(\frac{\epsilon/D}{3.7}+\frac{2.51}{Re\sqrt{f}}\right)=0.$$
Approach
Re = np.logspace(3,8,1000)
epsD = 0.001
fguess = 0.01
def F(f, Re, epsD):
return 1/np.sqrt(f) + 2*np.log10(epsD/3.7 + 2.51/Re/np.sqrt(f))
# first, just get it to work with just one Re
f = fsolve(F, fguess, args=(5000, epsD))
print("for Re=5000, f =", f)
# now, make an array of f and fill in each value in a loop:
f = np.empty(len(Re))
for i in range(len(f)):
f[i] = fsolve(F, fguess, args=(Re[i], epsD)) # note the Re[i]
#--------- plot results
plt.rc('font', size=14)
plt.semilogx(Re, f)
plt.xlabel('Re')
plt.ylabel('f');
for Re=5000, f = [0.03849536]