%%%%		-*- Mode: LaTeX -*-
%%%% Author:    J. Robert Buchanan
%%%% Address:   Department of Mathematics
%%%%		Millersville University
%%%%            P.O. Box 1002
%%%%            Millersville, PA 17551-0302
%%%% Phone:     717-871-7305
%%%% FAX:	717-871-7948
%%%%
%%%% Email:     Robert.Buchanan@millersville.edu
%%%%
%%%% Date:      04/05/2025
%%%%
%%%% Purpose:   Iterative Methods
%%%%

\documentclass[10pt]{beamer}
\usepackage{mathtools}

\mode<presentation>
{
  %\usetheme{Warsaw}
  \usetheme{default}
  %\usetheme{CambridgeUS}
  % or ...

  %\setbeamercovered{transparent}
  % or whatever (possibly just delete it)
}


\usepackage[english]{babel}
% or whatever

\usepackage[latin1]{inputenc}
% or whatever

\usepackage{times}
\usepackage[T1]{fontenc}
% Or whatever. Note that the encoding and the font should match. If T1
% does not look nice, try deleting the line with the fontenc.


\title{Iterative Methods}
\subtitle{\textit{Partial Differential Equations}} % (optional)


\author % (optional, use only with lots of authors)
{J~Robert~Buchanan}
% - Use the \inst{?} command only if the authors have different
%   affiliation.

\institute[Millersville University of Pennsylvania]
% (optional, but mostly needed)
{
  Department of Mathematics
}
% - Use the \inst command only if there are several affiliations.
% - Keep it simple, no one is interested in your street address.

\date{ % (optional)
Fall
%Spring
%Summer
\the\year}

\subject{Talks}
% This is only inserted into the PDF information catalog. Can be left
% out. 



% If you have a file called "university-logo-filename.xxx", where xxx
% is a graphic format that can be processed by latex or pdflatex,
% resp., then you can add a logo as follows:

% \pgfdeclareimage[height=0.5cm]{university-logo}{university-logo-filename}
% \logo{\pgfuseimage{university-logo}}



% Delete this, if you do not want the table of contents to pop up at
% the beginning of each subsection:
\AtBeginSubsection[]
{
  \begin{frame}<beamer>
    \frametitle{Outline}
    \tableofcontents[currentsection,currentsubsection]
  \end{frame}
}


% If you wish to uncover everything in a step-wise fashion, uncomment
% the following command: 

%\beamerdefaultoverlayspecification{<+->}


\begin{document}

\begin{frame}
  \titlepage
\end{frame}

%\begin{frame}
%  \frametitle{Outline}
%  \tableofcontents
  % You might wish to add the option [pausesections]
%\end{frame}


% Since this a solution template for a generic talk, very little can
% be said about how it should be structured. However, the talk length
% of between 15min and 45min and the theme suggest that you stick to
% the following rules:  

% - Exactly two or three sections (other than the summary).
% - At *most* three subsections per section.
% - Talk about 30s to 2min per frame. So there should be between about
%   15 and 30 frames, all told.

%\section{Introduction}

\begin{frame}\frametitle{Objectives}

In this lesson we will:
\begin{itemize}
\item present and use the Jacobi method for approximating the solution
  to a linear system of equations,
\item present and use the Gauss-Seidel method for approximating the solution
  to a linear system of equations,
\item present and use the Successive Over-Relaxation (SOR) method for
  approximating the solution to a linear system of equations,
\end{itemize}

\end{frame}
%%

\begin{frame}\frametitle{Background}

If $A$ is an $n \times n$ matrix and $\mathbf{u}$ and $\mathbf{b}$ are
column vectors with $n$ components, then let $\mathbf{u}^{(0)}$ be an
approximation to $\mathbf{u}$ where $\mathbf{u}$ is the solution of
\[
A \, \mathbf{u} = \mathbf{b}	.
\]
Each of the three iterative methods generates a sequence of vectors $\{
\mathbf{u}^{(k)} \}_{k=0}^{\infty}$ which converges to $\mathbf{u}$.

\end{frame}
%%

\begin{frame}\frametitle{Jacobi Method}

