# Parabolic PDEs

A classic example of a parabolic partial differential equation (PDE) is the one-dimensional unsteady heat equation:
\begin{equation}
\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial t^2} 
\end{equation}
where $T(x,t)$ is the temperature varying in space and time, and $\alpha$ is the thermal diffusivity: $\alpha = k / (\rho c_p)$, which is a constant.

We can solve this using finite differences to represent the spatial derivatives and time derivatives separately.
First, let's rearrange the PDE slightly:
\begin{equation}
\frac{\partial^2 T}{\partial x^2} = \frac{1}{\alpha} \frac{\partial T}{\partial t}
\end{equation}

## Explicit scheme

Let's use a *central difference* for the spatial derivative with a spacing of $\Delta x$, and a *forward difference* for the time derivative with a time-step size of $\Delta t$. With these choices, we can obtain an approximation to the PDE that applies at time $t^k$ and spatial location $x_i$:
\begin{equation}
\frac{T_{i-1}^k - 2 T_i^k + T_{i+1}^k}{\Delta x^2} = \frac{1}{\alpha} \left( \frac{T_i^{k+1} - T_i^k}{\Delta t} \right)
\end{equation}
where $T_i^k$ is the temperature at time $t^k$ and spatial location $x_i$. The following figure shows the stencil of points involved in the PDE, for a domain with five points in the $x$-direction.

:::{figure-md} fig-stencil-explicit


Stencil for explicit solution to heat equation
:::

To solve the heat equation for a one-dimensional domain over $0 \leq x \leq L$, we will need both initial conditions at $t = 0$ and boundary conditions at $x=0$ and $x=L$ (for all time). In terms of our nodal values, this means we need $T_i^{k=1}$ for $i = 1 \ldots n$, where $n$ is the number of points, as well as information about $T_1^k$ and $T_n^k$ for all times $k$.

We can rearrange the above equation to obtain our recursion formula:
\begin{equation}
T_i^{k+1} = \left( T_{i+1}^k + T_{i-1}^k \right) \frac{\alpha \Delta t}{\Delta x^2} + T_i^k  \left( 1 - 2 \frac{\alpha \Delta t}{\Delta x^2} \right) \;.
\end{equation}
This is an **explicit** scheme in time, similar to the Forward Euler method we used for ordinary differential equations, and like that method it may have stability issues. The combination of terms we see repeated is also known as the Fourier number: $\text{Fo} = \frac{\alpha \Delta t}{\Delta x^2}$, and governs the stability.
We can rewrite the recursion formula using this:
\begin{equation}
T_i^{k+1} = \left( T_{i+1}^k + T_{i-1}^k \right) \text{Fo} + T_i^k  \left( 1 - 2 \text{Fo} \right) \;.
\end{equation}

The term in parentheses there must be greater than or equal to zero for stability ($1 - 2 \text{Fo} \geq 0$); if not, the solution may become unstable and blow up. This gives us some conditions on our choice of step sizes:
\begin{align}
1 - 2 \text{Fo} &\geq 0 \\
1 & \geq 2 \text{Fo} \\
\text{Fo} &\leq \frac{1}{2} \\
\frac{\alpha \Delta t}{\Delta x^2} &\leq \frac{1}{2}
\end{align}
This is the stability criterion for the explicit method: the Fourier number must be smaller than 0.5.
For a given thermal diffusivity and chosen spatial step size, this also gives us the limit on the time-step size: $\Delta t \leq \Delta x^2 / (2 \alpha)$.

Let's look at an example where the initial temperature is 200, the temperature at the boundaries are 50, the thermal diffusivity is $\alpha = 2.3 \times 10^{-1}$ m$^2 /$ s, and $L = 1$. 
In other words,
\begin{align}
T(x, t=0) &= 200 \\
T(x=0, t) &= 50 \\
T(x=L, t) &= 50
\end{align}
We'll integrate out to $t = 1$, using a Fourier number of 0.25 to be comfortably below the stability limit (`Fo = 0.25`):


