%%%%		-*- 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/04/2025
%%%%
%%%% Purpose:   Heat Equation
%%%%

\documentclass[10pt]{beamer}

\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{Heat Equation}
\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 develop an explicit method for numerically approximating the
  solution to a heat/diffusion equation initial boundary value
  problem, and
\item develop an implicit method (Crank-Nicolson) for numerically approximating the
  solution to a heat/diffusion equation initial boundary value
  problem.
\end{itemize}

\end{frame}
%%

\begin{frame}\frametitle{Initial Boundary Value Problem}

Numerically approximate the solution to the following Fourier problem:
\begin{align*}
u_{t}
	&= u_{xx}
\mbox{ for $0 < x < 1$ and $t > 0$}	\\
u(0, t)
	&= u(1,t) = 0
\mbox{ for $t > 0$}	\\
u(x, 0)
	&= 100
\mbox{ for $0 < x < 1$.}
\end{align*}
The spatial interval $[0, 1]$ must be subdivided into $N$ intervals of
length $h = 1/N$.
A forward time step of size $k > 0$ must be chosen and a differencing
scheme must be selected for the partial differential equation.
\[
%\label{eq:findiff:explicitheat}
\frac{u(x, t+k) - u(x, t)}{k}
	\approx \frac{u(x+h, t) - 2 u(x, t) + u(x-h, t)}{h^{2}}	.
\]
This finite difference approximation is classified as having truncation error
$O(k) + O(h^{2})$.

\end{frame}
%%

\begin{frame}\frametitle{Explicit Scheme}

\begin{align*}
\frac{1}{k} (u_{i}^{j+1} - u_{i}^{j})
	&= \frac{1}{h^{2}} (u_{i+1}^{j} - 2 u_{i}^{j} + u_{i-1}^{j})	\\
%\label{eq:findiff:explicitr}
u_{i}^{j+1}
	&= r u_{i-1}^{j} + (1-2r) u_{i}^{j} + r u_{i+1}^{j}
\end{align*}
for $i = 1, 2, \ldots, N-1$, where $r = k/h^{2}$.
This is an example of an
\textbf{explicit}
finite difference scheme.

\vspace*{0.10in}

For $i = 0$ and $i = N$ the homogeneous Dirichlet boundary conditions
provide $u_{0}^{j} = u_{N}^{j} = 0$ for all $j \in \mathbb{N}$.

\vspace*{0.10in}

Note that the explicit scheme gives the value of $u_{i}^{j+1}$
as a function of three values of $u$ one time step earlier.

\end{frame}
%%

\begin{frame}\frametitle{Stencil}

\begin{center}
\resizebox{1.00\textwidth}{!}{\includegraphics{explicitheat}}
\end{center}
The values of $u_{i}^{0}$ are the initial conditions and are thus
known.
This enables the solution to be approximated at $t = j \, k$ for $j
\in \mathbb{N}$ iteratively.

\end{frame}
%%

\begin{frame}\frametitle{Matrix/Vector Form}
\fontsize{8pt}{9.6}\selectfont

If the vector $\mathbf{u}^{(j)} = (u_{1}^{j}, u_{2}^{j}, \ldots,
u_{N-1}^{j})^{T}$ then the system of equations
may be written as
$\mathbf{u}^{(j+1)}  = A(r) \, \mathbf{u}^{(j)}$ or
\[
%\label{eq:findiff:explicitrmatrix}
\left[ \begin{array}{c}
u_{1}^{j+1} \\ u_{2}^{j+1} \\ u_{3}^{j+1} \\ \vdots \\ u_{N-3}^{j+1}
         \\ u_{N-2}^{j+1} \\ u_{N-1}^{j+1}
\end{array} \right]
	= \left[ \begin{array}{ccccccc}
1-2r & r & 0 & \cdots & 0 & 0 & 0	\\
r & 1-2r & r & \cdots & 0 & 0 & 0	\\
0 & r & 1-2r & \cdots & 0 & 0 & 0	\\
\vdots & \vdots & \vdots & {} & \vdots & \vdots & \vdots	\\
0 & 0 & 0 & \cdots & 1-2r & r & 0	\\
0 & 0 & 0 & \cdots & r & 1-2r & r	\\
0 & 0 & 0 & \cdots & 0 & r & 1-2r
\end{array} \right]
\left[ \begin{array}{c}
u_{1}^{j} \\ u_{2}^{j} \\ u_{3}^{j} \\ \vdots \\ u_{N-3}^{j} \\ u_{N-2}^{j} \\ u_{N-1}^{j}
\end{array} \right]	.
\]

