diff options
Diffstat (limited to 'tutorials/module_3/2_roots_optimization.md')
| -rw-r--r-- | tutorials/module_3/2_roots_optimization.md | 500 |
1 files changed, 0 insertions, 500 deletions
diff --git a/tutorials/module_3/2_roots_optimization.md b/tutorials/module_3/2_roots_optimization.md deleted file mode 100644 index 49c6273..0000000 --- a/tutorials/module_3/2_roots_optimization.md +++ /dev/null @@ -1,500 +0,0 @@ -# Root Finding Methods - -By now you've learned the quadratic formula to find the roots of a second degree polynomial. -$$ -x=\frac{-b \pm \sqrt{b^2-4ac}}{2a} -$$ -This works nicely if we're given the function $f(x)$. However, how would you solve for $x$ if you were given experimental data? In this section we will cover bracketing and open methods to find roots of a function. - -Let us take a step back and think *algorithmically*. To find a root numerically, we can implement a search algorithm. Let's consider a small phone book - a database of names and numbers sorted alphabetically by last name. Our goal is to find the phone number of James Watt (our root). We could make an algorithm that starts at the the first page, checks to see if the current name equals our target (Watt), if it doesn't, check the next entry. This is method is called *incremental search*. As you may see, this can take very long if the database is large. - -> [!NOTE] Example Phonebook Database Entries -> **A** – Archimedes -> **B** – Bernoulli -> **C** – Curie -> **D** – Darwin -> **E** – Einstein -> **F** – Feynman -> **G** – Galileo -> **H** – Hubble -> **I** – Ibn Sina -> **J** – Joule -> **K** – Kelvin -> **L** – Leonardo -> **M** – Maxwell -> **N** – Newton -> **O** – Oppenheimer -> **P** – Pascal -> **Q** – Quetelet -> **R** – Röntgen -> **S** – Schrödinger -> **U** – Urey -> **V** – Volta -> **W** – Watt -> **X** – Xenophon -> **Y** – Yukawa -> **Z** – Zuse - -Another approach may be to index the book and jump straight to the 'W' names however this requires two steps. - -When - -# Bracketing Method -## Incremental Search -The incremental search method is the simplest of the bracketing methods. The basic idea is to start from an initial point and move forward in small increments along the function. At each step, the sign of the function value is checked. If a change of sign occurs between two consecutive points, then a root must lie between them (by the **Intermediate Value Theorem**). This gives us a bracket — an interval that contains the root — which can then be refined using more powerful methods such as bisection or false position. - -Although incremental search is not efficient for accurately finding roots, it is useful for _detecting the presence of roots_ and identifying approximate intervals where they exist. Once such an interval is located, other methods can be applied to converge quickly to the actual root. In this sense, incremental search often serves as a **pre-processing step** for bracketing or open methods, especially when the root’s location is not known in advance. -## Bisection -This method start by sectioning the database into two halves. We can do this with a phone book by simply opening it in the middle of the book. In our data set, we have find the middle entry to be 'M - Maxwell'. Based on the premises that all names are sorted alphabetically we can conclude that we will find our target in the second half, essentially eliminating the first half. Doing this a second time gives us our new midpoint 'S – Schrödinger' and a third time get us to our target 'W - Watt'. In only 3 iterations we have reached out target goal. This is a huge improvement from the Incremental search (22 iterations) versus our new 3 iterations. - -This exact same concept can be applied to finding roots of functions. Instead of searching for strings we are searching for floats. For computers this is a lot easier to do. Take a moment to think about how you would write psuedo-code to program this. - -> The **Intermediate Value Theorem** says that if f(x) is a continuous function between a and b, and sign(f(a))≠sign(f(b)), then there must be a c, such that a<c<b and f(c)=0. - -Let's consider a continuous function $f(x)$ with an unknown root $x_r$ . Using the intermediate value theorem we can bracket a root with the points $x_{lower}$ and $x_{upper}$ such that $f(x_{lower})$ and $f(x_{upper})$ have opposite signs. This idea is visualized in the graph below. -<img - style="display: block; - margin-left: auto; - margin-right: auto; - width: 50%;" - src="https://pythonnumericalmethods.studentorg.berkeley.edu/_images/19.03.01-Intermediate-value-theorem.png" - alt="Intermediate Value Theorem"> -Once we bisect the interval and found we set the new predicted root to be in the middle. We can then compare the two sections and see if there is a sign change between the bounds. Once the section with the sign change has been identified, we can repeat this process until we near the root. -<img - style="display: block; - margin-left: auto; - margin-right: auto; - width: 50%;" - src="bisection.png" - alt="Bisection"> -As you the figure shows, the predicted root $x_r$ get's closer to the actual root each iteration. In theory this is an infinite process that can keep on going. In practice, computer precision may cause error in the result. A work-around to these problems is setting a tolerance for the accuracy. As engineers it is our duty to determine what the allowable deviation is. - -So let's take a look at how we can write this in python. -```python -import numpy as np - -def my_bisection(f, x_l, x_u, tol): - # approximates a root, R, of f bounded - # by a and b to within tolerance - # | f(m) | < tol with m the midpoint - # between a and b Recursive implementation - - # check if a and b bound a root - if np.sign(f(x_l)) == np.sign(f(x_u)): - raise Exception( - "The scalars x_l and x_u do not bound a root") - - # get midpoint - x_r = (x_l + x_u)/2 - - if np.abs(f(x_r)) < tol: - # stopping condition, report m as root - return x_r - elif np.sign(f(x_l)) == np.sign(f(x_r)): - # case where m is an improvement on a. - # Make recursive call with a = m - return my_bisection(f, x_r, x_u, tol) - elif np.sign(f(x_l)) == np.sign(f(x_r)): - # case where m is an improvement on b. - # Make recursive call with b = m - return my_bisection(f, x_l, x_r, tol) -``` - -### Problem 1 -Call the *my_bisection* function to solve the equation $f(x) = x^2 -2$ -```python -f = lambda x: x**2 - 2 - -r1 = my_bisection(f, 0, 2, 0.1) -print("r1 =", r1) -r01 = my_bisection(f, 0, 2, 0.01) -print("r01 =", r01) - -print("f(r1) =", f(r1)) -print("f(r01) =", f(r01)) -``` - - -### Problem 2 -Write an algorithm that uses a combination of the bisection and incremental search to find multiple roots in the first 8 minutes of the data set. -```python -import pandas as pd -import matplotlib.pyplot as plt - -# Read the CSV file into a DataFrame called data -data = pd.read_csv('your_file.csv') - -# Plot all numeric columns -data.plot() - -# Add labels and show plot -plt.xlabel("Time (min)") -plt.ylabel("Strain (kPa)") -plt.title("Data with multiple roots") -plt.legend() -plt.show() -``` - -## False Position - -# Open Methods -So far we looked at methods that require us to bracket a root before finding it. However, let's think outside of the box and ask ourselves if we can find a root with only 1 guess. Can we improve the our guess by if we have some knowledge of the function itself? The answer - derivatives. - -## Newton-Raphson -Possibly the most used root finding method implemented in algorithms is the Newton-Raphson method. By using the initial guess we can draw a tangent line at this point on the curve. Where the tangent intercepts the x-axis is usually an improved estimate of the root. It can be more intuitive to see this represented graphically. -<img - style="display: block; - margin-left: auto; - margin-right: auto; - width: 50%;" - src="https://pythonnumericalmethods.studentorg.berkeley.edu/_images/19.04.01-Newton-step.png" - alt="Newton-Raphson illustration"> -Given the function $f(x)$ we make our first guess at $x_0$. Using our knowledge of the function at that point we can calculate the derivative to obtain the tangent. Using the equation of the tangent we can re-arrange it to solve for $x$ when $f(x)=0$ giving us our first improved guess $x_1$. - -As seen in the figure above. The derivative of the function is equal to the slope. -$$ -f'(x_0) = \frac{f(x_0)-0}{x_0-x_1} -$$ -Re-arranging for our new guess we get: -$$ -x_{1} = x_0 - \frac{f(x_0)}{f'(x_0)} -$$ -Since $x_0$ is our current guess and $x_0$ is our next guess, we can write these symbolically as $x_i$ and $x_{i+1}$ respectively. This gives us the *Newton-Raphson formula*. -$$ -\boxed{x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}} \tag{3.1} -$$ -### Problem 1 -From experimental data we extrapolated a function f. Write a python function called *newtonraphson* which as the following input parameters - - `f` - function - - `df` - first derivative of the function - - `x0` - initial guess - - `tol` - the tolerance of error - That finds the root of the function $f(x)=e^{\left(x-6.471\right)}\cdot1.257-0.2$ and using an initial guess of 8. - -Using Recursion -```python -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: - return newtonraphson(f,df,x0) - -f = lambda x: x**2-1.5 -f_prime = labmda x: 2*x - -``` - -Using Iteration -```python -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. - -This method uses the finite derivative to guess where the root lies. We can then iterate this guessing process until we're within tolerance. - -$$ -f'(x_i) \simeq \frac{f(x_{i-1})-f(x_i)}{x_{i-1}x_i} -$$ -Then substituting this in equation (3.1) we get the equation for the secant method: - -$$ -x_{i+1} = x_i - \frac{f(x_i) (x_{i-1}x_i)}{f(x_{i-1})-f(x_i)} \tag {3.2} -$$ - -We can then *modify* the equation to instead of choosing two arbitrary numbers we use a small value $\delta$ in the independent variable to get a better reading. -$$ -x_{x+1} = x_i - \frac{f(x_i)(x_{i-1}-x_i)}{f(x_i+\delta x)-f(x_i)} -$$ -$$ - -\boxed{x_{x+1} = x_i - \frac{f(x_i)\delta x_i}{f(x_i+\delta x)-f(x_i)}} \tag{3.3} -$$ -Note that if $\delta$ is too small the method can fail by round-off error due to the subractive cancellation in the denominator in eqn. 3.3. If $\delta$ is too large, the method becomes inefficient or even divergent. - -## Issues with open methods - -Open methods may also have it's problems. Let's consider the following function: $f(x)=3.2*tan^{-1}(x-4.3)$ if we were to apply the Newton-Raphson to find the root of this function we can see that the results diverge. - -{Divergence Demo} - - -## 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$ - -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 -So far we've solved system of linear equations where the equations come in the form: -$$ -f(x) = a_1x_1 + a_2x_2 + \ldots + a_nx_n-b = 0 -$$ -where the $a$'s and $b$ are constants. Equations that do not fit this format are called *non linear* equations. For example: - -$$ -x^2 + xy = 10 -$$ -$$ -y+3xy^2=57 -$$ -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. - -$$ - \begin{align*} -f_{1}(x_{1},x_{2},\ldots ,x_{n}) &= 0 \\ -f_{2}(x_{1},x_{2},\ldots ,x_{n}) &= 0 \\ -&\ \vdots \\ -f_{n}(x_{1},x_{2},\ldots ,x_{n}) &= 0 -\end{align*} -$$ - - - -We've applied the Newton-Raphson to - - - -## Problem 1: Newton-Raphson for a Nonlinear system - -Use the multiple-equation Newton-Raphson method to determine the roots of the following system of non-linear equations -$$ -\begin{cases} -u(x,y) = x^2 + xy - 10 = 0 \\ -v(x,y) = y + 3xy^2 - 57 = 0 -\end{cases} -$$ -after - a) 1 iteration - b) Write a python function that iterates until the tolerance is within 0.002 - -Note the solution to $x$ =2 and $y$ = 3. Initiate the computation with guesses of $x$ = 1.5 and $y$ = 3.5. - -Solution: -a) -```python -u = lambda x, y: x**2 + x*y - 10 -v = lambda x, y: y + 3*x*y**2 - 57 - -# Initial guesses -x_0 = 1.5 -y_0 = 3.5 - -# Evaluate partial derivatives at initial guesses -dudx = 2*x_0+y_0 -dudy = x_0 - -dvdx = 3*y_0**2 -dvdy = 1 + 6*x_0*y_0 - -# Find determinant of the Jacobian -det = dudx * dvdy - dudy * dvdx - -# Values of functions valuated at initial guess -u_0 = u(x_0,y_0) -v_0 = v(x_0,y_0) - -# Substitute into newton-raphson equation -x_1 = x_0 - (u_0*dvdy-v_0*dudy) / det -y_1 = y_0 - (v_0*dudx-u_0*dvdx) / det - -print(x_1) -print(y_1) -``` - -b) -```python -import numpy as np - -def newton_raphson_system(f, df, g0, tol=1e-8, max_iter=50, verbose=False): - """ - Newton-Raphson solver for a 2x2 system of nonlinear equations. - - Parameters - ---------- - f : array of callables - f[i](x, y) should evaluate the i-th equation. - df : 2D array of callables - df[i][j](x, y) is the partial derivative of f[i] wrt variable j. - g0 : array-like - Initial guess [x0, y0]. - tol : float - Convergence tolerance. - max_iter : int - Maximum iterations. - verbose : bool - If True, prints iteration details. - - Returns - ------- - (x, y) : solution vector - iters : number of iterations - """ - x, y = g0 - - for k in range(1, max_iter + 1): - # Evaluate system of equations - F1 = f[0](x, y) - F2 = f[1](x, y) - - # Evaluate Jacobian entries - J11 = df[0][0](x, y) - J12 = df[0][1](x, y) - J21 = df[1][0](x, y) - J22 = df[1][1](x, y) - - # Determinant - det = J11 * J22 - J12 * J21 - if det == 0: - raise ZeroDivisionError("Jacobian is singular") - - # Solve for updates using 2x2 inverse - dx = ( F1*J22 - F2*J12) / det - dy = ( F2*J11 - F1*J21) / det - - # Update variables - x_new = x - dx - y_new = y - dy - - if verbose: - print(f"iter {k}: x={x_new:.6f}, y={y_new:.6f}, " - f"|dx|={abs(dx):.2e}, |dy|={abs(dy):.2e}") - - # Convergence check - if max(abs(dx), abs(dy)) < tol: - return np.array([x_new, y_new]), k - - x, y = x_new, y_new - - raise ValueError("Newton-Raphson did not converge") - - -# ------------------- -# Example usage -# ------------------- - -# Define system: u(x,y), v(x,y) -u = lambda x, y: x**2 + x*y - 10 -v = lambda x, y: y + 3*x*y**2 - 57 - -# Partial derivatives -dudx = lambda x, y: 2*x + y -dudy = lambda x, y: x -dvdx = lambda x, y: 3*y**2 -dvdy = lambda x, y: 1 + 6*x*y - -f = np.array([u, v]) -df = np.array([[dudx, dudy], - [dvdx, dvdy]]) - -g0 = np.array([1.5, 3.5]) - -sol, iters = newton_raphson_system(f, df, g0, tol=1e-10, verbose=True) -print("Solution:", sol, "in", iters, "iterations") - -``` - - - - |