diff options
Diffstat (limited to 'tutorials/module_3/2_roots_optimization.md')
| -rw-r--r-- | tutorials/module_3/2_roots_optimization.md | 103 |
1 files changed, 96 insertions, 7 deletions
diff --git a/tutorials/module_3/2_roots_optimization.md b/tutorials/module_3/2_roots_optimization.md index 73bf56c..49c6273 100644 --- a/tutorials/module_3/2_roots_optimization.md +++ b/tutorials/module_3/2_roots_optimization.md @@ -176,7 +176,7 @@ Using Recursion def newtonraphson(f, df, x0, tol): # Output of this function is the estimated root of f # using the Newton-Raphson method - + # recursive method if abs(f(x0)) < tol: return x0 else: @@ -189,12 +189,47 @@ f_prime = labmda x: 2*x Using Iteration ```python -def newtonraphson(f, df, x0, tol): - for i in -``` +import numpy as np + +def newton_raphson(f, f_prime, x0, tol=1e-8): + """ + Newton-Raphson root finding method. + Parameters: + f : function, the equation f(x) = 0 + f_prime : function, derivative of f(x) + x0 : float, initial guess + tol : float, stopping tolerance + Returns: + x : float, approximate root + history : list, values of x at each iteration + """ + + x = x0 + roots = [x] + + while True: + x_new = x - f(x) / f_prime(x) + roots.append(x_new) + if abs(x_new - x) < tol: + sol = x_new + total_iter = len(roots) + return sol, total_iter, roots + + x = x_new + +# Example: solve x^2 - 2 = 0 +f = lambda x: x**2 - 2 +f_prime = lambda x: 2*x +sol, total_iter, root = newton_raphson(f, f_prime, x0=1.4) +actual=np.sqrt(2) +print("Solution=", sol) +print("No. iterations:", total_iter) +print("Iteration history:", root) +print("Actual sqrt(2) =", np.sqrt(2)) +``` ## Modified Secant A possible issue with the Newton-Raphson method is that we are required to know the derivative of the function we are finding the root for. Sometimes this may be extremely difficult or impossible to find analytically. Thus we can use the modified secant method. A numerical variant of the Newton-Raphson. @@ -226,12 +261,66 @@ Open methods may also have it's problems. Let's consider the following function: {Divergence Demo} -# Pipe Friction Example +## Problem 2: Pipe friction +The flow of fluids through pipes is a common problem in engineering. In mechanical and aerospace fields, it often arises in applications such as liquid and gas transport for cooling systems. The resistance to this flow is described by a dimensionless parameter known as the friction factor. For turbulent flow, the Colebrook equation provides a method to determine this factor. +$$ +0=\frac{1}{\sqrt{f}}+2log(\frac{\epsilon}{3.7D}+\frac{2.51}{Re\sqrt{f}}) \tag{2.2} +$$ +where +$$ +Re=\frac{\rho VD}{\mu} +$$ +Suppose you're working on a carburetor which mixes air and fuel before. You are tasked with finding the friction factor of the pipes before entering the carburetor. A lab technician measure and provides the following data of the tube to you: +$D=0.005 m$ +$\epsilon=0.0015mm$ +$V=40 m/s$ +We know we are working with air and therefore: +$\rho= 1.23 \frac{kg}{m^3}$ +$\mu=1.79*10^{-5}Ns/m^3$ -Numerical methods for Engineers 7th Edition Case study 8.4. +If we define a function $g(f)$ equal to equation above. The root of this function will give us the friction factor of this tube. +a) Plot the function of $g(f)$ to select a good initial guess. +b) Determine the root by using the Newton-Raphson method. +```python +import numpy as np +import matplotlib.pyplot as plt +# Given parameters +rho = 1.23 # kg/m^3 +mu = 1.79e-5 # Pa·s (N·s/m^2) +D = 0.005 # m +V = 40 # m/s +e = 0.0015 / 1000 # convert mm to m + +# Reynolds number +Re = rho * V * D / mu + +# Define g(f) function (Colebrook residual) +g = lambda f: 1/np.sqrt(f) + 2.0 * np.log10((e/(3.7*D)) + (2.51/(Re*np.sqrt(f)))) + +# Range of f values to check +f_values = np.linspace(0.008, 0.08, 500) +g_values = g(f_values) + +# Plot +plt.figure(figsize=(8,6)) +plt.plot(f_values, g_values) +plt.axhline(0, color='red', linestyle='--') # root line +plt.grid(True) +plt.xlabel('f') +plt.ylabel('g(f)') +plt.title('Colebrook Residual Function g(f)') +plt.show() +``` + +```python +g_prime = lambda x: 2*x +x= = 1.5 + +sol, total_iter, root = newton_raphson(g, g_prime, x0) +``` # Systems of Non-linear equations @@ -247,7 +336,7 @@ $$ $$ y+3xy^2=57 $$ -these are two nonlinear equations with two unknowns. To solve this system we can ^ +these are two nonlinear equations with two unknowns. To solve this system we can This roots of single equations and solved system of equations, but what if the system of equations are dependent on non-linear. |