Consider a linear system of equations expressed in matrix/vector form
as $A \, \mathbf{u} = \mathbf{b}$.
Suppose the diagonal entries of matrix $A$ are all nonzero.
The $i$th equation of this system can be expressed as
\[
a_{i1} u_{1} + a_{i2} u_{2} + \cdots + a_{ii} u_{i} + \cdots + a_{in}
u_{n}
	= b_{i}	.
\]
Solving this equation for $u_{i}$ yields
\[
u_{i}
	= \frac{1}{a_{ii}} \left( b_{i} - \sum_{\mathclap{j=1, j \ne i}}^{n} a_{ij} u_{j} \right)	.
\]
Thus if $\mathbf{u}^{(k)} = (u_{1}^{(k)}, u_{2}^{(k)}, \ldots,
u_{n}^{(k)})^{T}$ for some $k = 0, 1, \ldots$ is an approximation to
the solution of $A \, \mathbf{u} = \mathbf{b}$ then the next approximation is
$\mathbf{u}^{(k+1)}$ where
\[
%\label{eq:findiff:Jacobi}
u_{i}^{(k+1)}
	= \frac{1}{a_{ii}} \left( b_{i} - \sum_{\mathclap{j=1, j
      \ne i}}^{n} a_{ij} u_{j}^{(k)} \right)	,
\]
for $i = 1, 2, \ldots, n$.

\end{frame}
%%

\begin{frame}\frametitle{Jacobi Method}

If matrix $A = D + L + U$ where $D$ is a
diagonal matrix who entries are the diagonal entries of $A$, $L$ is a
strictly lower triangular matrix (zeros on the diagonal) whose entries
are the entries of $A$ strictly below the diagonal, and $U$ is a
strictly upper triangular matrix whose entries are the entries of $A$
strictly above the diagonal, then
\[
A \, \mathbf{u}
	= (D + L + U) \mathbf{u}
	= \mathbf{b}	.
\]
If the diagonal entries of $A$ are all nonzero then $D$ is invertible
and the Jacobi method may be written in matrix/vector form as
\[
%\label{eq:findiff:Jacobimatrix}
\mathbf{u}^{(k+1)}
	= -D^{-1} (L+U) \mathbf{u}^{(k)} + D^{-1} \mathbf{b}
	= T_{J} \mathbf{u}^{(k)} + D^{-1} \mathbf{b}
\]
where $T_{J} = -D^{-1} (L+U)$ is called the
\textbf{Jacobi iteration matrix}.

\end{frame}
%%

\begin{frame}\frametitle{Stopping Condition}

A stopping condition should indicate that further iteration is
unnecessary (either because the solution to the linear system has been
found or because the results of further iteration will not appreciably
change).
Many choices of stopping condition are available and a generally
effective choice is that the iteration should cease when the relative
change in successive approximations falls below a specified threshold.
Thus if the threshold is given by $\epsilon > 0$, calculation of
further iterates of the sequence $\{ \mathbf{u}^{(k)} \}$ can stop
when
\[
%\label{eq:findiff:Jacobistop}
\frac{\| \mathbf{u}^{(k)} -\mathbf{u}^{(k+1)} \|_{2}}{\|
  \mathbf{u}^{(k+1)} \|_{2}}
	< \epsilon	.
\]
The subscript of $2$ indicates the $l_{2}$-norm is used to measure the
vectors.

\end{frame}
%%

\begin{frame}\frametitle{Example}

%\label{xmp:findiff:Jacobi}
Use a maximum of five iterations of the Jacobi method to approximate a
solution to the following system of linear equations.
Let $\epsilon = 10^{-4}$ and $\mathbf{u}^{(0)} = (1, 1, 1, 1)^{T}$.
\begin{align*}
4 u_{1} + u_{2} - u_{3} + u_{4}
	&= 6	\\
u_{1} + 4 u_{2} - u_{3} - u_{4}
	&= 6	\\
-u_{1} -  u_{2} + 4 u_{3} + u_{4}
	&= 3	\\
u_{1} -  u_{2} + u_{3} + 3 u_{4}
	&= -2
\end{align*}
For the purposes of comparison the exact solution is $\mathbf{u} =
\frac{1}{37} (80, 38, 73, -63)^{T}$.

\end{frame}
%%

\begin{frame}\frametitle{Solution}

