Periodic boundary condition for the heat equation in ]0,1[

Question

Let us consider a smooth initial condition and the heat equation in one dimension : $$ \partial_t u = \partial_{xx} u$$ in the open interval $]0,1[$, and let us assume that we want to solve it numerically with finite differences.

I know that for my problem to be well-posed, I need to endow it with boundary conditions at $x=0$ and $x=1$. I know that Dirichlet or Neumann work well.

If I have in the first case $N$ interior points $x_k=\frac{k}{N+1}$ for $k=1,\cdots,N$, then I have $N$ unknowns: $u_k=u(x_k)$ for $k=1,\cdots,N$, because $u$ is prescribed at the boundaries.

In the second case I really have $N+2$ unknowns $u_0,\cdots,u_{N+1}$, and I know how to use the (homogeneous) Neumann BC to discretize the laplacian at the border, for example with the adjunction of two fictious points $x_{-1}$ and $x_{N+2}$ and the equalities :

$$\frac{u_1-u_{-1}}{2 h} = 0 = \frac{u_{N+2}-u_N}{2 h}$$

My question is about periodic BC. I have the feeling that I could use one equation, namely $$u(0) = u(1)$$ but maybe two, and then I would use $$\partial_x u(0) = \partial_x u(1)$$

but I am not sure. I don't know either how many unknowns I should have. Is it $N+1$ ?

Answer 1

The best way to do this is (as you said) to just use the definition of periodic boundary conditions and set up your equations correctly from the start using the fact that $u(0)=u(1)$. In fact, even more strongly, periodic boundary conditions identify $x=0$ with $x=1$. For this reason, you should only have one of these points in your solution domain. An open interval does not make sense when using periodic boundary conditions since there is no boundary.

This fact means that you should not place a point at $x=1$ since it is the same as $x=0$. Discretizing with $N+1$ points, you then use the fact that by definition, the point to the left of $x_0$ is $x_N$ and the point to the right of $x_N$ is $x_0$.

Your PDE can then be discretized in space as $$ \frac{\partial}{\partial t} \left[\begin{array}{c} x_0 \\ x_1 \\ \vdots \\ x_N \end{array}\right] = \frac{1}{\Delta x^2} \left[\begin{array}{c} x_N - 2x_0 + x_1 \\ x_0 - 2x_1 + x_2 \\ \vdots \\ x_{N-1} - 2x_N + x_0 \end{array}\right] $$

This can be written in matrix form as $$ \frac{\partial}{\partial t}\vec{x} = \frac{1}{\Delta x^2} \mathbf{A} \vec{x} $$ where $$ \mathbf{A} = \left[\begin{array}{c} -2 & 1 & 0 & \cdots & 0 & 1 \\ 1 & -2 & 1 & 0 & \cdots & 0 \\ & \ddots & \ddots & \ddots \\ && \ddots & \ddots & \ddots \\ 0 & \cdots & 0 & 1 & -2 & 1 \\ 1 & 0 & \cdots & 0 & 1 & -2 \end{array}\right]. $$

Of course there is no need to actually create or store this matrix. The finite differences should be computed on the fly, taking care to handle the first and last points specially as needed.

As a simple example, the following MATLAB script solves $$ \partial_t u = \partial_{xx}u + b(t,x) $$ with periodic boundary conditions on the domain $x\in[-1,1)$. The manufactured solution $u_\text{Ref}(t,x) = \exp(-t)\cos(5\pi x)$ is used, meaning $b(t,x) = (25\pi^2-1)\exp(-t)\cos(5\pi x)$. I used forward Euler time discretization for simplicity and computed the solution both with and without forming the matrix. The results are shown below.

clear

% Solve: u_t = u_xx + b
% with periodic boundary conditions

% analytical solution:
uRef = @(t,x) exp(-t)*cos(5*pi*x);
b = @(t,x) (25*pi^2-1)*exp(-t)*cos(5*pi*x);

% grid
N = 30;
x(:,1) = linspace(-1,1,N+1);

% leave off 1 point so initial condition is periodic
% (doesn't have a duplicate point)
x(end) = [];
uWithMatrix = uRef(0,x);
uNoMatrix = uRef(0,x);

dx = diff(x(1:2));
dt = dx.^2/2;

%Iteration matrix:
e = ones(N,1);
A = spdiags([e -2*e e], -1:1, N, N);
A(N,1) = 1;
A(1,N) = 1;
A = A/dx^2;

%indices (left, center, right) for second order centered difference
iLeft = [numel(x), 1:numel(x)-1]';
iCenter = (1:numel(x))';
iRight = [2:numel(x), 1]';

%plot
figure(1)
clf
hold on
h0=plot(x,uRef(0,x),'k--','linewidth',2);
h1=plot(x,uWithMatrix);
h2=plot(x,uNoMatrix,'o');
ylim([-1.2, 1.2])
legend('Analytical solution','Matrix solution','Matrix-free solution')
ht = title(sprintf('Time t = %0.2f',0));
xlabel('x')
ylabel('u')
drawnow

for t = 0:dt:1
    uWithMatrix = uWithMatrix + dt*( A*uWithMatrix + b(t,x) );
    uNoMatrix = uNoMatrix + dt*(  ( uNoMatrix(iLeft) ...
                                - 2*uNoMatrix(iCenter) ...
                                  + uNoMatrix(iRight) )/dx^2 ...
                                + b(t,x) );
    set(h0,'ydata',uRef(t,x))
    set(h1,'ydata',uWithMatrix)
    set(h2,'ydata',uNoMatrix)
    set(ht,'String',sprintf('Time t = %0.2f',t))
    drawnow
end

Answer 2

According to this you should impose periodic boundary conditions as:

\begin{equation} u(0, t) = u(1, t) \\ u_x(0, t) = u_x(1, t) \end{equation}

One way of discretising the Heat Equation implicitly using backward Euler is

\begin{equation} \frac{u^{n+1}-u_{n}}{\Delta t}= \frac{u^{n+1}_{i+1}-2u^{n+1}_{i}+u^{n+1}_{i+1}}{\Delta x^2} \end{equation}

Solving the system of equations

\begin{equation} \begin{bmatrix} I-\frac{\Delta t}{\Delta x^2}A \end{bmatrix} \begin{bmatrix} u^{n+1}_1 \\ u^{n+1}_1 \\ \vdots \\ \\ \\ u^{n+1}_{N} \end{bmatrix} = \begin{bmatrix} u^n_1 \\ u^n_2 \\ \vdots \\ \\ \\ u^n_{N} \end{bmatrix} \end{equation}

Where \begin{equation} A= \begin{bmatrix} -2 & 1 & 0 & 0 & 0 & \ldots & 0 \\ 1 & -2 & 1 & 0 & 0 &\ldots & 0 \\ 0 & 1 & -2& 1 & 0 & \ldots & 0 \\ 0 & 0 & 1& -2 & 1 & \ldots & 0 \\ 0 & 0 & 0& 1 & -2 & \ldots & 0 \\ \vdots &\vdots &\vdots &\vdots & \ddots &\ddots &\vdots \\ 0 & 0 & 0 & \ldots &0 &1 &-2 \end{bmatrix} \end{equation}

Your periodic boundary conditions can be included by adding two more equations and two additional(ghost) cells $u_{0}$ and $u_{N+1}$ such that:

\begin{equation} u_1 - u_{N} = 0 \\ \frac{u_{2}-u_{0}}{2 \Delta x} - \frac{u_{N+1}-u_{N-1}}{2 \Delta x} = 0 \end{equation}

According to Section 2.11 LeVeque this gives you a 2nd order accuracy for $u_x$

Finally your system of equations will look like: \begin{equation} \begin{bmatrix} 0 & 1 & 0 & 0 & \ldots & 0 & -1 & 0 \\ -1 & 0 & 1 & 0 & \ldots & 1 & 0 & -1 \\ 0 & & & & & & & 0\\ 0 & & & & & & & 0\\ 0 & & & I-\frac{\Delta t}{\Delta x^2}A & & & & 0\\ 0 & & & & & & & 0\\ 0 & & & & & & & 0\\ 0 & & & & & & & 0 \end{bmatrix} \begin{bmatrix} u^{n+1}_{0}\\ u^{n+1}_1 \\ u^{n+1}_2 \\ \\ \\ \\ \\ u^{n+1}_{N} \\ u^{n+1}_{N+1} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ u^n_1 \\ u^n_2 \\ \\ \\ \\ \\ u^n_{N} \end{bmatrix} \end{equation}

Which gives you N+2 equations and N+2 unknowns.

You can also get rid of the first to equations and the ghost cells and arrive at a system of N equations and N unknowns.