summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorChristian Kolset <christian.kolset@gmail.com>2025-09-29 08:24:55 -0600
committerChristian Kolset <christian.kolset@gmail.com>2025-09-29 08:24:55 -0600
commit0678159d2801c99437e728823dfccccb9d608174 (patch)
treebe7fd43a92eee9331413950e73be553e05ffc4bd
parent695a279f26bdf526a2ba90f0a65501067051dbe3 (diff)
Added work and newton-raphson example
-rw-r--r--README.md26
-rw-r--r--tutorials/module_3/2_roots_optimization.md103
-rw-r--r--tutorials/module_3/4_numerical_integration.md50
3 files changed, 144 insertions, 35 deletions
diff --git a/README.md b/README.md
index d85070c..6a610a8 100644
--- a/README.md
+++ b/README.md
@@ -50,8 +50,30 @@ See the tutorials under the /tutorials/ directory or use the links below.
- Verification and Validation
- [Error](tutorials/module_2/error.md)
#### Module 3: Applications of Computational Mathematics in ME
-- Systems of Equations and LU Decomposition
-- Nonlinear Equation Solver
+- Linear Equations
+ - [Linear Equations](1_linear_equations.md)
+ - Linear Algebra
+ - LU Decomposition
+- Non-Linear Equations
+ - Root fining methods
+ - Bracketing Methods
+ - Open Methods
+ - Systems of Non-linear equations
+ - Systems of Non-linear equations II
+- Numerical differentiation and integration application
+ - Differentiation
+ - Integration
+ - ODE's
+ - Explicit Methods
+ - Implicit Methods
+ - Systems of ODE's
+ - PDE's
+
+
+
+
+
+
- [[numerical_differentiation]]
- [[numerical_integration]]
- [Ordinary Differential Equations](tutorials/module_3/ode)
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$).