Added work and newton-raphson example

2 files changed, 120 insertions, 33 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.
diff --git a/tutorials/module_3/4_numerical_integration.md b/tutorials/module_3/4_numerical_integration.md
index 0c2e755..4448985 100644
--- a/tutorials/module_3/4_numerical_integration.md
+++ b/tutorials/module_3/4_numerical_integration.md
@@ -145,7 +145,7 @@ while (iter_ < iter_max) and (err > tol):
 ```
 
 ### Gaussian Quadrature
-Gaussian quadrature chooses evaluation points and weights optimally (based on Legendre polynomials) to maximize accuracy. With $n$ evaluation points, Gaussian quadrature is exact for polynomials of degree up to $2n-1$. This makes it extremely efficient for smooth integrands.
+Gaussian quadrature chooses evaluation points and weights optimally to maximize accuracy. With $n$ evaluation points, Gaussian quadrature is exact for polynomials of degree up to $2n-1$. This makes it extremely efficient for smooth integrands.
 
 ```python
 import numpy as np
@@ -205,29 +205,40 @@ print(f"Error (3-point): {err3:.2e}")
 
 
 
-
-
-
-
----
-
-
 ## More integral
 
 ### Newton-Cotes Algorithms for Equations
 ### Adaptive Quadrature
 
 
-
-
 # Problems 
 
 ## Numerical Integration to Compute Work
-
+In physics we've learned that work is computed
+$$
+Work = force * distance
+$$
+This can be written in int's integral form:
+$$
+W = \int{F(x)dx}
+$$
+If F(x) is easy to integrate, we could solve this problem analytically. However, a realistic problem the force may not be available to you as a function, but rather, tabulated data. Suppose some measurements were take of when a weighted box was pulled with a wire. If we data of the force on the wire and the angle of the wire from the horizontal plane.
+
+| x (ft) | F(x) (lb) | θ (rad) | F(x) cos θ |
+| ------ | --------- | ------- | ---------- |
+| 0      | 0.0       | 0.50    | 0.0000     |
+| 5      | 9.0       | 1.40    | 1.5297     |
+| 10     | 13.0      | 0.75    | 9.5120     |
+| 15     | 14.0      | 0.90    | 8.7025     |
+| 20     | 10.5      | 1.30    | 2.8087     |
+| 25     | 12.0      | 1.48    | 1.0881     |
+| 30     | 5.0       | 1.50    | 0.3537     |
+|        |           |         |            |
+Use the trapezoidal rule to compute the work done on the box.
 
 ## Implementing the Composite Trapezoidal Rule
 
-**Objective:** Implement a Python function to approximate integrals using the trapezoidal rule.
+Implement a Python function to approximate integrals using the trapezoidal rule.
 
 ```python
 import numpy as np
@@ -245,17 +256,4 @@ for n in [4, 8, 16, 32]:
     print(n, trapz(f1, 0, np.pi, n))
 ```
 
-Students should compare results for increasing $n$ and observe how the error decreases with $O(h^2$).
-
----
-## Gaussian Quadrature
-
-Write a Python function for two-point and three-point Gauss–Legendre quadrature over an arbitrary interval $[a,b]$. Verify exactness for polynomials up to the appropriate degree and compare performance against the trapezoidal rule on oscillatory test functions.
-
----
-
-## Simpson’s 1/3 Rule
-
-Code the composite Simpson’s 1/3 rule. Test its accuracy on smooth functions and compare its performance to the trapezoidal rule and Gaussian quadrature. Document error trends and discuss cases where Simpson’s method is preferable.
-
----
+Compare the results for increasing $n$ and observe how the error decreases with $O(h^2$).