As a amateur in applied mathematics, I find that my biggest frustration is not knowing how to implement the models presented in papers I read. The pretty figures in the papers seemed so compelling but at the same time unreachable. In this section I want to try and step through some of the key ideas behind the code that was used to generate the gifs scattered around the website. This is an extremely basic introduction to the numerical methods to solve PDE's, but should be a good starting point.

First, I will go through the (explicit and implicit) numerical solution of the 1D Heat Equation, and then I will step through the explicit numerical solution for a 2D Heat Equation.

Before we begin let's recall what the Heat Equation is. This should be familiar from our discussion of diffusion as they essentially the same.

The Heat equation (also the diffusion equation) is: \begin{equation} \label{heatEqn} \frac{\partial u}{\partial t}\, - \,D_u \frac{\partial ^2 u}{\partial x^2}\, =\, 0, \end{equation} where $u$ is temperature, and $D_u$ is a constant that determines how fast heat spreads out. Imagine we have an insulated rod with initial condition $u(x,0) = u_0$ boundary conditions $u(x,t) = u(L,t) = 0$. We want to solve the PDE in Equation [\ref{heatEqn}].

The trick is to assume that the solutions are separable, so \begin{equation} \label{eqn:heatEqnsep1} u(x,t) = T(t)X(x) \end{equation}

We can substitute this solution into Equation [\ref{heatEqn}] to get \begin{equation} \label{eqn:heatEqnsep2} \underbrace{T'(t)X(x)}_{\Large =\frac{\partial u}{\partial t}} = D_u\underbrace{ T(t)X''(x)}_{\Large =\frac{\partial^2 u}{\partial x^2}}. \end{equation}

Rearranging terms we get

\begin{equation} \label{eqn:heatEqnsep3} \frac{ T'(t)}{D_u T(t)} = \frac{X''(x)}{X(x)}. \end{equation}

Note that the LHS (left hand side) is a function of time where the RHS (right hand side) is a function of position. In order for this equality to be satisfied for all $x$ and $t$, we must have that they are both equal to a constant. We let

\begin{equation} \label{eqn:heatEqnsep4} \frac{ T'(t)}{D_u T(t)} = \frac{X''(x)}{X(x)} = -\lambda^2 \text{, } \quad \lambda > 0. \end{equation}

Solving for each variable we get that \begin{align} T(t) &= Ae^{-\lambda^2 D_u t}\\ X(x) &= Bcos(\lambda x) + Csin(\lambda x) \end{align}

Here, $A$ is a positive real constant. We can apply our boundary conditions to find what $B$ and $C$ should be. Since the edges at $x=0$ is clamped at zero ($u(0,t)=0$), then $B=0$. The rod is also clamped at zero at $x=L$, and since $C\neq 0$, then $sin(\lambda L) =0$ This means that $\lambda_n = \frac{n\pi}{L}$ where $n\in \mathbb{Z}$.

Note that:

• $n=0$ give a trivial solution
• $n\in \mathbb{Z}^-$ just makes all the coefficients negative, so we can just consider the positive values of n.

As a result, we get the following eigenfunctions:

\begin{equation} \label{eqn:heatEqneigen1} u_n(x,t) =A_n e^{-\lambda_n^2 D_u t} sin(\lambda_n x) \text{, } \quad n=1,2,..., \end{equation} with eigenvalues

\begin{equation} \label{eqn:heatEqneigen2} \lambda_n = \frac{n\pi}{L} \end{equation}

We are not done yet, because we still need to find the conditions s.t. we satisfy the initial conditions. Since linear combinations of solutions will still be solutions, then \begin{equation} \label{eqn:heatEqnfour1} u(x,t) = \sum_{n=1}^\infty A_n e^{-\lambda_n^2 D_u t} sin(\lambda_n x), \end{equation} is still a solution. We would like it then, if we could express $u_0$ as a sum of sine functions as shown in Equation [\ref{eqn:heatEqnfour1}]. This is true for most $u_0$ and it is the theory of Fourier Series.

Although I will not get into the details here, it can be shown that the Fourier sine coefficients are:

\begin{equation} \label{eqn:heatEqneigen3} A_n = \frac{2}{L} \int_0^L u_0(x)sin(\lambda_n x) dx \end{equation}

However, in general, we look for numerical solutions to PDE's and this can be done using the finite-difference method where we replace every derivative by a finite difference approximation. We consider a few examples for the Heat Equation below. (For a more complete discussion and more Matlab examples and exercises see [(kharab_2002)]).

The key idea behind the forwrd finite difference method is that we can rewrite: \begin{gather} \label{eqn:uxx} \frac{\patial^2 u}{\partial x^2} \approx

[(:kharab_2002» author: Kharab A and Guenther R.B. ref-author: Kharab and Guenther title: An Introduction to Numerical Methods: A Matlab Approach publisher: Chapman and Hall year: 2002 )]