\begin{center}
\begin{tabular}{c|rrrr}
$k$ & $u_{1}^{(k)}$ & $u_{2}^{(k)}$ & $u_{3}^{(k)}$ & $u_{4}^{(k)}$
  \\
\hline
0 & $1.0000$ & $1.0000$ & $1.0000$ & $1.0000$ \\
1 & $1.2500$ & $1.7500$ & $1.0000$ & $-1.0000$ \\
2 & $1.5625$ & $1.1875$ & $1.7500$ & $-0.8333$ \\
3 & $1.8490$ & $1.3385$ & $1.6458$ & $-1.3750$ \\
4 & $1.9206$ & $1.1055$ & $1.8906$ & $-1.3854$ \\
5 & $2.0426$ & $1.1462$ & $1.8529$ & $-1.5686$
\end{tabular}
\end{center}
The $l_{2}$-norm error in $\mathbf{u}^{(5)}$ is
\[
\| \mathbf{u}^{(5)} - \mathbf{u} \|_{2} \approx 0.2468
\]
The stopping condition was not achieved in the first five
iterations, thus emphasizing the need for practical implementations of
the Jacobi method (and other iterative methods) to specify a maximum
number of iterations to compute.

\end{frame}
%%

\begin{frame}\frametitle{Gauss-Seidel Method}

The Gauss-Seidel method modifies the Jacobi method by using the
already calculated entries of $\mathbf{u}^{(k+1)}$ in place of the
corresponding entries of $\mathbf{u}^{(k)}$.
\[
%\label{eq:findiff:GS}
u_{i}^{(k+1)}
	= \frac{1}{a_{ii}} \left( b_{i} -
\sum_{j=1}^{i-1} a_{ij} u_{j}^{(k+1)} -
\sum_{j=i+1}^{n} a_{ij} u_{j}^{(k)}
\right)
\]

\end{frame}
%%

\begin{frame}\frametitle{Gauss-Seidel Method}

If $A = D + L + U$ then
\begin{align*}
A \, \mathbf{u}
	&= \mathbf{b}	\\
(D + L + U) \mathbf{u}
	&= \mathbf{b}	\\
\mathbf{u}
	&= D^{-1} (\mathbf{b} - L \, \mathbf{u} - U \, \mathbf{u}).
\end{align*}
If $a_{ii} \ne 0$ for all $i = 1, 2, \ldots, n$ then the inverse of
diagonal matrix $D$ exists.
The Gauss-Seidel formula can then be expressed
as
\[
\mathbf{u}^{(k+1)}
	= D^{-1} (\mathbf{b} - L \, \mathbf{u}^{(k+1)} - U
\, \mathbf{u}^{(k)})	.
\]
Collecting the $k+1$st iterates on the left-hand side of the equation
enables the Gauss-Seidel iteration formula to be written as
\[
%\label{eq:findiff:GSmatrix}
\mathbf{u}^{(k+1)}
	= -(D+L)^{-1} U \mathbf{u}^{(k)} + (D+L)^{-1} \mathbf{b}
	= T_{GS} \mathbf{u}^{(k)} + (D+L)^{-1} \mathbf{b}
\]
where $T_{GS} = -(D+L)^{-1} U$ is called the
\textbf{Gauss-Seidel iteration matrix}.

\end{frame}
%%

\begin{frame}\frametitle{Example}

%\label{xmp:findiff:GS}
Use the Gauss-Seidel method to approximate a
solution to the following system of linear equations.
Let $\epsilon = 10^{-4}$ and $\mathbf{u}^{(0)} = (1, 1, 1, 1)^{T}$.
\begin{align*}
4 u_{1} + u_{2} - u_{3} + u_{4}
	&= 6	\\
u_{1} + 4 u_{2} - u_{3} - u_{4}
	&= 6	\\
-u_{1} -  u_{2} + 4 u_{3} + u_{4}
	&= 3	\\
u_{1} -  u_{2} + u_{3} + 3 u_{4}
	&= -2
\end{align*}
For the purposes of comparison the exact solution is $\mathbf{u} =
\frac{1}{37} (80, 38, 73, -63)^{T}$.

\end{frame}
%%

\begin{frame}\frametitle{Solution}