\end{frame}
%%

\begin{frame}\frametitle{Results}

If $h = 1/10$ and $k = 1/1000$ then $r = 1/10$ and the values
of $u(x, t)$ are shown for $t = 0.012$, $t = 0.036$, and $t = 0.060$.
\begin{center}
\resizebox{0.90\textwidth}{!}{\includegraphics{explicitr10}}
\end{center}
The solid curves provide a comparison with the Fourier series solution.

\end{frame}
%%

\begin{frame}\frametitle{Remarks}

\begin{itemize}
\item Since the forward time step size is $k = 0.001$, one thousand
iterations of the explicit formula are required to approximate the
solution for $0 < t \le 1$.
\item To reduce the number of iterations required to reach the $t = 1$
milestone, $k$ could be increased.
\item If $k = 0.004$ only two hundred fifty iterations are required.
\end{itemize}

\end{frame}
%%

\begin{frame}\frametitle{Results}

If $h = 1/10$ and $k = 0.004$ then $r = 0.4$ and the values
of $u(x, t)$ are shown for $t = 0.012$, $t = 0.036$, and $t = 0.060$.
\begin{center}
\resizebox{0.90\textwidth}{!}{\includegraphics{explicitr2}}
\end{center}
The solid curves provide a comparison with the Fourier series solution.
The numerical solution does not appear to agree as well with the
infinite series solution when $r = 0.4$.

\end{frame}
%%

\begin{frame}\frametitle{Results}

If $h = 1/10$ and $k = 0.006$ then $r = 0.6$ and the values
of $u(x, t)$ are shown for $t = 0.012$, $t = 0.036$, and $t = 0.060$.
\begin{center}
\resizebox{0.90\textwidth}{!}{\includegraphics{explicitr6}}
\end{center}
The solid curves provide a comparison with the Fourier series solution.
The numerical solution does not agree well with the
infinite series solution when $r = 0.6$.

\end{frame}
%%

\begin{frame}\frametitle{Crank-Nicolson Method}

Suppose $f$ is a twice continuously differentiable function on a
neighborhood of $t + k/2$, and expand $f$ as a Taylor polynomial
about $t + k/2$.
\[
f(s)
	= f\left( t + \frac{k}{2} \right) + f'\left( t + \frac{k}{2}
\right) \left(s - t - \frac{k}{2} \right) + O\left( \left(s - t -
    \frac{k}{2} \right)^{2} \right)	,
\]
for $s$ close to $t + k/2$.
Replace $s$ by $t$ and again by $t+k$:
\begin{align*}
f(t)
	&= f\left( t + \frac{k}{2} \right) - f'\left( t + \frac{k}{2}
\right) \left(\frac{k}{2} \right) + O(k^{2})	\\
f(t+k)
	&= f\left( t + \frac{k}{2} \right) + f'\left( t + \frac{k}{2}
\right) \left(\frac{k}{2} \right) + O(k^{2})	.
\end{align*}
Add these two equations and solve for $f(t + k/2)$.
\[
f\left( t + \frac{k}{2} \right) = \frac{f(t) + f(t+k)}{2} + O(k^{2})	,
\]
which has an $O(k^{2})$ truncation error.

\end{frame}
%%

\begin{frame}\frametitle{Developing an Implicit Scheme}

Choose $f$ to be $u_{xx}(x, t)$, then
\begin{align*}
u_{xx}(x, t + k/2)
	= {} & \frac{u_{xx}(x, t) + u_{xx}(x, t+k)}{2} + O(k^{2})	\\
	= {} & \frac{u(x+h, t) - 2 u(x, t) + u(x-h, t)}{2 h^{2}}	\\
	&+ \frac{u(x+h, t+k) - 2 u(x, t+k) + u(x-h, t+k)}{2 h^{2}}	\\
	&+ O(k^{2}) + O(h^{2})	.
\end{align*}
A formula for $u_{xx}(x, t+k/2)$ has been found by averaging the
three-point formulas for the second derivative at times $t =
k \, j$ and $t = k (j+1)$.

\vspace*{0.10in}

This formula has $O(k^{2}) + O(h^{2})$ accuracy as $h, \, k \to 0$.

\end{frame}
%%

\begin{frame}\frametitle{Stencil}

\begin{center}
\resizebox{1.00\textwidth}{!}{\includegraphics{cranknicheat}}
\end{center}

\end{frame}
%%

