diff options
Diffstat (limited to 'tutorials/module_3/1_linear_equations.md')
| -rw-r--r-- | tutorials/module_3/1_linear_equations.md | 540 |
1 files changed, 0 insertions, 540 deletions
diff --git a/tutorials/module_3/1_linear_equations.md b/tutorials/module_3/1_linear_equations.md deleted file mode 100644 index 9f50b97..0000000 --- a/tutorials/module_3/1_linear_equations.md +++ /dev/null @@ -1,540 +0,0 @@ -# Linear Equations -Let's consider an linear equation -$$ -y = m x + b -$$ -where $a$ and $b$ are two known constants we can solve for $x$ easily. - -# Linear Algebra -Although this isn't a course in linear algebra we are going to use some fundamental concepts from linear algebra to solve systems of equations. If you haven't taken linear algebra before, it is the study of linear equations. These equations can be represented in the form of matrices. Matrices are rectangular arrays of numbers arranged in rows and columns. They are widely used in both engineering and computer science. Let's say we have the following system of equation. -$$ -\begin{cases} -a_{11} x_1 + a_{12} x_2 + a_{13} x_3 = b_1 \\ -a_{21} x_1 + a_{22} x_2 + a_{23} x_3 = b_2 \\ -a_{31} x_1 + a_{32} x_2 + a_{33} x_3 = b_3 -\end{cases} -\tag{1} - -$$ -Where both the $a_{nn}$ coefficients and the $b_n$ variables are all known. We can re-write this system into a matrix form by creating an $A$ matrix of all $a_{nn}$ values and a vector in the form of: -$$ - \left[ {\begin{array}{cc} - a_{11} & a_{12} & a_{13}\\ - a_{21} & a_{22} & a_{23}\\ - a_{31} & a_{32} & a_{33}\\ - \end{array} } \right] - \left[ {\begin{array}{cc} - x_{1}\\ - x_{2}\\ - x_{3}\\ - \end{array} } \right] -= - \left[ {\begin{array}{cc} - b_{1}\\ - b_{2}\\ - b_{3}\\ - \end{array} } \right] - \tag{2} -$$ -where -$$ - \left[ {\begin{array}{cc} - a_{11} & a_{12} & a_{13}\\ - a_{21} & a_{22} & a_{23}\\ - a_{31} & a_{32} & a_{33}\\ - \end{array} } \right] - = A, - \left[ {\begin{array}{cc} - b_{1}\\ - b_{2}\\ - b_{3}\\ - \end{array} } \right] - = b -$$ -we can also write the equation in it's simplified form -$$ -Ax=b -\tag{3} -$$ - -Take a close look at equation 1 and 2 and look at their similarity and differences. Notice how the pattern of the $A$ matrix and $b$ vector and follow the same patterns as however the elements of the $x$ vector ($x_1$,$x_2$,$x_3$) are flipped vertically. This is because of how we multiple matrices - -## Matrix operations in python - -- Matrices of the same size can be added or subtracted element-wise: -- Multiply a matrix by a scalar yields a new matrix where every element is multiplied by the scalar. -- The determinant of a matrix is number of a square matrix that summarizes certain properties of the matrix. -- **Identity Matrix $I$**: Square matrix with 1’s on the diagonal and 0’s elsewhere. -- **Inverse Matrix $A^{-1}$**: Satisfies $AA^{-1} = I$. (Only square, non-singular matrices have inverses.) - -## Matrix Math in Python -Let's create two square matrices $A$ and $B$ to perform some matrix operations on: -```python -import numpy as np -np.set_printoptions(suppress=True) - -A = np.array([[1, 3], - [2, 4]]) # 2x2 -B = np.array([[5, 6], - [7, 8]]) # 2x2 -A.shape, B.shape - -b = np.array([3, 5]) -b.shape -``` - -Test and see how it works. - -```python -# Addition and Subtraction of same shapes -A + B -A - B - -# Scalar Multiplication -3*A - -5*B - -# Matrix Multiplication -mat_mult = A @ B - -# Transpose -B.T - -# Identity and Inverse matrix - -I = np.eye(2) -A_inv = np.linalg.inv(A) - -A @ A_inv, A_inv @ A - -# Determinants - -detA = np.linalg.det(A) -detB = np.linalg.det(B) -detA, detB -``` - -### Problem 1 -Re-write the following systems of equation in matrix form to normal algebraic equations. -$$ - \left[ {\begin{array}{cc} - 3 & 5 & 7\\ - 2 & 1 & 3\\ - 2 & 3 & 6\\ - \end{array} } \right] - \left[ {\begin{array}{cc} - x_{1}\\ - x_{2}\\ - x_{3}\\ - \end{array} } \right] -= - \left[ {\begin{array}{cc} - 24\\ - 52\\ - 12\\ - \end{array} } \right] -$$ -## Problem 2 -```python -import numpy as np - -# Coefficient matrix A -A = np.array([[2, 3], [4, 5]]) -# Right - hand side vector b -b = np.array([4, 6]) - -# Solve the system of equations -x = np.linalg.solve(A, b) -print(x) -``` -## Problem 3 - - - -Matrix Determinate - -Cramer's Rule - - - -# Systems of Linear Equations - -## Techniques to solve Systems of Equations -In this section our goals is to learn some algorithms that help solve system of equations in python. Let us try to solve system of equation think about how we solve a system of two linear equations by hand. -$$ -\begin{cases} 2x + y = 5 \\ 4x - y = 3 \end{cases} -$$ - -You’ve probably learned to solve a system of equations using methods such as substitution, graphical solutions, or elimination. These approaches work well for small systems, but they can become tedious as the number of equations grows. The question is: *can we find a way to carry out elimination in an algorithmic way?* - -When solving systems of equations by hand, one strategy is to eliminate variables step by step until you are left with a simpler equation that can be solved directly. Forward elimination is the algorithmic version of this same process. - -The idea is to start with the first equation, we use it to eliminate the first variable from all equations below it. Then we move to the second equation, using it to eliminate the second variable from all equations below it, and so on. This process continues until the coefficient matrix has been transformed into an upper-triangular form (all zeros below the diagonal). - -For example, suppose we have a 3×3 system of equations. After forward elimination, it will look like: -$$ -\begin{bmatrix} a_{11} & a_{12} & a_{13} & | & b_1 \\ 0 & a_{22} & a_{23} & | & b_2 \\ 0 & 0 & a_{33} & | & b_3 \\ \end{bmatrix} $$ -Notice how everything below the diagonal has been reduced to zero. At this point, the system is much easier to solve. - -Notice how everything below the diagonal has been reduced to zero. At this point, the system is much easier to solve because each row now involves fewer unknowns as you move downward. The last equation contains only one variable, which can be solved directly: -$$ -a_{nn} x_n = b_n \quad \Rightarrow \quad x_n = \frac{b_n}{a_{nn}} -$$ -Once $x_n$ is known, it can be substituted into the equation above it to solve for $x_{n-1}$: -$$ -a_{(n-1)(n-1)}x_{n-1} + a_{(n-1)n}x_n = b_{n-1} \quad \Rightarrow \quad x_{n-1} = \frac{b_{n-1} - a_{(n-1)n}x_n}{a_{(n-1)(n-1)}} -$$ -This process continues row by row, moving upward through the system. In general, the formula for the $m$-th variable is: -$$ -x_m = \frac{b_m - \sum_{j=m+1}^n a_{mj}x_j}{a_{mm}} -$$ -This step-by-step procedure is known as **backward substitution**, and together with forward elimination it forms the backbone of systematic solution methods like Gaussian elimination. -<img - style="display: block; - margin-left: auto; - margin-right: auto; - width: 60%;" - src="fw_elim_bw_sub.png" - alt="Forward Elimination and Backwards Substitution"> -Now that we have an algorithm to follow, let's code it. -```python -import numpy as np - -# Create the Augmented matrix -A = np.array([[ 2.2, 1.5, -1.3, 8.2], - [-3.4, -1.7, 2.8, -11.1], - [-2.0, 1.0, 2.4, -3.5]]) - - -``` -### Forward Elimination -```python -n = len(A) - -# Forward elimination loop -for i in range(n-1): - for j in range(i+1, n): - factor = A[j, i] / A[i, i] # multiplier to eliminate A[j,i] - A[j, i:] = A[j, i:] - factor * A[i, i:] - -print("Upper triangular matrix after forward elimination:") -print(A) -``` -### Back Substitution -```python -x = np.zeros(n) - -for i in range(n-1, -1, -1): - x[i] = (A[i, -1] - np.dot(A[i, i+1:n], x[i+1:n])) / A[i, i] - -print("Solution:", x) - -``` - -### Naive Gauss Elimination -We can write a function combining these two steps: -```python -def gaussian_elimination_naive(A, b): - A = A.astype(float).copy() - b = b.astype(float).copy() - n = len(b) - # Forward elimination - for k in range(n-1): - if np.isclose(A[k,k], 0.0): - raise ZeroDivisionError("Zero pivot encountered (naïve). Use pivoting.") - for i in range(k+1, n): - m = A[i,k]/A[k,k] - A[i, k:] -= m * A[k, k:] - b[i] -= m * b[k] - # Back substitution - x = np.zeros(n) - for i in range(n-1, -1, -1): - s = A[i, i+1:] @ x[i+1:] - x[i] = (b[i]-s)/A[i,i] - return x -``` -### Gauss Elimination -```python -def gaussian_elimination_pp(A, b): - A = A.astype(float).copy() - b = b.astype(float).copy() - n = len(b) - # Forward elimination with partial pivoting - for k in range(n-1): - piv = k + np.argmax(np.abs(A[k:, k])) - if np.isclose(A[piv, k], 0.0): - raise np.linalg.LinAlgError("Matrix is singular or nearly singular.") - if piv != k: # row swap - A[[k, piv]] = A[[piv, k]] - b[[k, piv]] = b[[piv, k]] - for i in range(k+1, n): - m = A[i,k]/A[k,k] - A[i, k:] -= m * A[k, k:] - b[i] -= m * b[k] - # Back substitution - x = np.zeros(n) - for i in range(n-1, -1, -1): - s = A[i, i+1:] @ x[i+1:] - x[i] = (b[i]-s)/A[i,i] - return x -``` - - -# LU Decomposition -Imagine you’re designing the heat shield of a spacecraft. Engineers provide you with a structural model of the heat shield, represented by a coefficient matrix $A$. To evaluate how the design performs under different operating conditions, you need to test many different load cases, each represented by a right-hand side vector $\mathbf{b}$. - -Now suppose your model has a 50x50 matrix $A$, and you want to simulate 100 different load cases. If you rely on Gaussian elimination alone, you would need to **repeat the elimination process 100 times** - once for each $\mathbf{b}$. This quickly becomes computationally expensive, especially for larger systems. - -This is where **LU decomposition** comes in. Instead of redoing elimination for every load case, we can factorize the matrix $A$ **just once** into two simpler pieces: -$$ -A = LU -$$ -expanding it: -$$ -\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix} -= -\begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \\ \end{bmatrix} -\begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \\ \end{bmatrix} -$$ -Once we have this factorization, solving $A\mathbf{x} = \mathbf{b}$ becomes a two-step process: - -1. Solve $L\mathbf{d} = \mathbf{b}$ by **forward substitution**. -$$ -d_m \;=\; \frac{\,b_m \;-\; \sum_{n=1}^{m-1} a_{mn}\, d_n\,}{a_{mm}} -$$ -2. Solve $U\mathbf{x} = \mathbf{d}$ by **backward substitution**. - $$ - x_m=\frac{\,d_m-\sum_{n=m+1}^{N} u_{mn}\,x_n\,}{u_{mm}} - $$ - -The big advantage is that the factorization is reused. For 100 different load vectors, we only do the expensive decomposition once, and then solve two triangular systems (fast operations) for each new $\mathbf{b}$. -<img - style="display: block; - margin-left: auto; - margin-right: auto; - width: 40%;" - src="https://upload.wikimedia.org/wikipedia/commons/thumb/b/bd/Tadeusz_Banachiewicz_%28NAC%29.jpg/1024px-Tadeusz_Banachiewicz_%28NAC%29.jpg" - alt="Tadeusz Banchiewicz"> -First formalized in 1938 by Tadeusz Banachiewicz, based on Gauss’s elimination ideas. - -Example - - -<img - style="display: block; - margin-left: auto; - margin-right: auto; - width: 75%;" - src="lu_decomp.png" - alt="LU Decomposition"> -## Problem 1 -```python -""" -LU Decomposition Solver (with partial pivoting) - -Solves A x = b by computing P, L, U such that P A = L U, -then performing: - 1) y from L y = P b (forward substitution) - 2) x from U x = y (back substitution) - -- Works for a single RHS (b shape: (n,)) or multiple RHS (b shape: (n, m)) -- Raises a ValueError if a near-zero pivot is encountered. -""" - -import numpy as np - -def lu_decompose(A, pivot_tol=1e-12): - """ - Compute LU factorization with partial pivoting. - Returns P, L, U such that P @ A = L @ U. - - Parameters - ---------- - A : (n, n) ndarray - pivot_tol : float - If the absolute value of the pivot falls below this, the matrix is - treated as singular. - - Returns - ------- - P, L, U : ndarrays of shape (n, n) - """ - A = np.array(A, dtype=float, copy=True) - n = A.shape[0] - if A.shape[0] != A.shape[1]: - raise ValueError("A must be square") - - U = A.copy() - L = np.eye(n) - P = np.eye(n) - - for k in range(n - 1): - # Partial pivoting: find index of largest |U[i,k]| for i >= k - pivot_idx = np.argmax(np.abs(U[k:, k])) + k - pivot_val = U[pivot_idx, k] - - if abs(pivot_val) < pivot_tol: - raise ValueError(f"Numerically singular matrix: pivot ~ 0 at column {k}") - - # Swap rows in U and P (and the part of L constructed so far) - if pivot_idx != k: - U[[k, pivot_idx], :] = U[[pivot_idx, k], :] - P[[k, pivot_idx], :] = P[[pivot_idx, k], :] - if k > 0: - L[[k, pivot_idx], :k] = L[[pivot_idx, k], :k] - - # Gaussian elimination to build L and zero out below the pivot in U - for i in range(k + 1, n): - L[i, k] = U[i, k] / U[k, k] - U[i, k:] -= L[i, k] * U[k, k:] - U[i, k] = 0.0 - - # Final pivot check on U[-1, -1] - if abs(U[-1, -1]) < pivot_tol: - raise ValueError(f"Numerically singular matrix: pivot ~ 0 at last diagonal") - - return P, L, U - -def forward_substitution(L, b): - """ - Solve L y = b for y, where L is unit lower-triangular. - b can be (n,) or (n, m). - """ - L = np.asarray(L) - b = np.asarray(b, dtype=float) - n = L.shape[0] - - if b.ndim == 1: - b = b.reshape(n, 1) - - y = np.zeros_like(b, dtype=float) - for i in range(n): - # L has ones on the diagonal (Doolittle), so we don't divide by L[i,i] - y[i] = b[i] - L[i, :i] @ y[:i] - return y if y.shape[1] > 1 else y.ravel() - -def back_substitution(U, y): - """ - Solve U x = y for x, where U is upper-triangular. - y can be (n,) or (n, m). - """ - U = np.asarray(U) - y = np.asarray(y, dtype=float) - n = U.shape[0] - - if y.ndim == 1: - y = y.reshape(n, 1) - - x = np.zeros_like(y, dtype=float) - for i in range(n - 1, -1, -1): - if abs(U[i, i]) < 1e-15: - raise ValueError(f"Zero (or tiny) pivot encountered at U[{i},{i}]") - x[i] = (y[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i] - return x if x.shape[1] > 1 else x.ravel() - -def lu_solve(A, b, pivot_tol=1e-12): - """ - Solve A x = b using LU with partial pivoting. - Returns x and the factors (P, L, U). - """ - P, L, U = lu_decompose(A, pivot_tol=pivot_tol) - Pb = P @ np.asarray(b, dtype=float) - y = forward_substitution(L, Pb) - x = back_substitution(U, y) - return x, (P, L, U) - -def main(): - # Example system: - # 3x + 2y - z = 1 - # 2x - 2y + 4z = -2 - # -x + 0.5y - z = 0 - A = np.array([ - [ 3.0, 2.0, -1.0], - [ 2.0, -2.0, 4.0], - [-1.0, 0.5, -1.0] - ]) - b = np.array([1.0, -2.0, 0.0]) - - x, (P, L, U) = lu_solve(A, b) - print("Solution x:\n", x) - print("\nP:\n", P) - print("\nL:\n", L) - print("\nU:\n", U) - - # Residual check - r = A @ x - b - print("\nResidual A@x - b:\n", r) - print("Residual norm:", np.linalg.norm(r, ord=np.inf)) - -if __name__ == "__main__": - main() -``` - -## Problem 2 -Ohm's and Kirchhoff's laws have been used to find a set of equations to model a circuit. Using the built-in Numpy solvers. Plot the function of current at different voltages. - - -```python -import numpy as np -import matplotlib.pyplot as plt - -# Unknown ordering -VR1, VR2, VR3, VR4, VR5, I1, I2, IM1, IM2, IM3 = range(10) - -def build_system_matrices(R1, R2, R3, R4, R5): - A = np.zeros((10,10), float) - b_const = np.zeros(10, float) - - # KVL/KCL + Ohm's law equations - A[0, VR1] = -1; A[0, VR2] = -1 # -VR1 - VR2 = -Vs + 5 - A[1, VR2] = 1; A[1, VR3] = -1; A[1, VR4] = -1; b_const[1] = -2 # VR2 - VR3 - VR4 = -2 - A[2, VR4] = 1; A[2, VR5] = -1; b_const[2] = 2 # VR4 - VR5 = 2 - - A[3, VR1] = 1; A[3, I1 ] = -R1 # VR1 = I1 R1 - A[4, VR2] = 1; A[4, IM1] = -R2 # VR2 = IM1 R2 - A[5, VR3] = 1; A[5, I2 ] = -R3 # VR3 = I2 R3 - A[6, VR4] = 1; A[6, IM2] = -R4 # VR4 = IM2 R4 - A[7, VR5] = 1; A[7, IM3] = -R5 # VR5 = IM3 R5 - - A[8, I1 ] = 1; A[8, IM1] = -1; A[8, I2 ] = -1 # I1 = IM1 + I2 - A[9, I2 ] = 1; A[9, IM2] = -1; A[9, IM3] = -1 # I2 = IM2 + IM3 - return A, b_const - -def b_vector(Vs, b_const): - b = b_const.copy() - b[0] = -Vs + 5.0 - return b - -# Example component values -R1, R2, R3, R4, R5 = 10.0, 5.0, 8.0, 6.0, 4.0 - -A, b_const = build_system_matrices(R1, R2, R3, R4, R5) - -# Sanity: ensure A is non-singular -condA = np.linalg.cond(A) -condA - -Vs_vals = np.arange(0.0, 48.0 + 0.1, 0.1) - -# Option 1 (straightforward): loop with numpy.linalg.solve (LAPACK under the hood) -IM1_vals = np.empty_like(Vs_vals) -IM2_vals = np.empty_like(Vs_vals) -IM3_vals = np.empty_like(Vs_vals) - -for k, Vs in enumerate(Vs_vals): - b = b_vector(Vs, b_const) - x = np.linalg.solve(A, b) # Uses LU factorization internally - IM1_vals[k], IM2_vals[k], IM3_vals[k] = x[IM1], x[IM2], x[IM3] - -# Plot -plt.figure() -plt.plot(Vs_vals, IM1_vals, label="IM1 (solve)") -plt.plot(Vs_vals, IM2_vals, label="IM2 (solve)") -plt.plot(Vs_vals, IM3_vals, label="IM3 (solve)") -plt.xlabel("Supply Voltage Vs (V)") -plt.ylabel("Motor Current (A)") -plt.title(f"Motor Currents vs Supply Voltage (cond(A) ≈ {condA:.1f}, diff≈{max_diff:.2e})") -plt.legend() -plt.show() - -```
\ No newline at end of file |