\documentclass{report}
\pagestyle{empty}

\usepackage{amsfonts}


\renewcommand{\dot}{{\cdot}}
\renewcommand{\ddot}{{:}}
\newcommand{\cross}{{\times}}
\newcommand{\Grad}{\mathrm{grad}}
\newcommand{\Div}{\mathrm{div}}
\newcommand{\Curl}{\mathrm{curl}}
\newcommand{\grad}{{\nabla}}
\renewcommand{\div}{{\nabla\cdot}}
\newcommand{\curl}{{\nabla\times}}
\newcommand{\Laplacian}{{\Delta}}
\newcommand{\laplacian}{{\Delta}}


\begin{document}

\noindent
A. J. Meir\hfill\break

\centerline{\bf Math-5620/6620}
\centerline{Mathematical Computation and Scientific Visualization}
\centerline{Fall 2006}
\vskip5pt
%\centerline {{\Huge DRAFT} \qquad Final Project \qquad {\Huge DRAFT}}
\centerline {Final Project}

\vskip40pt
The project is due on or before Dec. 13.

You are to submit a complete project write-up (short report) which you are proud of....
The report should include a brief introduction, description of the problem, methods you
use for solution (use to construct approximate solutions), and discussion of
parallelizability, scalability, speedup, and parallel efficiency, if applicable.

You are to write a computer program to approximate solutions of a heat equation posed
on a 2-dimensional spatial domain. A copy of the program, as well as, visualizations of
solutions of representative examples, should be submitted.

The problem we consider is a heat equation posed on a 2-dimensional spatial domain
(on the unit square).

\[
  \frac{\partial u}{\partial t}  -\div(K \grad u ) = f \qquad \mathrm{in\ }  Q = \Omega\times (0, T) \,.
\]
Here \( \Omega = (0,1)\times (0,1) \) and \( T \) is some given final time, \( u \) is the temperature,
\( K \) is the heat conduction coefficient (this can be a tensor, e.g. ), and \( f \) is a heat source, or sink.
You may consider problems with a constant heat conduction coefficient \( K \), ones
where this coefficient is position dependent \( K = K(x) \), as well as one where this
is also temperature dependent \( K = K(x, u) \).

We denote the boundary of \( \Omega \) by \( \partial\Omega \). In order to describe the possible
boundary conditions, assume \( \partial\Omega_i \) for \( i=1\ldots 3 \) are open, mutually
disjoint subsets of \( \partial\Omega \) and that
\( \partial\Omega = \cup_{i=1}^3 \overline{\partial\Omega_i} \).
Also denote by \( \Gamma_i \) part of the lateral boundary \( \partial\Omega_i\times (0, T) \).

The partial differential equation is supplemented by an initial condition
\[
  u(x,0) = u_0(x) \qquad \mathrm{in\ } \Omega
\]
and boundary conditions
\[
  u(x,t)|_{\Gamma_1} = g(x,t)
\]
\[
  -K\grad u(x,t)\cdot\mathbf{n}|_{\Gamma_2} = h(x,t)
\]

\[
  -K\grad u(x,t)\cdot\mathbf{n}|_{\Gamma_3} = H{u(x,t)|_{\Gamma_3} - u_\mathrm{a}(x,t)}
\]
where \( u_\mathrm{a} \) is the ambient temperature (you may assume it is constant, if you
choose to use this condition), \( H \) is a given constant and \( g \) and \( h \) are given functions.

Note, you may be able to concisely write the three possible boundary conditions as
\[
  \alpha(x) K\grad u(x,t)\cdot\mathbf{n}|_{\Gamma} + \beta(x) u(x,t)|_{\Gamma} = \gamma(x,t)
\]
for an appropriate choice of \( \alpha \), \( \beta \), and \( \gamma \).

We can discretize the p.d.e. as follows, for \( 0 \le \theta \le 1 \)
\begin{eqnarray}
  \frac{u^{n+1}_{i,j} - u^n_{i,j}}{k} & - & \theta \Bigg [  
  {K_1}_{i,j}\frac{u^{n+1}_{i+1,j} - 2 u^{n+1}_{i,j} + u^{n+1}_{i-1, j}}{h_1^2}
  + \frac{{K_1}_{i+1,j} - {K_1}_{i-1,j}}{2h_1}\frac{u^{n+1}_{i+1,j} - u^{n+1}_{i-1,j}}{2 h_1} \nonumber \\
  & + &  {K_2}_{i,j}\frac{u^{n+1}_{i,j+1} - 2 u^{n+1}_{i,j} + u^{n+1}_{i, j-1}}{h_2^2}
  + \frac{{K_2}_{i,j+1} - {K_2}_{i,j-1}}{2h_2}\frac{u^{n+1}_{i,j+1} - u^{n+1}_{i,j-1}}{2 h_2} \Bigg ] \nonumber \\
  & - & (1 - \theta) \Bigg [  
  {K_1}_{i,j}\frac{u^n_{i+1,j} - 2 u^n_{i,j} + u^n_{i-1, j}}{h_1^2}
  + \frac{{K_1}_{i+1,j} - {K_1}_{i-1,j}}{2h_1}\frac{u^n_{i+1,j} - u^n_{i-1,j}}{2 h_1} \nonumber \\
  & + &  {K_2}_{i,j}\frac{u^n_{i,j+1} - 2 u^n_{i,j} + u^n_{i, j-1}}{h_2^2}
  + \frac{{K_2}_{i,j+1} - {K_2}_{i,j-1}}{2h_2}\frac{u^n_{i,j+1} - u^n_{i,j-1}}{2 h_2} \Bigg ] \nonumber \\
  = \theta f^{n+1}_{i,j} & + & (1 - \theta) f^{n}_{i,j}\nonumber \,.
\end{eqnarray}

For \( \theta = 0 \) the scheme is explicit, for \( \theta = 1\) it is fully implicit, for \( \theta = 1/2 \) it is
known as the Crank-Nicholson scheme. The explicit scheme is stable (and hence, since
consistent also convergent) if
\( kK_1/h_1^2 + kK_2/h_2^2 \le 1/2 \). When \( 0 \le \theta < 1/2 \) the scheme is stable if and
only is \( kK_1/h_1^2 + kK_2/h_2^2 \le 1/(2 - 4\theta) \), when \( 1/2 \le \theta < 1 \) the scheme
is stable for all \( kK_1/h_1^2 + kK_2/h_2^2 \).


The boundary conditions may be discretized as follows. First order discretization
\[
  \frac{\partial u}{\partial x_1}(x_{i,j}, t_n) \approx \frac{u^n_{i+1,j} - u^n_{i,j}}{h_1}
\]
one sided second order discretization
\[
  \frac{\partial u}{\partial x_1}(x_{i,j}, t_n) \approx \frac{-u^n_{i+2,j} +4 u^n_{i+1,j} - 3u^n_{i, j}}{2h_1}
\]
or a centered second order discretization (requires the introduction of ghost points)
\[
  \frac{\partial u}{\partial x_1}(x_{i,j}, t_n) \approx \frac{u^n_{i+1,j} - u^n_{i-1, j}}{2h_1} \,.
\]

For simplicity you may want to consider only the constant coefficient case, with Dirichlet boundary
conditions. Also for simplicity you may want to limit attention to the explicit case. This can be solved
without constructing matrices, and is more easily parallelizable (though in this case there is a
constraint on the time step that can be used, see stability condition, above).
\end{document}