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)
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 = fsolve(f, 2)
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.$$
def f(x):
return x**2 - 5
x0 = -2
x = fsolve(f, x0)
print(x, f(x))
[-2.23606798] [-1.77635684e-15]
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.$$
import numpy as np
def hg(yz):
y = yz[0]
z = yz[1]
h = y + 2*z
g = np.sin(y)/z
return np.array([h,g])
yz0 = np.array([1,2])
yz = fsolve(hg, yz0)
yz
array([ 9.42477796, -4.71238898])
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])
xyz0 = np.array([2.0,2,2])
xyz = fsolve(f, xyz0)
x = xyz[0]
y = xyz[1]
z = xyz[2]
print(x,y,z, f(xyz))
0.10056088509586326 -0.9949309063391776 1.0050435272221339 [ 1.60760294e-13 -1.38555833e-13 5.32907052e-15]
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))
fsolve(F, fguess, args=())