Feat: hw2 done

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "docs/template"]
+	path = docs/template
+	url = ssh://gitea@konchin.com:2222/ytshih/typst-template.git

HW2_111550013/HW2_111550013_code/problem1.cpp

@@ -0,0 +1,49 @@
+// Problem 1: Gaussian elimination with partial pivoting
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+using namespace std;
+
+const int N = 4;
+
+int main() {
+    double A[N][N+1] = {
+        {4, 2, -2, -1, 7},
+        {0, 4, 1, 2, 10},
+        {3, -2, 1, 2, 2},
+        {2, 0, 3, -5, 3}
+    };
+
+    // Forward elimination with partial pivoting
+    for (int k = 0; k < N - 1; k++) {
+        // Find pivot
+        int maxRow = k;
+        for (int i = k + 1; i < N; i++)
+            if (fabs(A[i][k]) > fabs(A[maxRow][k])) maxRow = i;
+
+        if (maxRow != k) {
+            cout << "Column " << (k+1) << ": swap row " << (k+1) << " with row " << (maxRow+1) << "\n";
+            for (int j = 0; j <= N; j++) swap(A[k][j], A[maxRow][j]);
+        }
+
+        // Eliminate
+        for (int i = k + 1; i < N; i++) {
+            double m = A[i][k] / A[k][k];
+            for (int j = k; j <= N; j++) A[i][j] -= m * A[k][j];
+        }
+    }
+
+    // Back substitution
+    double x[N];
+    for (int i = N - 1; i >= 0; i--) {
+        x[i] = A[i][N];
+        for (int j = i + 1; j < N; j++) x[i] -= A[i][j] * x[j];
+        x[i] /= A[i][i];
+    }
+
+    cout << "\nSolution:\n";
+    for (int i = 0; i < N; i++)
+        cout << "x" << (i+1) << " = " << fixed << setprecision(6) << x[i] << "\n";
+
+    return 0;
+}

HW2_111550013/HW2_111550013_code/problem2.cpp

@@ -0,0 +1,62 @@
+// Problem 2: Scaled partial pivoting
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <vector>
+using namespace std;
+
+double roundSig(double val, int sig) {
+    if (val == 0) return 0;
+    double d = ceil(log10(fabs(val)));
+    double p = pow(10, sig - d);
+    return round(val * p) / p;
+}
+
+void solve(vector<vector<double>> A, int n, int sig) {
+    cout << "\n" << sig << " significant digits:\n";
+
+    vector<double> s(n);
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < n; j++) s[i] = max(s[i], fabs(A[i][j]));
+
+    vector<int> idx(n);
+    for (int i = 0; i < n; i++) idx[i] = i;
+
+    for (int k = 0; k < n - 1; k++) {
+        int maxRow = k;
+        double maxRatio = fabs(A[idx[k]][k]) / s[idx[k]];
+        for (int i = k + 1; i < n; i++) {
+            double ratio = fabs(A[idx[i]][k]) / s[idx[i]];
+            if (ratio > maxRatio) { maxRatio = ratio; maxRow = i; }
+        }
+        swap(idx[k], idx[maxRow]);
+
+        for (int i = k + 1; i < n; i++) {
+            double m = roundSig(A[idx[i]][k] / A[idx[k]][k], sig);
+            for (int j = k; j <= n; j++)
+                A[idx[i]][j] = roundSig(A[idx[i]][j] - m * A[idx[k]][j], sig);
+        }
+    }
+
+    vector<double> x(n);
+    for (int i = n - 1; i >= 0; i--) {
+        x[i] = A[idx[i]][n];
+        for (int j = i + 1; j < n; j++)
+            x[i] = roundSig(x[i] - A[idx[i]][j] * x[j], sig);
+        x[i] = roundSig(x[i] / A[idx[i]][i], sig);
+    }
+
+    for (int i = 0; i < n; i++)
+        cout << "x" << (i+1) << " = " << fixed << setprecision(6) << x[i] << "\n";
+}
+
+int main() {
+    vector<vector<double>> A = {
+        {4.13, -2.20, 0.95, 3.02},
+        {6.14, 4.45, -1.45, -4.02},
+        {1.03, 1.86, 0.44, 5.22}
+    };
+    solve(A, 3, 6);
+    solve(A, 3, 3);
+    return 0;
+}

HW2_111550013/HW2_111550013_code/problem3.cpp