\begin{frame}\frametitle{Developing an Implicit Scheme}

The heat equation can be discretized as
\[
\frac{u_{i}^{j+1}-u_{i}^{j}}{k}
	= \frac{u_{i-1}^{j} - 2 u_{i}^{j} + u_{i+1}^{j}}{2 h^{2}} + 
\frac{u_{i-1}^{j+1} - 2 u_{i}^{j+1} + u_{i+1}^{j+1}}{2 h^{2}}
\]
or equivalently as,
\[
%\label{eq:findiff:cranknic}
-r u_{i-1}^{j+1} + 2 (1+r) u_{i}^{j+1} - r u_{i+1}^{j+1}
	= r u_{i-1}^{j} + 2 (1-r) u_{i}^{j} + r u_{i+1}^{j}
\]
for $i = 1, 2, \ldots, N-1$ and where $r = k/h^{2}$.
\begin{itemize}
\item Relative to the point at $(x, t) = (i \, h, (j+1/2) k)$, the first
derivative in $t$ and the second derivative in $x$ have been replaced
with central difference approximations which, in general, are more
accurate than forward (or backward) difference approximations.
\item The three values $u_{i-1}^{j+1}$, $u_{i}^{j+1}$, and $u_{i+1}^{j+1}$
are unknown and must be solved for simultaneously and depend on values
found on the $t = j \, k$ row (which were either given as initial
conditions or were calculated previously).
\end{itemize}

\end{frame}
%%

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

If the boundary conditions are of homogeneous Dirichlet type, the
system of equations can be written in matrix/vector form as
%\label{eq:findiff:tridiagonal}
\[
A(r) \mathbf{u}^{(j+1)}
	= A(-r) \mathbf{u}^{(j)}
\]
where
\[
A(r)
	= \left[ \begin{array}{ccccccc}
2(1+r) & -r & 0 & \cdots & 0 & 0 & 0	\\
-r & 2(1+r) & -r & \cdots & 0 & 0 & 0	\\
\vdots & \vdots & \vdots & {} & \vdots & \vdots & \vdots	\\
0 & 0 & 0 & \cdots & -r & 2(1+r) & -r 	\\
0 & 0 & 0 & \cdots & 0 & -r & 2(1+r)
\end{array} \right]	.
\]

\end{frame}
%%

\begin{frame}\frametitle{Remarks}

\begin{itemize}
\item A unique solution exists if and only
if $A(r)$ is invertible.
\item The inverse of $A(r)$ always exists for any $r > 0$ since the
matrix is \textbf{strictly diagonally
dominant}.
\item The placement of the nonzero entries in this matrix gives
it a special structure known as a
\textbf{tridiagonal matrix}
which is easily inverted.
\item A system of equations of the form $A \, \mathbf{x} = \mathbf{b}$ where
$A$ is an $n \times n$ tridiagonal matrix can be solved by a Gaussian
elimination algorithm requiring $O(n)$ operations (multiplications and
additions).
\item The matrix is
\textbf{symmetric} and the off-diagonal
elements of the matrix are all identical.
\end{itemize}

\end{frame}
%%

\begin{frame}\frametitle{Example}

Consider the following initial boundary value problem for the heat
equation:
\begin{align*}
u_{t}
	&= u_{xx}
\mbox{ for $0 < x < 1$ and $t > 0$}	\\
u(0, t)
	&= u(1, t) = 0
\mbox{ for $t > 0$}	\\
u(x, 0)
	&= x (1-x)^{2}
\mbox{ for $0 < x < 1$.}
\end{align*}
Approximate the solution for $0 \le t \le 1/10$ using the
Crank-Nicolson method with $h = 1/10$ and
$k = 1/100$.

\end{frame}
%%

\begin{frame}\frametitle{Solution (1 of 2)}
\fontsize{8pt}{9.6}\selectfont

