$$x_{k+1} = g(x_k)$$ * That is: guess $x_0$. * Evaluate $g(x_0)$. * The result is $x_1$. * Repeat.
Geometrically $x_{k+1} = g(x_k)$ finds the intersection of the $y=g(x)$ and $y=x$ lines.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x = np.linspace(-2,3,100)
y1 = x**2 - 2.0
y2 = 1.0+2./x
y3 = (x+2.)**0.5
y4 = (x**2 + 2)/(2*x -1)
plt.rc("font", size=16)
plt.plot(x,y1,label=r'$x^2-2$')
plt.plot(x,y2,label=r'$1+2/x$')
plt.plot(x,y3,label=r'$(x+2)^{1/2}$')
plt.plot(x,y4,label=r'$(x^2+2)/(2x-1)$')
plt.plot(x,x,'k--',label=r'$x$')
plt.legend(loc='upper left', frameon=False, fontsize=12)
plt.ylim([-1,8]); # plt.ylim([-3,8])
plt.xlim([1,3]); # plt.xlim([-2,3])
plt.xlabel('x')
plt.ylabel('g(x)');
import numpy as np
def FP(g, x, tol=1E-5, maxit=1000):
for niters in range(1,maxit+1):
xnew = g(x)
err = np.linalg.norm(xnew-x)/np.linalg.norm(x)
if err <= tol:
return xnew, niters
x = xnew
print(f"Warning no converged in {maxit} iterations")
return xnew, niters
def FP(g, x, tol=1E-5, maxit=1000):
for niter in range(1,maxit+1):
xnew = g(x)
if np.linalg.norm(xnew-x)/np.linalg.norm(x) <=tol:
return xnew, niter
x = xnew
print('warning, reached maxit =', maxit)
return x, niter
#-----------
def g(x):
return np.sqrt(x+2)
#-----------
xguess = 4.0
x, nit = FP(g, xguess)
print("x, nit = ", x, nit)
x, nit = 2.000006616131052 9
Equations:
$$ Re = \frac{\rho Dv}{\mu}.$$$$ \frac{1}{\sqrt{f}} = -\log_{10}\left(\frac{\epsilon/D}{3.7} + \frac{2.51}{Re\sqrt{f}}\right),$$$$ \frac{\Delta P}{\rho} = \frac{fLv^2}{2D} \rightarrow v = \sqrt{\frac{2D\Delta P}{\rho f L}}.$$Approach:
Note, this calls FP
from example 1
import numpy as np
#-------------------
def g_fluids(vfhat):
v = vfhat[0]
fhat = vfhat[1]
ρ = 1000 # kg/m3
μ = 1E-3 # Pa*s = kg/m*s
D = 0.1 # m
L = 100 # m
ϵ = D/100 # m
ΔP = 101325 # Pa
Re = ρ*D*v/μ
fhat = -np.log10(ϵ/D/3.7 + 2.51*fhat/Re)
v = np.sqrt(2*D*ΔP/ρ/L)*fhat
return np.array([v, fhat])
#-------------------
vfhat_g = np.array([1.0, 0.1])
vfhat, nit = FP(g_fluids, vfhat_g)
print(f"v = {vfhat[0]:.4f} (m/s)")
print(f"f = {1/vfhat[1]**2:.4f}")
print(f"# iterations = {nit}")
v = 1.1521 (m/s) f = 0.1527 # iterations = 3
x = 1
print(f"x = {x:.15f}")
for k in range(6):
x = x - (x**2 - np.pi**2)/(2*x)
print(f"x = {x:.15f}")
x = 1.000000000000000 x = 5.434802200544679 x = 3.625401431921964 x = 3.173874724746142 x = 3.141756827069927 x = 3.141592657879261 x = 3.141592653589793