\begin{center}
\begin{tabular}{c|rrrr}
$k$ & $u_{1}^{(k)}$ & $u_{2}^{(k)}$ & $u_{3}^{(k)}$ & $u_{4}^{(k)}$
  \\
\hline
0 & $1.0000$ & $1.0000$ & $1.0000$ & $1.0000$ \\
1 & $1.2500$ & $1.6875$ & $1.2344$ & $-0.9323$	\\
2 & $1.6198$ & $1.1706$ & $1.6807$ & $-1.3766$	\\
3 & $1.9717$ & $1.0831$ & $1.8578$ & $-1.5821$	\\
4 & $2.0892$ & $1.0466$ & $1.9295$ & $-1.6574$	\\
5 & $2.1351$ & $1.0343$ & $1.9567$ & $-1.6858$
\end{tabular}
\end{center}
The error in the last approximation to the solution is
\[
\| \mathbf{u}^{(5)} - \mathbf{u} \|_{2} \approx 0.0002
\]
which is clearly an improvement over the Jacobi method's result after
five iterations.

\end{frame}
%%

\begin{frame}\frametitle{Example}

Consider the following boundary value problem:
\begin{align*}
u_{xx} + u_{yy}
	&= 0
\mbox{ for $0 < x < 1$ and $0 < y < 1$}	\\
u(x, 0)
	&= e^{x}
\mbox{ for $0 < x < 1$}	\\
u(x, 1)
	&= \sin (\pi x/2)
\mbox{ for $0 < x < 1$}	\\
u(0, y)
	&= 1 - y^{2}
\mbox{ for $0 < y < 1$}	\\
u(1, y)
	&= e^{1-y}
\mbox{ for $0 < y < 1$.}
\end{align*}
If $h = 1/10$, then the linear system $A \, \mathbf{u} = \mathbf{b}$
consists of a $81 \times 81$ matrix $A$ and a column vector
$\mathbf{b}$ with 81 elements.
Setting a numerical tolerance of $\epsilon = 10^{-6}$
and using the Gauss-Seidel method
requires on the order of 100 iterations from a starting approximation
of $\mathbf{u}^{(0)} = (1, 1, \ldots, 1)^{T}$.

\end{frame}
%%

\begin{frame}\frametitle{Solution}

\begin{center}
\resizebox{0.75\textwidth}{!}{\includegraphics{exampleLaplace}}
\end{center}

\end{frame}
%%

\begin{frame}\frametitle{Successive Over-Relaxation Method}

Since $\mathbf{u}^{(k+1)} - \mathbf{u}^{(k)}$ can be thought of as the
step to take from approximation $\mathbf{u}^{(k)}$ to the better
approximation $\mathbf{u}^{(k+1)}$, a natural improvement to the
method is to take a longer step in the same direction.

\vspace*{0.10in}

If $\omega > 1$ then $\omega (\mathbf{u}^{(k+1)} - \mathbf{u}^{(k)})$
is a vector in the same direction as $\mathbf{u}^{(k+1)} -
\mathbf{u}^{(k)}$ but with greater magnitude.
By taking this longer corrective step an algorithm can ``accelerate''
toward the exact solution.
\[
%\label{eq:findiff:relaxation}
\mathbf{u}^{(k+1)}
	= \mathbf{u}^{(k)} + \omega (\mathbf{u}^{(k+1)} - \mathbf{u}^{(k)})	.
\]
When $\omega > 1$ this iterative method is called an
\textbf{over-relaxation method}.
When $0 < \omega < 1$ this iterative method is called an
\textbf{under-relaxation method}.

\vspace*{0.10in}

The \textbf{successive over-relaxation
method} (abbreviated as SOR)
is an accelerated modification of the Gauss-Seidel method.

\end{frame}
%%

\begin{frame}\frametitle{Matrix/Vector Form}

\begin{align*}
\mathbf{u}^{(k+1)}
	&= \mathbf{u}^{(k)} + \omega \left[
D^{-1} (\mathbf{b} - L \, \mathbf{u}^{(k+1)} - D \, \mathbf{u}^{(k)} -
   U \, \mathbf{u}^{(k)}) \right]	\\
%\label{eq:findiff:SORmatrix}
(D + \omega L) \mathbf{u}^{(k+1)}
	&= \left[ (1 - \omega) D - \omega U \right] \mathbf{u}^{(k)} +