```matlab
clear all

dx = 0.1;
alpha = 2.3e-1;

% for stability, set the Fourier number at 0.25 (half the stability limit of 0.5)
Fo = 0.25;
% then choose the time step based on the Fourier number
dt = Fo * dx^2 / alpha;

x = [0 : dx : 1]; n = length(x);
t = [0 : dt : 1]; m = length(t);

T = zeros(m, n);

% initial conditions
T(1,:) = 200;

plot(x, T(1,:))
axis([0 1 50 200]);
xlabel('Distance'); ylabel('Temperature');
F(1) = getframe(gcf);

for k = 1 : m - 1
    for i = 1 : n
        if i == 1
            T(k+1, 1) = 50;
        elseif i == n
            T(k+1, n) = 50;
        else
            T(k+1, i) = (T(k,i+1) + T(k,i-1))*Fo + T(k,i)*(1 - 2*Fo);
        end
    end
    plot(x, T(k+1,:))
    axis([0 1 50 200]);
    xlabel('Distance'); ylabel('Temperature');
    F(k+1) = getframe(gcf);
end
close

%% If you are working interactively, you can use this to make a movie in Matlab
%fig = figure;
%movie(fig, F, 2)

%% This generates a GIF of the results (for use in Jupyter Notebook)
filename = 'parabolic_animated.gif';
for i = 1 : length(F)
    im = frame2im(F(i)); 
    [imind,cm] = rgb2ind(im,256); 
    % Write to GIF
    if i == 1 
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3); 
    else 
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3); 
    end
end
```

:::{figure-md} fig-transient-heat


Animated solution to 1D transient heat transfer PDE
:::

This shows the temperature decaying exponentially from the initial conditions, constrained by the boundary conditions.

What happens if we tried to use a Fourier number larger than 0.5, or arbitrarily chose a time-step size that was too large (and resulted in $\text{Fo} > 0.5$)?


```matlab
clear all

dx = 0.1;
alpha = 2.3e-1;

%% Purposely choose a Fourier number that is past the stability limit:
Fo = 0.75;
dt = Fo * dx^2 / alpha;

x = [0 : dx : 1]; n = length(x);
t = [0 : dt : 1]; m = length(t);

T = zeros(m, n);

% initial conditions
T(1,:) = 200;

plot(x, T(1,:))
axis([0 1 50 200]);
xlabel('Distance'); ylabel('Temperature');
F(1) = getframe(gcf);

for k = 1 : m - 1
    for i = 1 : n
        if i == 1
            T(k+1, 1) = 50;
        elseif i == n
            T(k+1, n) = 50;
        else
            T(k+1, i) = (T(k,i+1) + T(k,i-1))*Fo + T(k,i)*(1 - 2*Fo);
        end
    end
    plot(x, T(k+1,:))
    xlabel('Distance'); ylabel('Temperature');
    F(k+1) = getframe(gcf);
end
close

%% If you are working interactively, you can use this to make a movie in Matlab
%fig = figure;
%movie(fig, F, 2)

%% This generates a GIF of the results (for use in Jupyter Notebook)
filename = 'parabolic_unstable_animated.gif';
for i = 1 : length(F)
    im = frame2im(F(i)); 
    [imind,cm] = rgb2ind(im,256); 
    % Write to GIF
    if i == 1 
        imwrite(imind,cm,filename,'gif', 'Loopcount',2, 'DelayTime',1e-3); 
    else 
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3); 
    end
end
```

:::{figure-md} fig-unstable-solution


Animated unstable solution
:::

In this case, the solution becomes unstable and blows up, leading to unphysical results.

For this explicit scheme, the choice of $\Delta t$ is limited by the stability criterion. This means that we may be stuck using a small time-step size.

Rather than being forced to use a very small time-step size, we can also explore *implicit schemes* that are unconditionally stable.

## Implicit scheme




```matlab
clear all

alpha = 2.3e-1;
dx = 0.1;
x = [0 : dx : 1]; n = length(x);

% choose a Fourier number that is deliberately past the explicit method stability limit
Fo = 0.75;
dt = Fo * dx^2 / alpha;

t = [0 : dt : 1]; m = length(t);

T = zeros(m, n);

% Initial conditions
T(1,:) = 200;

plot(x, T(1,:))
axis([0 1 50 200]);
xlabel('Distance'); ylabel('Temperature');
F(1) = getframe(gcf);

for k = 1 : m - 1
    A = zeros(n,n);
    b = zeros(n,1);
    for i = 1 : n
        if i == 1
            A(1,1) = 1;
            b(1) = 50;
        elseif i == n
            A(n,n) = 1;
            b(n) = 50;
        else
            A(i,i-1) = Fo;
            A(i,i) = -2*Fo - 1;
            A(i,i+1) = Fo;
            b(i) = -T(k,i);
        end
    end
    
    T(k+1, :) = A \ b;
    plot(x, T(k+1,:))
    axis([0 1 50 200]);
    xlabel('Distance'); ylabel('Temperature');
    F(k+1) = getframe(gcf);
end
close

%% This generates a GIF of the results (for use in Jupyter Notebook)
filename = 'parabolic_implicit_animated.gif';
for i = 1 : length(F)
    im = frame2im(F(i)); 
    [imind,cm] = rgb2ind(im,256); 
    % Write to GIF
    if i == 1 
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3); 
    else 
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3); 
    end
end
```

:::{figure-md} fig-implicit-solution


Solution to 1D heat equation with implicit method. Fo = 0.75
:::
