188 lines
9.0 KiB
TeX
188 lines
9.0 KiB
TeX
\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}
|