\omega \mathbf{b}	\\
\mathbf{u}^{(k+1)}
	&= T_{SOR} \mathbf{u}^{(k)} +
\omega (D + \omega L)^{-1} \mathbf{b}
\end{align*}
where $T_{SOR} = (D + \omega L)^{-1} \left[ (1 - \omega) D - \omega
  U \right]$ is called the
\textbf{SOR iteration matrix}.

\vspace*{0.10in}

The $i$th component of $\mathbf{u}^{(k+1)}$ can be written as
\[
%\label{eq:findiff:SOR}
u_{i}^{(k+1)}
	= u_{i}^{(k)} + \frac{\omega}{a_{ii}} \left(
b_{i} - \sum_{j=1}^{i-1} a_{ij} u_{j}^{(k+1)} - \sum_{j=i}^{n} a_{ij} u_{j}^{(k)}
\right)	.
\]
Note that if $\omega = 1$ then the successive over-relaxation method
is just the Gauss-Seidel method.

\end{frame}
%%

\begin{frame}\frametitle{Example}

Compare the results of the Gauss-Seidel and SOR methods with $\omega =
1.21539$ for approximating the solution to the following linear system of
equations.
For both methods use $\epsilon = 10^{-4}$ in the stopping criterion.
\[
\left[ \begin{array}{rrrrr}
-4 & 1 & 0 & 0 & 0	\\
1 & -4 & 1 & 0 & 0	\\
0 & 1 & -4 & 1 & 0	\\
0 & 0 & 1 & -4 & 1	\\
0 & 0 & 0 & 1 & -4
\end{array} \right]
\left[ \begin{array}{c}
x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5}
\end{array} \right]
	= \left[ \begin{array}{c}
1 \\ 2 \\ 3 \\ 4 \\ 5
\end{array} \right]
\]

\end{frame}
%%

\begin{frame}\frametitle{Gauss-Seidel Results}

The initial approximation to the solution will be
$(0, 0, 0, 0, 0)^{T}$.
The first ten iterations of the Gauss-Seidel method yield a close
approximation of the solution to the equation.
The exact solution is $\mathbf{x} =
(-129/260, -64/65, -75/52, -116/65, -441/260)^{T}$.
\begin{center}
\begin{tabular}{c|rrrrr}
$k$ & $x_{1}^{(k)}$ & $x_{2}^{(k)}$ & $x_{3}^{(k)}$ & $x_{4}^{(k)}$ & $x_{5}^{(k)}$	\\
\hline
0 & $0.0000$ & $0.0000$ & $0.0000$ & $0.0000$ & $0.0000$	\\
1 & $-0.2500$ & $-0.5625$ & $-0.8906$ & $-1.2227$ & $-1.5557$	\\
2 & $-0.3906$ & $-0.8203$ & $-1.2607$ & $-1.7041$ & $-1.6760$	\\
3 & $-0.4551$ & $-0.9290$ & $-1.4083$ & $-1.7711$ & $-1.6928$	\\ 
4 & $-0.4822$ & $-0.9726$ & $-1.4359$ & $-1.7822$ & $-1.6955$	\\
5 & $-0.4932$ & $-0.9823$ & $-1.4411$ & $-1.7842$ & $-1.6960$	\\
6 & $-0.4956$ & $-0.9842$ & $-1.4421$ & $-1.7845$ & $-1.6961$	\\
7 & $-0.4960$ & $-0.9845$ & $-1.4423$ & $-1.7846$ & $-1.6961$	\\
8 & $-0.4961$ & $-0.9846$ & $-1.4423$ & $-1.7846$ & $-1.6962$	\\
9 & $-0.4961$ & $-0.9846$ & $-1.4423$ & $-1.7846$ & $-1.6962$	\\
10& $-0.4962$ & $-0.9846$ & $-1.4423$ & $-1.7846$ & $-1.6962$	\\
\hline
\end{tabular}
\end{center}

\end{frame}
%%

\begin{frame}\frametitle{SOR Results}