For these choices of $h$ and $k$ the parameter $r = 1$.
\[
\left[
\begin{array}{rrrrrr}
4 & -1 & 0 & \cdots & 0 & 0	\\
-1 & 4 & -1 & \cdots & 0 & 0	\\
0 & -1 & 4 & \cdots & 0 & 0	\\
\vdots & \vdots & \vdots & {} & \vdots & \vdots	\\
0 & 0 & 0 & \cdots & 4 & -1 	\\
0 & 0 & 0 & \cdots & -1 & 4
\end{array}
\right]
\left[ \begin{array}{c}
u_{1}^{j+1} \\
u_{2}^{j+1} \\
u_{3}^{j+1} \\
\vdots	\\
u_{N-2}^{j+1}	\\
u_{N-1}^{j+1}
\end{array} \right]
	= \left[ \begin{array}{cccccc}
0 & 1 & 0 & \cdots & 0 & 0	\\
1 & 0 & 1 & \cdots & 0 & 0	\\
0 & 1 & 0 & \cdots & 0 & 0	\\
\vdots & \vdots & \vdots & {} & \vdots & \vdots	\\
0 & 0 & 0 & \cdots & 0 & 1 	\\
0 & 0 & 0 & \cdots & 1 & 0
\end{array} \right]
\left[ \begin{array}{c}
u_{1}^{j} \\
u_{2}^{j} \\
u_{3}^{j} \\
\vdots	\\
u_{N-2}^{j}	\\
u_{N-1}^{j}
\end{array} \right]
\]
Let $u_{i}^{0} = (i h) (1 - i h)^{2}$ with $h = 1/10$ for $i = 1, 2,
\ldots, 9$.

\end{frame}
%%

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

\begin{center}
\resizebox{0.85\textwidth}{!}{\includegraphics{example01}}
\end{center}

\end{frame}
%%

\begin{frame}\frametitle{Robin Boundary Conditions}

Consider the initial boundary value problem,
\begin{align*}
u_{t}
	&= u_{xx}
\mbox{ for $0 < x < 1$ and $t > 0$}	\\
u_{x}(0, t)
	&= u(0, t) - e^{-t}
\mbox{ for $t > 0$}	\\
-u_{x}(1, t)
	&= 2 u(1, t) - t
\mbox{ for $t > 0$}	\\
u(x, 0)
	&= \frac{1}{3} e^{-2 x}
\mbox{ for $0 < x < 1$.}
\end{align*}
Introduce fictitious points $u_{-1}^{j}$ and the discretized form of
the boundary condition at $x = 0$ can be written as
\[
\frac{u_{1}^{j} - u_{-1}^{j}}{2 h}
	= u_{0}^{j} - e^{-j k}	,
\]
which implies the values at the fictitious points are
\[
u_{-1}^{j}
	= -2h u_{0}^{j} + u_{1}^{j} + 2h e^{-j k}	.
\]

\end{frame}
%%

\begin{frame}\frametitle{Crank-Nicolson Formulation}

The boundary condition at $x = 0$ can be discretized as
\[
%\label{eq:findiff:cnRobin0}
(1 + r + h r) u_{0}^{j+1} - r u_{1}^{j+1}
	= (1 - r - h r) u_{0}^{j} + r u_{1}^{j} + h r (e^{-(j+1)k} + e^{-j k})	.
\]
Similarly the boundary condition at $x = 1$ can be approximated using
the central difference formula as
\[
-\frac{u_{N+1}^{j} - u_{N-1}^{j}}{2 h}
	= 2 u_{N}^{j} - j k	.
\]
Solving for the fictitious point where $i = N+1$ yields
\[
u_{N+1}^{j}
	= u_{N-1}^{j} - 4 h u_{N}^{j} + 2 h j k	.
\]
Hence the boundary condition at $x = 1$ can be discretized as
\[
%\label{eq:findiff:cnRobin1}
-r u_{N-1}^{j+1} + (1 + r + 2 h r) u_{N}^{j+1}
	= r u_{N-1}^{j} + (1 - r - 2 h r) u_{N}^{j} + h r (2 j + 1) k	.
\]

\end{frame}
%%

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

The linear system can be expressed as
\[
A(r) \, \mathbf{u}^{j+1}
	= A(-r) \, \mathbf{u}^{j} +
\left[ \begin{array}{c}
h r (e^{-(j+1)k} + e^{-j k})	\\
0	\\
\vdots	\\
0	\\
h r (2j+1) k
\end{array} \right]
\]
where
\[
A(r)
	= \left[ \begin{array}{ccccccc}
1+r+h r & -r & 0 & \cdots & 0 & 0 & 0	\\
-r & 2(1+r) & -r & \cdots & 0 & 0 & 0	\\
\vdots & \vdots & \vdots & {} & \vdots & \vdots & \vdots	\\
0 & 0 & 0 & \cdots & -r & 2(1+r) & -r	\\
0 & 0 & 0 & \cdots & 0 & -r & 1+r+2h r
\end{array} \right]	.
\]

\end{frame}
%%

%\begin{frame}\frametitle{Homework}

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

%\end{frame}
%%

%%%
\end{document}

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