@@ -0,0 +1,32 @@
+// Problem 3: Symmetric tridiagonal system
+#include <iostream>
+#include <iomanip>
+#include <vector>
+using namespace std;
+
+int main() {
+    int n = 6;
+    vector<double> d = {4, 4, 4, 4, 4, 4};      // diagonal
+    vector<double> a = {-1, -1, -1, -1, -1};    // off-diagonal
+    vector<double> b = {100, 200, 200, 200, 200, 100};
+
+    // Forward elimination
+    for (int i = 1; i < n; i++) {
+        double m = a[i-1] / d[i-1];
+        d[i] -= m * a[i-1];
+        b[i] -= m * b[i-1];
+    }
+
+    // Back substitution
+    vector<double> x(n);
+    x[n-1] = b[n-1] / d[n-1];
+    for (int i = n - 2; i >= 0; i--)
+        x[i] = (b[i] - a[i] * x[i+1]) / d[i];
+
+    cout << "Solution:\n";
+    for (int i = 0; i < n; i++)
+        cout << "x" << (i+1) << " = " << fixed << setprecision(6) << x[i] << "\n";
+
+    cout << "\nOperations: 7N - 6 = " << (7*n - 6) << " for N = " << n << "\n";
+    return 0;
+}

HW2_111550013/HW2_111550013_code/problem4.cpp

@@ -0,0 +1,57 @@
+// Problem 4: Jacobi and Gauss-Seidel methods
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+using namespace std;
+
+const double TOL = 1e-6;
+const int MAX_ITER = 100;
+
+int main() {
+    // Rearranged system (weakly diagonally dominant):
+    // 7x1 - 3x2 + 4x3 = 6
+    // 2x1 + 5x2 + 3x3 = -5
+    // -3x1 + 2x2 + 6x3 = 2
+
+    double b[3] = {6, -5, 2};
+
+    // Jacobi
+    cout << "Jacobi Method:\n";
+    vector<double> x = {0, 0, 0};
+    for (int iter = 1; iter <= MAX_ITER; iter++) {
+        vector<double> x_old = x;
+        x[0] = (b[0] + 3*x_old[1] - 4*x_old[2]) / 7;
+        x[1] = (b[1] - 2*x_old[0] - 3*x_old[2]) / 5;
+        x[2] = (b[2] + 3*x_old[0] - 2*x_old[1]) / 6;
+
+        double maxDiff = max({fabs(x[0]-x_old[0]), fabs(x[1]-x_old[1]), fabs(x[2]-x_old[2])});
+        if (maxDiff < TOL) {
+            cout << "Converged in " << iter << " iterations\n";
+            break;
+        }
+    }
+    for (int i = 0; i < 3; i++)
+        cout << "x" << (i+1) << " = " << fixed << setprecision(6) << x[i] << "\n";
+
+    // Gauss-Seidel
+    cout << "\nGauss-Seidel Method:\n";
+    x = {0, 0, 0};
+    for (int iter = 1; iter <= MAX_ITER; iter++) {
+        vector<double> x_old = x;
+        x[0] = (b[0] + 3*x[1] - 4*x[2]) / 7;
+        x[1] = (b[1] - 2*x[0] - 3*x[2]) / 5;
+        x[2] = (b[2] + 3*x[0] - 2*x[1]) / 6;
+
+        double maxDiff = max({fabs(x[0]-x_old[0]), fabs(x[1]-x_old[1]), fabs(x[2]-x_old[2])});
+        if (maxDiff < TOL) {
+            cout << "Converged in " << iter << " iterations\n";
+            break;
+        }
+    }
+    for (int i = 0; i < 3; i++)
+        cout << "x" << (i+1) << " = " << fixed << setprecision(6) << x[i] << "\n";
+
+    return 0;
+}

docs/main.typ