The SOR method has similar performance though to the four decimal
places reported in the table, the 9th and 10th iterates are the same.
\begin{center}
\begin{tabular}{c|rrrrr}
$k$ & $x_{1}^{(k)}$ & $x_{2}^{(k)}$ & $x_{3}^{(k)}$ & $x_{4}^{(k)}$ & $x_{5}^{(k)}$	\\
\hline
0 & $0.0000$ & $0.0000$ & $0.0000$ & $0.0000$ & $0.0000$	\\
1 & $-0.3038$ & $-0.7000$ & $-1.1242$ & $-1.5570$ & $-1.9923$	\\
2 & $-0.4511$ & $-0.9356$ & $-1.4268$ & $-1.9189$ & $-1.6732$	\\
3 & $-0.4910$ & $-0.9889$ & $-1.4878$ & $-1.7625$ & $-1.6944$	\\
4 & $-0.4986$ & $-0.9982$ & $-1.4299$ & $-1.7851$ & $-1.6967$	\\
5 & $-0.4998$ & $-0.9790$ & $-1.4434$ & $-1.7850$ & $-1.6962$	\\
6 & $-0.4937$ & $-0.9854$ & $-1.4424$ & $-1.7846$ & $-1.6961$	\\
7 & $-0.4969$ & $-0.9847$ & $-1.4423$ & $-1.7846$ & $-1.6962$	\\
8 & $-0.4960$ & $-0.9845$ & $-1.4423$ & $-1.7846$ & $-1.6962$	\\
9 & $-0.4962$ & $-0.9846$ & $-1.4423$ & $-1.7846$ & $-1.6962$	\\
10 & $-0.4962$ & $-0.9846$ & $-1.4423$ & $-1.7846$ & $-1.6962$	\\
\hline
\end{tabular}
\end{center}

\end{frame}
%%

\begin{frame}\frametitle{Remarks}

\begin{itemize}
\item There is no known formula for the optimal value of $\omega$ to
accelerate the SOR method for a general $n \times n$ matrix $A$.
\item For certain banded matrices (for example tridiagonal matrices) the
optimal value of $\omega$ can be determined \textit{a priori}.
\item For example, if $\omega \approx 1.52786$ then the SOR method
approximates the solution to a boundary value problem earlier in
the previous section in 25 iterations as opposed to nearly 100
iterations required for the Gauss-Seidel method.
\end{itemize}

\end{frame}
%%

\begin{frame}\frametitle{Example}

Consider the boundary value problem:
\begin{align*}
\Delta u
	&= e^{x - y}
\mbox{ on $R$}	\\
u(x, y)
	&= x \, y
\mbox{ on $\partial R$.}
\end{align*}
Let the domain of the boundary value problem be $R = \{ (x, y) \, | \,
0 < x < 4, \, 0 < y < 3 \}$.

\vspace*{0.10in}

Use the successive over-relaxation method with $\omega = 1.81621$, the
nine-point approximation to the Laplacian operator,
and $h = 1/10$ to approximate the solution to the boundary value problem.
Use $\epsilon = 10^{-6}$ in the stopping criterion.

\end{frame}
%%

\begin{frame}\frametitle{Solution (1 of 2)}

The rectangle over which the solution is sought has dimensions of
width $4$ and height $3$.
With Dirichlet boundary conditions and choosing $h = 1/10$ there are
\[
\left( \frac{4}{1/10} - 1 \right) \left( \frac{3}{1/10} - 1 \right)
	= 1131
\]
interior points at which the solution to Poisson's equation must be
approximated.
The linear system to be solved is of the form $A \, \mathbf{u} =
\mathbf{b}$ where $A$ is a banded tridiagonal matrix of size $1131
\times 1131$.
Despite the size of the linear system, the SOR technique solves it
in 92 iterations from a
starting approximation of $\mathbf{u}^{(0)} = \mathbf{0}$.

\end{frame}
%%

\begin{frame}\frametitle{Solution (2 of 2)}

\begin{center}
\resizebox{0.75\textwidth}{!}{\includegraphics{9ptPoissonSOR}}
\end{center}

\end{frame}
%%

%\begin{frame}\frametitle{Homework}

%\begin{itemize}
%\item Read Section $x.y$
%\item Exercises: 
%\end{itemize}

%\end{frame}
%%

%%%
\end{document}

%%%%
%%%% EOF
%%%%
