diff --git a/.gitignore b/.gitignore index 6d0ee45..da5ce1f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ .vscode -.DS_Store \ No newline at end of file +.DS_Store +__pycache__ +*.pyc diff --git a/project2/iterative_solvers.py b/project2/iterative_solvers.py new file mode 100644 index 0000000..944276f --- /dev/null +++ b/project2/iterative_solvers.py @@ -0,0 +1,154 @@ +import numpy as np +from typing import Callable, Any, List, Tuple + +# Tolerance and maximum iterations +TOL = 1e-10 +MAX_ITER = 20000 + +def stationary_method( + b: np.ndarray, + method: str, + get_diag_inv: Callable[[int, Any], float], + get_off_diag_sum: Callable[[np.ndarray, int, Any], float], + mat_vec_prod: Callable[[np.ndarray, Any], np.ndarray], + omega: float = 1.0, + **kwargs + ) -> Tuple[np.ndarray, int, List[float]]: + """ + Solves Ax=b using stationary methods (Jacobi, Gauss-Seidel, SOR). + """ + n = len(b) + x = np.zeros(n) + r0 = b.copy() + r0_norm = np.linalg.norm(r0) + if r0_norm == 0: + return x, 0, [0.0] + + residual_history = [] + + for k in range(MAX_ITER): + # 1. Calculate current residual and its norm + r = b - mat_vec_prod(x, **kwargs) + res_norm = np.linalg.norm(r) + + # 2. Append to history + residual_history.append(res_norm) + + # 3. Check for convergence + if res_norm / r0_norm < TOL: + return x, k, residual_history + + # 4. If not converged, perform the next iteration update + if method == 'jacobi': + x_old = x.copy() + for i in range(n): + off_diag_val = get_off_diag_sum(x_old, i, **kwargs) + diag_inv_val = get_diag_inv(i, **kwargs) + x[i] = diag_inv_val * (b[i] - off_diag_val) + else: # Gauss-Seidel and SOR + x_k_iter = x.copy() # Needed for SOR + for i in range(n): + off_diag_val = get_off_diag_sum(x, i, **kwargs) + diag_inv_val = get_diag_inv(i, **kwargs) + x_gs = diag_inv_val * (b[i] - off_diag_val) + + if method == 'gauss_seidel': + x[i] = x_gs + elif method == 'sor': + x[i] = (1 - omega) * x_k_iter[i] + omega * x_gs + + # If loop finishes (MAX_ITER reached), calculate and append final residual + r_final = b - mat_vec_prod(x, **kwargs) + residual_history.append(np.linalg.norm(r_final)) + return x, MAX_ITER, residual_history + + +def steepest_descent( + mat_vec_prod: Callable[[np.ndarray, Any], np.ndarray], + b: np.ndarray, + **kwargs + ) -> Tuple[np.ndarray, int, List[float]]: + """ + Solves Ax=b using Steepest Descent, with corrected residual logging. + """ + x = np.zeros_like(b) + r = b.copy() + + r0_norm = np.linalg.norm(r) + if r0_norm == 0: + return x, 0, [0.0] + + residual_history = [] + + for k in range(MAX_ITER): + # 1. Current residual norm is from previous step (or initial) + res_norm = np.linalg.norm(r) + + # 2. Append to history + residual_history.append(res_norm) + + # 3. Check for convergence BEFORE the update + if res_norm / r0_norm < TOL: + return x, k, residual_history + + # 4. If not converged, compute update + Ar = mat_vec_prod(r, **kwargs) + r_dot_Ar = np.dot(r, Ar) + if r_dot_Ar == 0: # Should not happen for SPD if r!=0 + break + alpha = res_norm**2 / r_dot_Ar + + x += alpha * r + r -= alpha * Ar + + # If loop finishes, append final residual + residual_history.append(np.linalg.norm(r)) + return x, MAX_ITER, residual_history + + +def conjugate_gradient( + mat_vec_prod: Callable[[np.ndarray, Any], np.ndarray], + b: np.ndarray, + **kwargs + ) -> Tuple[np.ndarray, int, List[float]]: + """ + Solves Ax=b using Conjugate Gradient, with corrected residual logging. + """ + x = np.zeros_like(b) + r = b.copy() + p = r.copy() + + r0_norm = np.linalg.norm(r) + if r0_norm == 0: + return x, 0, [0.0] + + res_norm_old = r0_norm + residual_history = [] + + for k in range(MAX_ITER): + # 1. Current residual norm is from previous step + # 2. Append to history + residual_history.append(res_norm_old) + + # 3. Check for convergence + if res_norm_old / r0_norm < TOL: + return x, k, residual_history + + # 4. If not converged, compute update + Ap = mat_vec_prod(p, **kwargs) + p_dot_Ap = np.dot(p, Ap) + if p_dot_Ap == 0: + break + alpha = res_norm_old**2 / p_dot_Ap + + x += alpha * p + r -= alpha * Ap + + res_norm_new = np.linalg.norm(r) + beta = res_norm_new**2 / res_norm_old**2 + + p = r + beta * p + res_norm_old = res_norm_new + # If loop finishes, append final residual + residual_history.append(res_norm_old) + return x, MAX_ITER, residual_history diff --git a/project2/main.tex b/project2/main.tex new file mode 100644 index 0000000..fcc0dc7 --- /dev/null +++ b/project2/main.tex @@ -0,0 +1,187 @@ +\documentclass{article} + +% Language setting +\usepackage[english]{babel} + +% Set page size and margins +\usepackage[letterpaper,top=2cm,bottom=2cm,left=3cm,right=3cm,marginparwidth=1.75cm]{geometry} + +% Useful packages +\usepackage{amsmath} +\usepackage{graphicx} +\usepackage{amsfonts} +\usepackage{float} % For better table and figure placement +\usepackage[colorlinks=true, allcolors=blue]{hyperref} +\usepackage{array} % For better table column definitions + +\title{MATH 6601 - Project 2 Report} +\author{Zhe Yuan (yuan.1435)} +\date{Nov 2025} + +\begin{document} +\maketitle + +\section*{Common Iterative Methods for Both Problems} +This project implements five iterative methods to solve large, sparse, symmetric positive-definite (SPD) linear systems of the form $Ax=b$. The implemented methods are: +\begin{itemize} + \item Jacobi method + \item Gauss-Seidel method + \item Successive Over-Relaxation (SOR) method + \item Steepest Descent method + \item Conjugate Gradient method +\end{itemize} +All five methods are implemented in Python. The core logic for each solver is contained in the `iterative\_solvers.py' file. The specific problems are set up and solved in `problem1.py' and `problem2.py', respectively. + +The initial guess was $x^{(0)}=0$, and the stopping condition was $\|r^{(k)}\|_2 / \|r^{(0)}\|_2 < 10^{-10}$. All $A$ were not explicitly constructed. + +\section{Problem 1: Ordinary Differential Equations} +This problem solves the boundary value problem $-u'' = t$ on $(0,1)$ with $u(0)=u(1)=0$. To solve the problem numerically, the interval $(0,1)$ is discretized into $n+1$ subintervals of width $h=1/(n+1)$. The second derivative $u''(t_i)$ is approximated using a central difference scheme: +\[ u''(t_i) = -t_i \approx \frac{u(t_{i-1}) - 2u(t_i) + u(t_{i+1})}{h^2} \quad \implies \quad -u_{i-1} + 2u_i - u_{i+1} = h^2 t_i \] +for the internal grid points $i=1, \dots, n$. Combined with the boundary conditions $u_0=0$ and $u_{n+1}=0$, this forms a linear system $Au=b$, where $u = [u_1, \dots, u_n]^T$. The matrix $A \in \mathbb{R}^{n \times n}$ is a tridiagonal matrix, and the vector $b \in \mathbb{R}^n$ is given by: +\[ +A = \begin{pmatrix} + 2 & -1 & & & \\ + -1 & 2 & -1 & & \\ + & \ddots & \ddots & \ddots & \\ + & & -1 & 2 & -1 \\ + & & & -1 & 2 +\end{pmatrix}, \quad +b = h^2 \begin{pmatrix} t_1 \\ t_2 \\ \vdots \\ t_n \end{pmatrix} +\] + +\subsection*{(a) Exact Solution} +The boundary value problem is given by $-u'' = t$ on the interval $(0, 1)$ with boundary conditions $u(0) = u(1) = 0$. +Integrating twice gives the general solution: +\[ u(t) = -\frac{t^3}{6} - C_1 t - C_2 \] +Applying the boundary conditions $u(0)=0$ and $u(1)=0$, we find the constants $C_2=0$ and $C_1 = -1/6$. Thus, the exact solution is: +\[ u(t) = \frac{t - t^3}{6} \] +The plot of this exact solution is shown in Figure~\ref{fig:p1_exact_solution}. +\begin{figure}[H] + \centering + \includegraphics[width=0.8\textwidth]{problem1_exact_solution.png} + \caption{Plot of the exact solution $u(t) = (t-t^3)/6$.} + \label{fig:p1_exact_solution} +\end{figure} + +\subsection*{(b, c) Iteration Solution} +According to the requirments from 1(b), The linear system was solved for $n=20$ and $n=40$. For SOR, the optimal relaxation parameter $\omega = \frac{2}{1+\sin(\pi h)}$ was used. The number of iterations required for each method to converge is summarized in Table~\ref{tab:p1_iterations}. + +\begin{table}[H] +\centering +\begin{tabular}{|l|r|r|} +\hline +\textbf{Method} & \textbf{Iterations (n=20)} & \textbf{Iterations (n=40)} \\ +\hline +Jacobi & 2031 & 7758 \\ +Gauss-Seidel & 1018 & 3882 \\ +SOR & 94 & 183 \\ +Steepest Descent & 2044 & 7842\\ +Conjugate Gradient & 20 & 40 \\ +\hline +\end{tabular} +\caption{Number of iterations for each method for Problem 1.} +\label{tab:p1_iterations} +\end{table} + +\noindent \textbf{Observations:} +\begin{enumerate} + \item The Conjugate Gradient method is the most efficient, requiring the fewest iterations. SOR is the second-best with a significant speedup over Gauss-Seidel. Gauss-Seidel converges about twice as fast as Jacobi. The Steepest Descent method performs poorly, rivaling Jacobi in its slow convergence. + \item As $n$ doubles, the number of iterations increases for all methods. The iteration count for Jacobi, Gauss-Seidel, and Steepest Descent increases by a factor of roughly 4, while for SOR and CG, the increase is closer to a factor of 2. +\end{enumerate} + +\subsection*{(d) Convergence History} +The convergence history is shown in Figures~\ref{fig:p1_conv_n20} and \ref{fig:p1_conv_n40}. The stationary methods and SD show linear convergence, while CG shows superlinear convergence. + +\begin{figure}[H] + \centering + \includegraphics[width=0.8\textwidth]{problem1_convergence_n20.png} + \caption{Convergence history for Problem 1 with n=20.} + \label{fig:p1_conv_n20} +\end{figure} + +\begin{figure}[H] + \centering + \includegraphics[width=0.8\textwidth]{problem1_convergence_n40.png} + \caption{Convergence history for Problem 1 with n=40.} + \label{fig:p1_conv_n40} +\end{figure} + +\newpage +\section{Problem 2: Partial Differential Equation} +This problem involves solving the PDE $-\nabla^2 u + u = 1$ (this should be $-\nabla^2 u + u = 1$ instead of $-\nabla u(\text{vector}) + u(\text{scalar}) = 1$ in question) on the unit square $[0,1]^2$ with zero boundary conditions. To solve the problem numerically, a finite difference approximation on an $n \times n$ grid of interior points leads to the discrete equations: +\[ (4+h^2)u_{i,j} - u_{i-1,j} - u_{i,j-1} - u_{i+1,j} - u_{i,j+1} = h^2, \quad i,j=1,\dots,n \] +To formulate this as a standard matrix system $Ax=b$, we reshape the 2D array of unknowns $u_{i,j}$ into a single 1D vector $x$ of length $N=n^2$. After flattening, $u_{i,j}$ maps to the index $k = (i-1)n + j$ in the vector $x$, the matrix $A \in \mathbb{R}^{N \times N}$ becomes a large, sparse, block tridiagonal matrix: +\[ +A = \begin{pmatrix} + T & -I & & & \\ + -I & T & -I & & \\ + & \ddots & \ddots & \ddots & \\ + & & -I & T & -I \\ + & & & -I & T +\end{pmatrix}_{n^2 \times n^2} +\] +where $I$ is the $n \times n$ identity matrix, and $T$ is an $n \times n$ tridiagonal matrix given by: +\[ +T = \begin{pmatrix} + 4+h^2 & -1 & & & \\ + -1 & 4+h^2 & -1 & & \\ + & \ddots & \ddots & \ddots & \\ + & & -1 & 4+h^2 & -1 \\ + & & & -1 & 4+h^2 +\end{pmatrix}_{n \times n} +\] +The right-hand side vector $b \in \mathbb{R}^{N}$ is a constant vector where every element is $h^2$. So: + +\[ +(A_{i,i}=T)_{j,j}x_{(i-1)n + j} = (4+h^2)u_{i,j}, (A_{i-1,i}=I)_{j,j}x_{(i-2)n + j} = u_{i-1,j},\dots +\] + +\subsection*{(a) Iterative Solution} +The system was solved for $n=16$ ($N=256$) and $n=32$ ($N=1024$) using the same five iterative methods. The number of iterations to reach convergence is reported in Table~\ref{tab:p2_iterations}. + +\begin{table}[H] +\centering +\begin{tabular}{|l|r|r|} +\hline +\textbf{Method} & \textbf{Iterations (n=16)} & \textbf{Iterations (n=32)} \\ +\hline +Jacobi & 1268 & 4792 \\ +Gauss-Seidel & 635 & 2397 \\ +SOR & 71 & 138 \\ +Steepest Descent & 1247 & 4807 \\ +Conjugate Gradient & 31 & 66 \\ +\hline +\end{tabular} +\caption{Number of iterations for each method for Problem 2.} +\label{tab:p2_iterations} +\end{table} + +\noindent \textbf{Observations:} +The relative performance of the methods remains consistent with Problem 1: CG is the fastest, followed by SOR, Gauss-Seidel, Jacobi, and Steepest Descent. Doubling $n$ from 16 to 32 quadruples the number of unknowns. Again, the iteration counts for slower methods increase by a factor of approximately 4, while for CG and SOR, the increase is only about a factor of 2. This reinforces the conclusion that CG and SOR scale much better with problem size. + +\subsection*{(b) Convergence History} +The convergence history plots are shown in Figures~\ref{fig:p2_conv_n16} and \ref{fig:p2_conv_n32}. + +\begin{figure}[H] + \centering + \includegraphics[width=0.8\textwidth]{problem2_convergence_n16.png} + \caption{Convergence history for Problem 2 with n=16.} + \label{fig:p2_conv_n16} +\end{figure} + +\begin{figure}[H] + \centering + \includegraphics[width=0.8\textwidth]{problem2_convergence_n32.png} + \caption{Convergence history for Problem 2 with n=32.} + \label{fig:p2_conv_n32} +\end{figure} + +\subsection*{(c) Advantages of Iterative Methods} +For large, sparse linear systems arising from PDE discretizations, iterative methods offer significant advantages over direct methods: +\begin{enumerate} + \item Iterative methods typically only require storing a few vectors of size $N$, leading to a memory complexity of $O(N)$. In contrast, direct methods require denser helper matrices like $Q$ and $R$. + \item The cost per iteration for these problems is dominated by the sparse matrix-vector product, which is $O(N)$. If the method converges in $K$ iterations, the total cost is $O(KN)$. If a good iterative method is used, $K$ can be much smaller than $N$ (like C.G.). Direct sparse solvers often have a higher complexity. When $N$ is very large, iterative methods can be much faster. + \item Implementing iterative methods is often simpler than implementing a direct solver for sparse matrix. +\end{enumerate} + +\end{document} diff --git a/project2/main.tex:Zone.Identifier b/project2/main.tex:Zone.Identifier new file mode 100644 index 0000000..d6c1ec6 Binary files /dev/null and b/project2/main.tex:Zone.Identifier differ diff --git a/project2/problem1.py b/project2/problem1.py new file mode 100644 index 0000000..3413fce --- /dev/null +++ b/project2/problem1.py @@ -0,0 +1,160 @@ +import numpy as np +import matplotlib.pyplot as plt +from iterative_solvers import * + +# (a): Exact Solution +def exact_solution(t: np.ndarray) -> np.ndarray: + """Computes the exact solution u(t) = t/6 * (1 - t^2).""" + return t / 6.0 * (1 - t**2) + +def plot_exact_solution(): + """Plots the exact solution of the BVP.""" + t_fine = np.linspace(0, 1, 500) + u_exact_fine = exact_solution(t_fine) + + plt.figure(figsize=(10, 6)) + plt.plot(t_fine, u_exact_fine, 'k-', label='Exact Solution $u(t) = t/6 \cdot (1 - t^2)$') + plt.title('Exact Solution of $-u\'\'=t$ with $u(0)=u(1)=0$') + plt.xlabel('$t$') + plt.ylabel('$u(t)$') + plt.grid(True) + plt.legend() + plt.savefig('problem1_exact_solution.png') + # plt.show() + print("Exact solution plot saved to 'problem1_exact_solution.png'") + + +# (b): Helper functions for iterative methods +def mat_vec_prod_p1(x: np.ndarray, n: int) -> np.ndarray: + """Computes the matrix-vector product Ax for Problem 1.""" + y = np.zeros(n) + # A is tridiagonal with [ -1, 2, -1 ] + # y_i = -x_{i-1} + 2x_i - x_{i+1} + for i in range(n): + term_prev = -x[i-1] if i > 0 else 0 + term_next = -x[i+1] if i < n - 1 else 0 + y[i] = term_prev + 2 * x[i] + term_next + return y + +def get_diag_inv_p1(i: int, n: int) -> float: + """Returns 1/A_ii for Problem 1. Diagonal elements are all 2, so inverse is 1/2.""" + return 1.0 / 2.0 + +def get_off_diag_sum_p1(x: np.ndarray, i: int, n: int) -> float: + """Computes sum_{j!=i} A_ij * x_j. For Problem 1, A_ij = -1 for j=i-1, i+1, else 0.""" + sum_val = 0.0 + if i > 0: + sum_val += -x[i-1] + if i < n - 1: + sum_val += -x[i+1] + return sum_val + +def solve_and_report(n: int): + """ + Solves the linear system for a given n and reports results. + """ + print(f"\n--- Solving for n = {n} ---") + h = 1.0 / (n + 1) + + # Setup the right-hand side vector b + t = np.linspace(h, 1.0 - h, n) + b = (h**2) * t + + # Calculate SOR parameter omega + omega = 2.0 / (1.0 + np.sin(np.pi * h)) + + # Store results + iteration_counts = [] + all_residuals = {} + all_solutions = {} + # Define methods and run solvers + methods_to_run = { + 'Jacobi': ('jacobi', {}), + 'Gauss-Seidel': ('gauss_seidel', {}), + 'SOR': ('sor', {'omega': omega}), + 'Steepest Descent': (steepest_descent, {}), + 'Conjugate Gradient': (conjugate_gradient, {}) + } + + print("\nStarting iterative solvers...") + for name, (method_func, params) in methods_to_run.items(): + print(f"Running {name}...") + common_args = {'b': b, 'n': n} + if isinstance(method_func, str): + solver_args = { 'method': method_func, 'get_diag_inv': get_diag_inv_p1, + 'get_off_diag_sum': get_off_diag_sum_p1, 'mat_vec_prod': mat_vec_prod_p1, + **common_args, **params } + solution, iters, res_hist = stationary_method(**solver_args) + else: + solver_args = {'mat_vec_prod': mat_vec_prod_p1, **common_args} + solution, iters, res_hist = method_func(**solver_args) + iteration_counts.append((name, iters)) + all_residuals[name] = res_hist + all_solutions[name] = solution if iters < MAX_ITER else None + # (c): Report iteration counts + print("\n--- Iteration Counts ---") + for name, iters in iteration_counts: + status = "converged" if iters < MAX_ITER else "did NOT converge" + print(f"{name:<20}: {iters} iterations ({status}, last residual = {all_residuals[name][-1]:.2e})") + # (d): Plot residue history + plt.figure(figsize=(12, 8)) + for name, res_hist in all_residuals.items(): + res_log = np.log10(np.array(res_hist) + 1e-20) + plt.plot(res_log, label=name, linestyle='-') + + plt.title(f'Convergence History for n = {n}') + plt.xlabel('Iteration Step (m)') + plt.ylabel('log10(||r^(m)||_2)') + plt.legend() + plt.grid(True) + + # Set the y-axis limit to make the plot readable + current_ylim = plt.ylim() + plt.ylim(bottom=max(current_ylim[0], -13), top=current_ylim[1]) + plt.savefig(f'problem1_convergence_n{n}.png') + # plt.show() + print(f"\nConvergence plot saved to 'problem1_convergence_n{n}.png'") + + if all_solutions: + plt.figure(figsize=(12, 8)) + + # Plot exact solution on a fine grid + t_fine = np.linspace(0, 1, 500) + u_exact_fine = exact_solution(t_fine) + plt.plot(t_fine, u_exact_fine, 'k-', label='Exact Solution', linewidth=1, zorder=10) # Black, thick, on top + + # Define some styles for the different numerical methods + styles = { + 'Jacobi': {'color': 'red', 'linestyle': '--'}, + 'Gauss-Seidel': {'color': 'blue', 'linestyle': '-.'}, + 'SOR': {'color': 'green', 'linestyle': ':'}, + 'Steepest Descent':{'color': 'purple','linestyle': '--'}, + 'Conjugate Gradient':{'color': 'orange','linestyle': '-.'}, + } + + # Plot each numerical solution + for name, solution in all_solutions.items(): + # Add boundary points for a complete plot + t_numerical_with_bounds = np.concatenate(([0], t, [1])) + u_numerical_with_bounds = np.concatenate(([0], solution, [0])) + + style = styles.get(name, {'color': 'gray', 'linestyle': '-'}) + plt.plot(t_numerical_with_bounds, u_numerical_with_bounds, + label=f'Numerical ({name})', **style, linewidth=2) + + plt.title(f'Comparison of Exact and Numerical Solutions for n = {n}') + plt.xlabel('$t$') + plt.ylabel('$u(t)$') + plt.legend() + plt.grid(True) + plt.savefig(f'problem1_all_solutions_comparison_n{n}.png') + # plt.show() + print(f"All solutions comparison plot saved to 'problem1_all_solutions_comparison_n{n}.png'") + +if __name__ == "__main__": + # Part (a) + plot_exact_solution() + + # Part (b), (c), (d) + solve_and_report(n=20) + solve_and_report(n=40) diff --git a/project2/problem1_all_solutions_comparison_n20.png b/project2/problem1_all_solutions_comparison_n20.png new file mode 100644 index 0000000..3ef0f11 Binary files /dev/null and b/project2/problem1_all_solutions_comparison_n20.png differ diff --git a/project2/problem1_all_solutions_comparison_n40.png b/project2/problem1_all_solutions_comparison_n40.png new file mode 100644 index 0000000..6911456 Binary files /dev/null and b/project2/problem1_all_solutions_comparison_n40.png differ diff --git a/project2/problem1_convergence_n20.png b/project2/problem1_convergence_n20.png new file mode 100644 index 0000000..a4bdeaf Binary files /dev/null and b/project2/problem1_convergence_n20.png differ diff --git a/project2/problem1_convergence_n40.png b/project2/problem1_convergence_n40.png new file mode 100644 index 0000000..afa2478 Binary files /dev/null and b/project2/problem1_convergence_n40.png differ diff --git a/project2/problem1_exact_solution.png b/project2/problem1_exact_solution.png new file mode 100644 index 0000000..8a0da2e Binary files /dev/null and b/project2/problem1_exact_solution.png differ diff --git a/project2/problem1_output.txt b/project2/problem1_output.txt new file mode 100644 index 0000000..e4f2a2f --- /dev/null +++ b/project2/problem1_output.txt @@ -0,0 +1,39 @@ +Exact solution plot saved to 'problem1_exact_solution.png' + +--- Solving for n = 20 --- + +Starting iterative solvers... +Running Jacobi... +Running Gauss-Seidel... +Running SOR... +Running Steepest Descent... +Running Conjugate Gradient... + +--- Iteration Counts --- +Jacobi : 2031 iterations (converged, last residual = 5.78e-13) +Gauss-Seidel : 1018 iterations (converged, last residual = 5.72e-13) +SOR : 94 iterations (converged, last residual = 4.52e-13) +Steepest Descent : 2044 iterations (converged, last residual = 5.74e-13) +Conjugate Gradient : 20 iterations (converged, last residual = 1.39e-18) + +Convergence plot saved to 'problem1_convergence_n20.png' +All solutions comparison plot saved to 'problem1_all_solutions_comparison_n20.png' + +--- Solving for n = 40 --- + +Starting iterative solvers... +Running Jacobi... +Running Gauss-Seidel... +Running SOR... +Running Steepest Descent... +Running Conjugate Gradient... + +--- Iteration Counts --- +Jacobi : 7758 iterations (converged, last residual = 2.15e-13) +Gauss-Seidel : 3882 iterations (converged, last residual = 2.16e-13) +SOR : 183 iterations (converged, last residual = 2.04e-13) +Steepest Descent : 7842 iterations (converged, last residual = 2.16e-13) +Conjugate Gradient : 40 iterations (converged, last residual = 5.40e-19) + +Convergence plot saved to 'problem1_convergence_n40.png' +All solutions comparison plot saved to 'problem1_all_solutions_comparison_n40.png' diff --git a/project2/problem2.py b/project2/problem2.py new file mode 100644 index 0000000..e64a8ce --- /dev/null +++ b/project2/problem2.py @@ -0,0 +1,147 @@ +import numpy as np +import matplotlib.pyplot as plt +from iterative_solvers import * + +def mat_vec_prod_p2(x_flat: np.ndarray, n: int, h_sq: float) -> np.ndarray: + """ + Computes the matrix-vector product y = Ax for the 2D PDE problem. + """ + x_2d = x_flat.reshape((n, n)) + x_padded = np.zeros((n + 2, n + 2)) + x_padded[1:-1, 1:-1] = x_2d + + y_2d = np.zeros((n, n)) + + # Apply the finite difference stencil: (4+h^2)u_ij - u_{i-1,j} - u_{i+1,j} - u_{i,j-1} - u_{i,j+1} + for i in range(n): + for j in range(n): + y_2d[i, j] = (4.0 + h_sq) * x_padded[i + 1, j + 1] \ + - x_padded[i, j + 1] \ + - x_padded[i + 2, j + 1] \ + - x_padded[i + 1, j] \ + - x_padded[i + 1, j + 2] + + return y_2d.flatten() + +def get_diag_inv_p2(i: int, n: int, h_sq: float) -> float: + """Returns 1/A_ii for Problem 2. Diagonal is constant: 4 + h^2.""" + return 1.0 / (4.0 + h_sq) + +def get_off_diag_sum_p2(x_flat: np.ndarray, i_flat: int, n: int, **kwargs) -> float: + """Computes the off-diagonal sum for a given row `i_flat` in Problem 2.""" + row = i_flat // n + col = i_flat % n + sum_val = 0.0 + + if row > 0: sum_val += -x_flat[i_flat - n] + if row < n - 1: sum_val += -x_flat[i_flat + n] + if col > 0: sum_val += -x_flat[i_flat - 1] + if col < n - 1: sum_val += -x_flat[i_flat + 1] + + return sum_val + +def solve_and_report_p2(n: int): + """Solves the 2D PDE linear system for a given n and reports results.""" + print(f"\n--- Solving for n = {n} (System size: {n*n} x {n*n}) ---") + h = 1.0 / (n + 1) + h_sq = h**2 + N = n * n + + b = np.full(N, h_sq) + omega = 2.0 / (1.0 + np.sin(np.pi * h)) + + iteration_counts = [] + all_residuals = {} + all_solutions = {} + + methods_to_run = { + 'Jacobi': ('jacobi', {}), + 'Gauss-Seidel': ('gauss_seidel', {}), + 'SOR': ('sor', {'omega': omega}), + 'Steepest Descent': (steepest_descent, {}), + 'Conjugate Gradient': (conjugate_gradient, {}) + } + + print("\nStarting iterative solvers...") + for name, (method_func, params) in methods_to_run.items(): + print(f"Running {name}...") + + solver_kwargs = {'n': n, 'h_sq': h_sq} + + if isinstance(method_func, str): + solver_args = {'b': b, 'method': method_func, 'get_diag_inv': get_diag_inv_p2, + 'get_off_diag_sum': get_off_diag_sum_p2, 'mat_vec_prod': mat_vec_prod_p2, + **solver_kwargs, **params} + solution, iters, res_hist = stationary_method(**solver_args) + else: + solver_args = {'b': b, 'mat_vec_prod': mat_vec_prod_p2, **solver_kwargs} + solution, iters, res_hist = method_func(**solver_args) + + iteration_counts.append((name, iters)) + all_residuals[name] = res_hist + all_solutions[name] = solution if iters < MAX_ITER else None + + # Report iteration counts + print("\n--- Iteration Counts ---") + for name, iters in iteration_counts: + status = "converged" if iters < MAX_ITER else "did NOT converge" + print(f"{name:<20}: {iters} iterations ({status})") + + # Plot residue history + plt.figure(figsize=(12, 8)) + for name, res_hist in all_residuals.items(): + res_log = np.log10(np.array(res_hist) + 1e-20) + plt.plot(res_log, label=name, linestyle='-') + + plt.title(f'Convergence History for PDE Problem, n = {n}') + plt.xlabel('Iteration Step (m)') + plt.ylabel('log10(||r^(m)||_2)') + plt.legend() + plt.grid(True) + + current_ylim = plt.ylim() + plt.ylim(bottom=max(current_ylim[0], -13), top=current_ylim[1]) + + plt.savefig(f'problem2_convergence_n{n}.png') + print(f"\nConvergence plot saved to 'problem2_convergence_n{n}.png'") + + if all_solutions: + print("Generating 3D surface plots of the solutions...") + + # Create the X, Y grid for the surface plot, including boundaries + plot_grid = np.linspace(0, 1, n + 2) + X, Y = np.meshgrid(plot_grid, plot_grid) + + # Create a figure for the subplots + fig = plt.figure(figsize=(18, 10)) + fig.suptitle(f'3D Surface Plots of Numerical Solutions for n = {n}', fontsize=16) + # Loop through the stored solutions and create a subplot for each + for i, (name, solution_flat) in enumerate(all_solutions.items()): + # Reshape the 1D solution vector into a 2D grid + solution_2d = solution_flat.reshape((n, n)) + + # Pad the solution with zero boundaries for plotting + u_padded = np.zeros((n + 2, n + 2)) + u_padded[1:-1, 1:-1] = solution_2d + + # Add a new 3D subplot + ax = fig.add_subplot(2, 3, i + 1, projection='3d') + + # Plot the surface + ax.plot_surface(X, Y, u_padded, cmap='viridis', edgecolor='none') + + # Set titles and labels for each subplot + ax.set_title(name) + ax.set_xlabel('x') + ax.set_ylabel('y') + ax.set_zlabel('u(x,y)') + ax.view_init(elev=30, azim=-60) # Set a nice viewing angle + plt.tight_layout(rect=[0, 0.03, 1, 0.95]) # Adjust layout to make room for suptitle + plt.savefig(f'problem2_3d_solutions_n{n}.png') + print(f"3D surface comparison plot saved to 'problem2_3d_solutions_n{n}.png'") + +if __name__ == "__main__": + # Solve Problem 2 for n=16 and n=32 + solve_and_report_p2(n=16) + solve_and_report_p2(n=32) + diff --git a/project2/problem2_3d_solutions_n16.png b/project2/problem2_3d_solutions_n16.png new file mode 100644 index 0000000..f155323 Binary files /dev/null and b/project2/problem2_3d_solutions_n16.png differ diff --git a/project2/problem2_3d_solutions_n32.png b/project2/problem2_3d_solutions_n32.png new file mode 100644 index 0000000..6d78d3b Binary files /dev/null and b/project2/problem2_3d_solutions_n32.png differ diff --git a/project2/problem2_convergence_n16.png b/project2/problem2_convergence_n16.png new file mode 100644 index 0000000..1d437e5 Binary files /dev/null and b/project2/problem2_convergence_n16.png differ diff --git a/project2/problem2_convergence_n32.png b/project2/problem2_convergence_n32.png new file mode 100644 index 0000000..abdeecf Binary files /dev/null and b/project2/problem2_convergence_n32.png differ diff --git a/project2/problem2_output.txt b/project2/problem2_output.txt new file mode 100644 index 0000000..3cf02af --- /dev/null +++ b/project2/problem2_output.txt @@ -0,0 +1,40 @@ + +--- Solving for n = 16 (System size: 256 x 256) --- + +Starting iterative solvers... +Running Jacobi... +Running Gauss-Seidel... +Running SOR... +Running Steepest Descent... +Running Conjugate Gradient... + +--- Iteration Counts --- +Jacobi : 1268 iterations (converged) +Gauss-Seidel : 635 iterations (converged) +SOR : 71 iterations (converged) +Steepest Descent : 1247 iterations (converged) +Conjugate Gradient : 31 iterations (converged) + +Convergence plot saved to 'problem2_convergence_n16.png' +Generating 3D surface plots of the solutions... +3D surface comparison plot saved to 'problem2_3d_solutions_n16.png' + +--- Solving for n = 32 (System size: 1024 x 1024) --- + +Starting iterative solvers... +Running Jacobi... +Running Gauss-Seidel... +Running SOR... +Running Steepest Descent... +Running Conjugate Gradient... + +--- Iteration Counts --- +Jacobi : 4792 iterations (converged) +Gauss-Seidel : 2397 iterations (converged) +SOR : 138 iterations (converged) +Steepest Descent : 4807 iterations (converged) +Conjugate Gradient : 66 iterations (converged) + +Convergence plot saved to 'problem2_convergence_n32.png' +Generating 3D surface plots of the solutions... +3D surface comparison plot saved to 'problem2_3d_solutions_n32.png'