@@ -0,0 +1,288 @@
+#import "template/module.typ": *
+
+#show: default.with(
+  title: "Numerical Methods Assignment 2",
+)
+
+= Gaussian Elimination with Partial Pivoting (15%)
+
+Solve the following system by Gaussian elimination with partial pivoting:
+$ mat(4, 2, -2, -1; 0, 4, 1, 2; 3, -2, 1, 2; 2, 0, 3, -5) vec(x_1, x_2, x_3, x_4) = vec(7, 10, 2, 3) $
+
+== Solution Process
+
+*Initial augmented matrix:*
+$ mat(augment: #4, 4, 2, -2, -1, 7; 0, 4, 1, 2, 10; 3, -2, 1, 2, 2; 2, 0, 3, -5, 3) $
+
+*Step 1 - Column 1:* The pivot element is $a_(11) = 4$. Checking column 1: $|4| >= |0|, |3|, |2|$. No row interchange needed.
+
+Multipliers: $m_(21) = 0$, $m_(31) = 3/4 = 0.75$, $m_(41) = 2/4 = 0.5$
+
+After elimination:
+$ mat(augment: #4, 4, 2, -2, -1, 7; 0, 4, 1, 2, 10; 0, -3.5, 2.5, 2.75, -3.25; 0, -1, 4, -4.5, -0.5) $
+
+*Step 2 - Column 2:* The pivot element is $a_(22) = 4$. Checking: $|4| >= |-3.5|, |-1|$. No row interchange needed.
+
+Multipliers: $m_(32) = -3.5/4 = -0.875$, $m_(42) = -1/4 = -0.25$
+
+After elimination:
+$ mat(augment: #4, 4, 2, -2, -1, 7; 0, 4, 1, 2, 10; 0, 0, 3.375, 4.5, 5.5; 0, 0, 4.25, -4, 2) $
+
+*Step 3 - Column 3:* Checking: $|3.375| < |4.25|$. *Row interchange required: swap row 3 with row 4.*
+
+After interchange:
+$ mat(augment: #4, 4, 2, -2, -1, 7; 0, 4, 1, 2, 10; 0, 0, 4.25, -4, 2; 0, 0, 3.375, 4.5, 5.5) $
+
+Multiplier: $m_(43) = 3.375/4.25 = 0.7941$
+
+After elimination:
+$ mat(augment: #4, 4, 2, -2, -1, 7; 0, 4, 1, 2, 10; 0, 0, 4.25, -4, 2; 0, 0, 0, 7.6765, 3.9118) $
+
+*Back substitution:*
+- $x_4 = 3.9118 / 7.6765 = 0.509579$
+- $x_3 = (2 - (-4)(0.509579)) / 4.25 = 0.950192$
+- $x_2 = (10 - (1)(0.950192) - (2)(0.509579)) / 4 = 2.007663$
+- $x_1 = (7 - (2)(2.007663) - (-2)(0.950192) - (-1)(0.509579)) / 4 = 1.348659$
+
+== Solution
+$ x_1 = 1.348659, quad x_2 = 2.007663, quad x_3 = 0.950192, quad x_4 = 0.509579 $
+
+*Column requiring row interchange:* Column 3
+
+#pagebreak()
+
+= Scaled Partial Pivoting (25%)
+
+Solve the system using scaled partial pivoting:
+$ mat(4.13, -2.20, 0.95, 3.02; 6.14, 4.45, -1.45, -4.02; 1.03, 1.86, 0.44, 5.22) $
+
+== Part (a): Six Significant Digits
+
+*Scale factors:*
+- $s_1 = max(|4.13|, |-2.20|, |0.95|) = 4.13$
+- $s_2 = max(|6.14|, |4.45|, |-1.45|) = 6.14$
+- $s_3 = max(|1.03|, |1.86|, |0.44|) = 1.86$
+
+*Step 1 - Column 1:*
+
+Scaled ratios:
+- Row 1: $|4.13|/4.13 = 1.000$
+- Row 2: $|6.14|/6.14 = 1.000$
+- Row 3: $|1.03|/1.86 = 0.554$
+
+Rows 1 and 2 are tied; keep current order.
+
+Multipliers: $m_(21) = 6.14/4.13 = 1.48668$, $m_(31) = 1.03/4.13 = 0.249395$
+
+After elimination:
+$ mat(augment: #3, 4.13, -2.20, 0.95, 3.02; 0, 7.7207, -2.86235, -8.50977; 0, 2.40867, 0.203075, 4.46683) $
+
+*Step 2 - Column 2:*
+
+Scaled ratios:
+- Row 2: $|7.7207|/6.14 = 1.257$
+- Row 3: $|2.40867|/1.86 = 1.295$
+
+Row 3 has larger ratio. *Interchange rows 2 and 3.*
+
+Multiplier: $m_(22) = 7.7207/2.40867 = 3.20538$
+
+After elimination:
+$ mat(augment: #3, 4.13, -2.20, 0.95, 3.02; 0, 2.40867, 0.203075, 4.46683; 0, 0, -3.51328, -22.8277) $
+
+*Back substitution:*
+- $x_3 = -22.8277 / (-3.51328) = 6.49755$
+- $x_2 = (4.46683 - 0.203075 times 6.49755) / 2.40867 = 1.30667$
+- $x_1 = (3.02 - 0.95 times 6.49755 - (-2.20) times 1.30667) / 4.13 = -0.0673$
+
+*Solution (6 significant digits):*
+$ x_1 = -0.0673, quad x_2 = 1.3067, quad x_3 = 6.4976 $
+
+== Part (b): Three Significant Digits
+
+Following the same process with 3-digit arithmetic:
+
+*Solution (3 significant digits):*
+$ x_1 = -0.0688, quad x_2 = 1.31, quad x_3 = 6.51 $
+
+== Comparison
+
+The solutions are close but show small differences due to rounding errors:
+- $x_1$: $-0.0673$ vs $-0.0688$ (difference $approx 2%$)
+- $x_2$: $1.3067$ vs $1.31$ (difference $approx 0.3%$)
+- $x_3$: $6.4976$ vs $6.51$ (difference $approx 0.2%$)
+
+The solution with 3 significant digits is slightly different but still reasonable. The accumulated rounding errors affect $x_1$ the most since it is computed last in back substitution.
+
+#pagebreak()
+
+= Symmetric Tridiagonal System (30%)
+
+The system:
+$ mat(4, -1, 0, 0, 0, 0; -1, 4, -1, 0, 0, 0; 0, -1, 4, -1, 0, 0; 0, 0, -1, 4, -1, 0; 0, 0, 0, -1, 4, -1; 0, 0, 0, 0, -1, 4) vec(x_1, x_2, x_3, x_4, x_5, x_6) = vec(100, 200, 200, 200, 200, 100) $
+
+== Part (a): Algorithm for Symmetric Tridiagonal System
+
+For a symmetric tridiagonal matrix, we only need to store:
+- $d_i$: diagonal elements ($n$ elements)
+- $a_i$: off-diagonal elements ($n-1$ elements, same above and below diagonal)
+- $b_i$: right-hand side ($n$ elements)
+
+*Algorithm (Modified Thomas Algorithm for Symmetric Case):*
+
+```
+// Input: d[1..n] (diagonal), a[1..n-1] (off-diagonal), b[1..n] (RHS)
+// Output: x[1..n] (solution)
+
+// Forward Elimination
+for i = 2 to n:
+    m = a[i-1] / d[i-1]
+    d[i] = d[i] - m * a[i-1]
+    b[i] = b[i] - m * b[i-1]
+
+// Back Substitution
+x[n] = b[n] / d[n]
+for i = n-1 downto 1:
+    x[i] = (b[i] - a[i] * x[i+1]) / d[i]
+```
+
+== Part (b): Solving the System
+
+*Initial compact representation:*
+#align(center)[
+  #table(
+    columns: 4,
+    [Row], [Sub-diagonal], [Diagonal], [RHS],
+    [1], [-], [4], [100],
+    [2], [-1], [4], [200],
+    [3], [-1], [4], [200],
+    [4], [-1], [4], [200],
+    [5], [-1], [4], [200],
+    [6], [-1], [4], [100],
+  )
+]
+
+*Forward Elimination:*
+- Step 1: $m = -1/4 = -0.25$, $d'_2 = 4 - (-0.25)(-1) = 3.75$, $b'_2 = 200 - (-0.25)(100) = 225$
+- Step 2: $m = -1/3.75 = -0.2667$, $d'_3 = 3.7333$, $b'_3 = 260$
+- Step 3: $m = -1/3.7333 = -0.2679$, $d'_4 = 3.7321$, $b'_4 = 269.64$
+- Step 4: $m = -1/3.7321 = -0.2679$, $d'_5 = 3.7321$, $b'_5 = 272.25$
+- Step 5: $m = -1/3.7321 = -0.2679$, $d'_6 = 3.7321$, $b'_6 = 172.95$
+
+*Back Substitution:*
+- $x_6 = 172.95/3.7321 = 46.34$
+- $x_5 = (272.25 - (-1)(46.34))/3.7321 = 85.37$
+- $x_4 = (269.64 - (-1)(85.37))/3.7321 = 95.12$
+- $x_3 = (260 - (-1)(95.12))/3.7333 = 95.12$
+- $x_2 = (225 - (-1)(95.12))/3.75 = 85.37$
+- $x_1 = (100 - (-1)(85.37))/4 = 46.34$
+
+*Solution:*
+$ x_1 = 46.34, quad x_2 = 85.37, quad x_3 = 95.12, quad x_4 = 95.12, quad x_5 = 85.37, quad x_6 = 46.34 $
+
+== Part (c): Arithmetic Operations Count
+
+For a system of $N$ equations:
+
+*Forward Elimination (loop runs $N-1$ times):*
+- $N-1$ divisions (computing $m$)
+- $N-1$ multiplications (computing $m dot a_(i-1)$)
+- $2(N-1)$ subtractions (updating $d_i$ and $b_i$)
+
+Total: $4(N-1)$ operations
+
+*Back Substitution:*
+- 1 division (computing $x_N$)
+- For $i = N-1$ down to 1: 1 multiplication + 1 subtraction + 1 division = $3(N-1)$ operations
+
+Total: $1 + 3(N-1) = 3N - 2$ operations
+
+*Total operations:* $4(N-1) + 3N - 2 = 7N - 6 = O(N)$
+
+This is much more efficient than general Gaussian elimination which requires $O(N^3)$ operations.
+
+#pagebreak()
+
+= Jacobi and Gauss-Seidel Methods (40%)
+
+Original system:
+$ mat(7, -3, 4; -3, 2, 6; 2, 5, 3) vec(x_1, x_2, x_3) = vec(6, 2, -5) $
+
+== Checking and Rearranging for Diagonal Dominance
+
+*Original matrix diagonal dominance check:*
+- Row 1: $|7| = 7$, $|-3| + |4| = 7$ (equal, not strictly dominant)
+- Row 2: $|2| = 2$, $|-3| + |6| = 9$ (NOT dominant)
+- Row 3: $|3| = 3$, $|2| + |5| = 7$ (NOT dominant)
+
+*Rearrangement:* Swap equations 2 and 3:
+$ mat(7, -3, 4; 2, 5, 3; -3, 2, 6) vec(x_1, x_2, x_3) = vec(6, -5, 2) $
+
+After rearrangement:
+- Row 1: $|7| = 7$, $|-3| + |4| = 7$ (weakly dominant)
+- Row 2: $|5| = 5$, $|2| + |3| = 5$ (weakly dominant)
+- Row 3: $|6| = 6 > |-3| + |2| = 5$ (strictly dominant)
+
+The system is weakly diagonally dominant. The methods may still converge.
+
+== Part (a): Jacobi Method
+
+*Iteration formulas:*
+$ x_1^((k+1)) &= (6 + 3x_2^((k)) - 4x_3^((k))) / 7 \
+x_2^((k+1)) &= (-5 - 2x_1^((k)) - 3x_3^((k))) / 5 \
+x_3^((k+1)) &= (2 + 3x_1^((k)) - 2x_2^((k))) / 6 $
+
+*Starting vector:* $bold(x)^((0)) = [0, 0, 0]^T$
+
+#align(center)[
+  #table(
+    columns: 4,
+    [Iteration], [$x_1$], [$x_2$], [$x_3$],
+    [0], [0], [0], [0],
+    [1], [0.8571], [-1.0000], [0.3333],
+    [2], [0.2381], [-1.5429], [1.0952],
+    [3], [-0.4299], [-1.7524], [0.9667],
+    [4], [-0.4463], [-1.4080], [0.7025],
+    [5], [-0.1477], [-1.2430], [0.5795],
+    [...], [...], [...], [...],
+    [37], [-0.1433], [-1.3746], [0.7199],
+  )
+]
+
+*Jacobi converged after 37 iterations.*
+
+*Solution:* $x_1 = -0.1433, quad x_2 = -1.3746, quad x_3 = 0.7199$
+
+== Part (b): Gauss-Seidel Method
+
+*Iteration formulas (using most recent values):*
+$ x_1^((k+1)) &= (6 + 3x_2^((k)) - 4x_3^((k))) / 7 \
+x_2^((k+1)) &= (-5 - 2x_1^((k+1)) - 3x_3^((k))) / 5 \
+x_3^((k+1)) &= (2 + 3x_1^((k+1)) - 2x_2^((k+1))) / 6 $
+
+#align(center)[
+  #table(
+    columns: 4,
+    [Iteration], [$x_1$], [$x_2$], [$x_3$],
+    [0], [0], [0], [0],
+    [1], [0.8571], [-1.3429], [1.2095],
+    [2], [-0.4095], [-1.5619], [0.6492],
+    [3], [-0.1832], [-1.3162], [0.6805],
+    [4], [-0.0958], [-1.3700], [0.7421],
+    [5], [-0.1540], [-1.3836], [0.7175],
+    [...], [...], [...], [...],
+    [16], [-0.1433], [-1.3746], [0.7199],
+  )
+]
+
+*Gauss-Seidel converged after 16 iterations.*
+
+*Solution:* $x_1 = -0.1433, quad x_2 = -1.3746, quad x_3 = 0.7199$
+
+== Comparison
+
+Both methods converge to the same solution:
+$ x_1 approx -0.143, quad x_2 approx -1.375, quad x_3 approx 0.720 $
+
+Gauss-Seidel converged in *16 iterations*, while Jacobi required *37 iterations*. This is expected because Gauss-Seidel uses updated values immediately, leading to faster convergence.
+