\documentclass{article}
 
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{float} % an [H] forces a float (such as algorithm) to appear HERE!
\usepackage{graphics}
\usepackage{listings}
\usepackage{tcolorbox}
\usepackage{verbatim}
 
\marginparwidth 0pt
\oddsidemargin  0pt
\evensidemargin  0pt
\marginparsep 0pt
\topmargin   0pt
\textwidth   6.5in
\textheight  8.5in
\setlength{\parindent}{0em}

\newenvironment{packed_itemize}{
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{itemize}}

\newenvironment{packed_enumerate}{
\begin{enumerate}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{enumerate}}

\definecolor{codeback}{RGB}{235,255,255}

\lstset{ %
language=Matlab,                % choose the language of the code
basicstyle=\footnotesize,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\footnotesize,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=10pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{codeback},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,           % adds a frame around the code
tabsize=2,          % sets default tabsize to 2 spaces
captionpos=b,           % sets the caption-position to bottom
breaklines=true,        % sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={(*@}{@*)}          % if you want to add a comment within your code
}

\begin{document}

\begin{center}
  \Huge{The steady 2D Poisson equation}\\
  \Large{MATH1091: ODE methods for a reaction diffusion equation}
  \normalsize{http://people.sc.fsu.edu/$\sim$jburkardt/classes/math1091\_2020/poisson\_2d\_steady/poisson\_2d\_steady.pdf}
\end{center}

\begin{center}
  \includegraphics[scale=0.40]{electric_field.png}
  \vskip 0.1in
  {\it{Electric field potential is an example of a Poisson problem.}}
\end{center}

\vskip 0.1in
\begin{tcolorbox}[colback=red!5,title=The steady Poisson Problem]
{\it{Given a Poisson equation on a 2D rectangular region, use
finite differences to create a model of the equation, set
up the corresponding linear system, plot the approximate
solution and compute the norm of the approximation error.
}}
\end{tcolorbox}

\section{Moving to 2D spatial domains}

For our initial value problems and boundary value problems, we have been looking at a
one dimensional domain:
\begin{packed_itemize}
  \item{ a time interval $\Omega = t_{init} \leq t \leq t_{final}$ for IVP's;}
  \item{ a space interval $\Omega = x_{left} \leq x \leq x_{right}$ for BVP's;}
\end{packed_itemize}
We'd like to look at some problems on more interesting domains.  In this discussion,
we will move to a 2D spatial domain.  Later, we will consider problems involving
variation over a domain that extends over both space and time.  To handle such problems,
we will generalize the methods that we developed for the 1D steady (time independen)
Poisson problem.
\vskip 0.1in
To keep things simple, our 2D domain $\Omega$ will usually be a rectangle.  This simplifies settiing up our
discretized geometry.  Although we don't have time to cover them, there are plenty of techniques 
for handling more complicated regions, such as a circle, a slice of pizza, or the outline of 
an arbitrary curve such as the border of France.

\section{The 2D Poisson equation}

The steady (time-independent) Poisson equation in $\Omega \subset \mathbb{R}^2$ can be written as:
\begin{displaymath}
  - \Delta u = f(x,y) \quad \text{ for all } x, y \in \Omega
\end{displaymath}
We will assume that $\Omega$ is a rectangle, $\Omega = \{(x,y): a \leq x \leq b, c \leq y \leq d \}$.
Here, in 2D, the mysterious symbol $\Delta u$ represents
\begin{displaymath}
  \Delta u \equiv \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = u_{xx} + u_yy
\end{displaymath}
In order to be able to determine how well our approximations are doing, we will always
assume that we know a formula $g(x,y)$ for the exact solution.  For simplicity, we will
generally assume that the boundary conditions will be of Dirichlet type, and that means
we can use the $g(x,y)$ function to define these boundary conditions:
\begin{align*}
  u(x,y) = g(x,y) \text{ for all } x,y \in \partial \Omega.
\end{align*}
Here, $\partial \Omega$ is the set of all boundary points of $\Omega$.  For a rectangle, 
$a \leq x \leq b, c \leq y \leq d$, this is simply the points on the perimeter, that is
\begin{displaymath}
\partial \Omega = \{ (x,y) \in \Omega: x=a {\text{ or }} x=b {\text{ or }} y=c {\text{ or }} y=d \}
\end{displaymath}
\vskip 0.1in
For this lecture, we will stick with this simple version of the Poisson equation.
Later we may see variations of the Poisson equation, involving new boundary conditions,
the inclusion of spatially dependent coefficients, or extra terms in the equation, involving $u$,
or $\frac{\partial u}{\partial x}$ or $\frac{\partial u}{\partial y}$.

\section{A Poisson equation on a 2D rectangle}

As we did for the 1D equation, we will replace the mathematical problem by a discrete
computational problem.  We start by making a discrete model of the the geometry.
We use a regular mesh of points; those along
the boundary will be associated with the boundary conditions, while points in the interior
will be associated with the Poisson equation.  At each point $(i,j)$, we will write
an equation involving the approximation $U_{i,j}$.  These equations will form a linear
system to be solved.  Our solution, then, will be a table of these values.
\vskip 0.1in
Because our problem is 2D, we have to be more careful about how we generate the mesh points,
how we form the linear system, and how we try to report the results.  
\vskip 0.1in
In particular, we will now be working with a mix of vectors (lists) and matrices (tables).
A useful MATLAB convention suggests using capital letters for matrices.  
\vskip 0.1in
To discretize the geometry, we might use a grid of {\tt{nx}} points in
the $x$ direction and {\tt{ny}} points in $y$.  In MATLAB, this can be done by:
\begin{lstlisting}
  x = linspace ( a, b, nx );
  y = linspace ( c, d, ny );
\end{lstlisting}
so that the node in column $i$ and row $j$ of the mesh has coordinates {\tt{(x(i),y(j))}}
\vskip 0.1in
We will seek an approximate solution only at this discrete set of points,
and we may use the notation $U_{i,j}$ to indicate the value of our
discrete solution at $(x_i,y_j)$.
\vskip 0.1in
Points corresponding to $i=1$, $i=nx$, $j=1$ or $j=ny$ are {\it{boundary
points}} where the Dirichlet boundary conditions will be imposed,
using the function {\tt{g(x,y)}}:
\begin{displaymath}
  U_{i,j} = g(x_i,y_j)
\end{displaymath}
At the remaining {\it{interior points}} we impose a
discretized version of the Poisson equation.  Recall that for the 1D equation,
we had a way to approximate the second derivative $u''$.  Now that $u$ is
a function of $x$ and $y$, we simply use the same approximate technique for
$\frac{\partial^2 u}{\partial x^2}$ and $\frac{\partial^2 u}{\partial y^2}$,
so that our discretized version of the Poisson equation at the point $(x_i,y_j)$ is:
\begin{displaymath}
  \frac{-U_{i-1,j}+2U_{i,j}-U_{i+1,j}}{hx^2} + 
  \frac{-U_{i,j-1}+2U_{i,j}-U_{i,j+1}}{hy^2} = f(x_i,y_j)
\end{displaymath}
\vskip 0.1in
This equation involves the solution values at 5 nodes, which we can think of as
{\tt{N=north}}, {\tt{S=south}}, {\tt{E=east}}, {\tt{W=west}} and {\tt{C=center}}.  
We think of this as equation number $k$ in our linear system, with coefficient matrix $A$.  
Then we can identify each coefficient $A(*,*)$ as occuring in matrix row $k$, the column 
corresponding to each of the five variables.  
\begin{verbatim}
A(k,W)*U(W) + A(k,E)*U(E) + A(k,N)*U(N) + A(k,S)*U(S) + A(k,C)*U(C) = rks(k)
\end{verbatim}
The values of the coefficients and right hand side will be:
\begin{verbatim}
  A(k,W) = -1/hx^2
  A(k,E) = -1/hx^2
  A(k,N) = -1/hy^2
  A(k,S) = -1/hy^2
  A(k,C) = 2/hx^2 + 2/hy^2
  rhs(k) = f(x(i),y(j) )
\end{verbatim}
Now we have a system of $nx \times ny$ equations, which we symbolize by
\begin{displaymath}
  A * u = rhs
\end{displaymath}
We expect our solution {\tt{U}} to be a matrix, but the MATLAB linear equation
solver is going to return {\tt{u}} as a vector, so that's another issue we are going to
have to deal with.  If we handle things correctly, then the command that turns that
vector into a table will be something like this:
\begin{lstlisting}
  u = A \ rhs;
  U = reshape ( u, nx, ny );
\end{lstlisting}
but we will explain what is going on a little later!

\section{Exercise \#1: Create a 4x3 discretized system}

Consider a Poisson problem defined on $\Omega$, the rectangle $[0,3] \times [0,2]$,
with $f(x,y)=-2*y$, exact solution $g(x,y)=x^2y$, and Dirichlet boundary conditions.
\vskip 0.1in
If we use a $4 \times 3$ mesh of points, we get 12 nodes, and 12 equations, 10 of
which are boundary equations, and 2 are discretized Poisson equations.  Suppose
we number the nodes (and the corresponding equations) as follows:
\begin{verbatim}
  2 | 9 10 11 12
Y=1 | 5  6  7  8
  0 | 1  2  3  4
    +-----------
      0  1  2  3
           X
\end{verbatim}
{\bf{Write the 12 linear equations of the discretized system.}}

\section{Assembling the linear system}

We have $nx \times ny$ equations to construct, and we need a way
to order them, just as I did in the diagram above.  We will use $i$ and $j$ 
to count $x$ and $y$ coordinates, and we want $k$ to count nodes and equations.  
A procedure to visit each node in order might look something like this:
\begin{lstlisting}[caption = Logic needed to visit each node.]
  k = 0;
  for j = 1 : ny
    for i = 1 : nx
      k = k + 1;
      if ( boundary)
        equation k is boundary equation at xi yj
      else
        equation k is Poisson equation at xi yj
      end
    end
  end
\end{lstlisting}
Assuming we have defined {\tt{function value = g(x,y)}}, a
boundary equation can be thought of as
\begin{lstlisting}
  U_{i,j} = g(x(i),y(j))
\end{lstlisting}
while, using {\tt{hx}} and {\tt{hy}} as the $x$ and $y$ spacings, and
{\tt{function value = f(x,y)}} for the right hand side, we have seen that 
the Poisson equation at node in column $i$ and row $j$ can be thought of as 
\begin{lstlisting}
   (-U(i-1,j) +2*U(i,j)-U(i+1,j))/hx^2 
  +(-U(i,j-1) +2*U(i,j)-U(i,j+1))/hy^2 = f(x(i),y(j))
\end{lstlisting}

\section{Relating column and row indices to nodes or (x,y) coordinates}

If we think of our nodes as a table or matrix, it is natural to
index them using $i$ and $j$.  However, in order to work with the linear
solver, we also need to index the nodes by a single index variable $k$, that
counts from 1 to $nx*ny$. 
\vskip 0.1in
The illustration below suggests an $nx=6 \times ny=5$ rectangle, for which 
we can identify a node by specifying {\tt{(i,j)}} coordinates, or a node index
which we might call {\tt{k}}.  In particular, when writing the discretized
Poisson equation at a node with coordinates {\tt{(i,j)}} and index {\tt{k}},
it is natural to describe the 4 neighbors as N=north at {\tt{(i,j+1)}}, S=south at {\tt{(i,j-1)}},
E=east at {\tt{(i+1,j)}}, W=west at {\tt{(i-1,j)}}.  However, our data is stored by node
index, so while we think of the North value in terms of its $i$ and $j$ indices,
we need to be able to produce the corresponding node index $k$.
\begin{verbatim}
     ^
   5 |   25 26 27 28 29 30
   4 |   19 20 21 22 23 24
   3 J   13 14 15 16 17 18
   2 |    7  8  9 10 11 12
   1 |    1  2  3  4  5  6
     |
     +-----------I-------->
          1  2  3  4  5  6
\end{verbatim}
Looking at the illustration, we can see that if we let the center variable be in 
column i=5, row j=3, and node index k=17, then its North node is in column i=5, row j=4,
and we need a way to automatically compute that the North node index is k=23.
\vskip 0.1in
The following one line function will take any pair $(i,j)$ in a grid and return
the corresponding $k$ index:
\begin{lstlisting}[caption = k\_from\_ij.m computes the node index.]
function k = k_from_ij ( i, j, nx )
  k = i + ( j - 1 ) * nx;
  return
end
\end{lstlisting}
In other words, since nx=6, the node with i=5, j=4 has k=23 because $5 + ( 4 - 1 ) * 6 = 23$.
\vskip 0.1in
Although we won't need it today, a similar function can be written to return the 
$(i,j)$ coordinates given the $k$ index:
\begin{lstlisting}[caption = ij\_from\_k.m computes the (ij) indices.]
function [ i, j ] = ij_from_k ( k, nx )
  i = mod ( k, nx );
  j = floor ( k / nx ) + 1;
  return
end
\end{lstlisting}

\section{Exercise \#2: Set up the linear system for the 4x3 problem}

Here is a sketch of a MATLAB function which sets up and solves the linear system 
for the 4x3 Poisson problem that you looked at in Exercise \#1.  
\vskip 0.1in
In order to create an equation for every node, we consider every combination of $i$ and $j$
using a pair of nested loops.  When we are ready to create a new equation, we increase the
$k$ index by 1.  The value $k$ identifies both the node we are considering, and the corresponding
equation.
\vskip 0.1in
We set the equation by defining the coefficients of node $k$, and its four neighbors, which
are the nodes to the east, west, north and south of it.  Since our center node $k$ has mesh
indices $(i,j)$, it's easy to generate the corresponding $i$ and $j$ indices of the four neighbors.
However, in order to put the coefficients into the matrix, we need to identify those nodes 
by their corresponding node index.  That's where we can use the {\tt{k\_from\_ij()}}function.
Here is a skeleton code that you need to fill in:
\begin{lstlisting}[caption = skeleton2.m]
  a = ?
  b = ?
  c = ?
  d = ?
  nx = ?
  ny = ?
  hx = ?
  x = ?
  hy = ?
  y = ?

  A = zeros(?,?);  
  rhs = zeros(?,1);

  k = 0;
  for j = 1 : ny
    for i = 1 : nx
      k = k + 1;
      if ( boundary )  % what values of i and j mean we are on the boundary?
        A(k,k) = ?
        rhs(k) = ?
      else
        CENTER = k;
        A(k,CENTER) = ?
        WEST = ?      % Work out the (i,j) coordinates of WEST, and use k_from_ij().
        A(k,WEST) = ?
        EAST = ?
        A(k,EAST) = ?
        SOUTH = ?
        A(k,SOUTH) = ?
        NORTH = ?
        A(k,NORTH) = ?
        rhs(k) = ?
      end
    end
  end
  
  u = A \ rhs

define function g(x,y) here
define function f(x,y) here
\end{lstlisting}
{\bf{Starting from {\it{skeleteon2.m}}, create a MATLAB function {\it{exercise2.m}} 
that sets up and solves the tiny Poisson linear system.}}
\vskip 0.1in
Your {\tt{u}} solution should be {\tt{[0,0,0,0,0,1,4,9,0,2,8,18]}}.

\section{Exercise \#3: Postprocessing by tabulating}

For small problems, we can report the results by simply printing the solution
vector.  However, it's helpful to include extra information, such as the
mesh indices and the spatial coordinates.
\vskip 0.1in
Here is one way to do this:
\begin{lstlisting}[caption=Tabulate solution]
  u = A \ rhs;
%
%  Tabulate the solution.
%
  fprintf ( 1, ' k   i   j    x       y       u       exact\n' );
  k = 0;  
  for j = 1 : ny
    fprintf ( 1, '\n' );
    for i = 1 : nx
      k = k + 1;
      fprintf ( 1, '%2d  %2d  %2d  %6.2f  %6.2f  %6.2f  %6.2f\n', ...
                      k,   i,   j, x(i),  y(j),  u(k),  g(x(i),y(j)) );
    end
  end
\end{lstlisting}
{\bf{Starting from {\it{exercise2.m}}, create a MATLAB function {\it{exercise3.m}} 
which tabulates your solution.}}

\section{Exercise \#4: Use a general solver}

Suppose we want to solve a different Poisson problem?  Since the program we wrote is fairly
complicated, is it possible to identify a part that doesn't have to change?  If so, we can
put that as a separate function, {\it{poisson\_solver.m}}, and just feed it the new information.
\vskip 0.1in
Here is how it might work, allowing the user to specify the size and discretization of
the rectangle, and the two special functions $f(x,y)$ and $g(x,y)$.  For this example,
we repeat the previous calculation.
\begin{lstlisting}[caption=exercise4 calls poisson\_solver.m]
a = 0.0;
b = 3.0;
c = 0.0;
d = 2.0;
nx = 4;
ny = 3;
%
%  Notice that data is now returned as TABLES (matrices).
%
[ X, Y, U ] = poisson_solver ( a, b, c, d, nx, ny, @f, @g );
%
%  You could do some post processing here using tables 
%  k, i, j, X(i,j), Y(i,j), U(i,j), g(X(i,j),Y(i,j)).
%

function value = f ( x, y )
  value = - 2.0 * y;
  return
end

function value = g ( x, y )
  value = x.^2 .* y;
  return
end
\end{lstlisting}
{\bf{Create this program {\it{exercise4.m}}, using a downloaded copy of {\it{poisson\_solver.m}}, 
and report the results in a table.}}

\section{Postprocessing: The RMS error norm}

Printing the data in a table is OK for small problems, but for large problems we want
something easier to understand quickly.  If we know the exact solution $g(x,y)$, which
we usually will for these sample problems, then we may be interested in reporting the
error in our approximation.  
\vskip 0.1in
We assume that {\tt{poisson\_solver()}} has returned the results {\tt{X, Y, U}} as matrices.
We can create a corresponding table of the exact solution at all the nodes by $V=g(X,Y)$.
(This assumes that the $g(x,y)$ function has been written using MATLAB's ``dot'' operators
for multiplication, division, and exponentiation!).  Then the error is the table
defined by {\tt{E=U-V}}.  To use MATLAB's {\tt{rms()}} function, we need to temporarily
make {\tt{E}} look like a vector, and so the error norm is calculated this way:
\begin{lstlisting}
  enorm = rms ( reshape ( E, nx*ny, 1 ) );
\end{lstlisting}
Once {\tt{enorm}} has been computed, it should be reported as the RMS norm of the error:
\begin{lstlisting}
  fprintf ( 1, '  nx = %d, ny = %d, rms error norm = %g\n', nx, ny, enorm );
\end{lstlisting}

\section{Exercise \#5: Compute RMS error norm for a new problem}

Our first example problem was too easy to solve.  We actually get zero error.  In order to
study convergence, we want a problem whose solution is harder to approximate.
\vskip 0.1in
Consider the Poisson problem on the rectangle $0 \leq x \leq 1$, $1 \leq y \leq 3$, with
right hand side $f(x,y)=\pi^2*(x^2+y^2)*\sin(\pi*x*y)$, for which the exact solution is
$g(x,y)=\sin(\pi*x*y)$.
\vskip 0.1in
{\bf{Create a MATLAB file {\it{exercise5.m}} which sets up this new problem, solves it
with {\tt{poisson\_solver()}}, and reports the error norm.  Use $nx=11$ and $ny=21$.}}
\vskip 0.1in
If you run your program again, with $nx=21$ and $ny=41$, the mesh spacing $h$ has decreased
by a factor of $\frac{1}{2}$.  By what factor does the error norm decrease?

\section{Homework: Plot the solution of the new problem}

The {\tt{poisson\_solver()}} code returns the solution data as tables {\tt{X, Y, U}}.
\vskip 0.1in
MATLAB uses exactly this type of data to make contour line, contour color, and surface plots.
Here are the commands to make this happen:
\begin{lstlisting}[caption=Contour and surface plots]
figure ( )
contour ( X, Y, U )
title ( 'Contour lines' );
print ( '-dpng', 'hw2_contour_lines.png' );

figure ( )
contourf ( X, Y, U )
title ( 'Contour colors' );
print ( '-dpng', 'hw2_contour_colors.png' );

figure ( )
surf ( X, Y, U );
title ( 'Surface' );
print ( '-dpng', 'hw2.png' );
\end{lstlisting}
{\bf{Starting from {\it{exercise5.m}}, make a MATLAB file {\it{hw2.m}} that creates
these three plots of the solution.}}
\vskip 0.1in
For the color contour and surface plots, you could call {\tt{colorbar()}} to display the
relationship between color and value.
\vskip 0.1in
By default, the surface plot includes black grid lines on the surface.  Especially when a
large number of mesh points are used, this can make the surface hard to see.  An alternative
command then is {\tt{surf ( X, Y, U, 'edgecolor', 'none' )}};
\vskip 0.1in
Notice that you can click in the surface plot window and ``move'' the surface so you can
study it.
\vskip 0.1in
{\bf{Send your surface plot {\it{hw2.png}} to me at {\tt{jvb25@pitt.edu}}.
I would like to see your work by Friday, 15 May 2020.}}

\end{document}
