Added Systems of non-linear equations

1 files changed, 173 insertions, 0 deletions

diff --git a/tutorials/module_3/2_roots_optimization.md b/tutorials/module_3/2_roots_optimization.md
index a91a464..73bf56c 100644
--- a/tutorials/module_3/2_roots_optimization.md
+++ b/tutorials/module_3/2_roots_optimization.md
@@ -235,4 +235,177 @@ Numerical methods for Engineers 7th Edition Case study 8.4.
 
 
 # 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")
+
+```
+
+
+