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

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$ ?

• Do you have Dirichlet or Neumann boundary conditions? The number of ghost cells depends on the order of the approximation for you Neumann boundary conditions. – ilciavo Jul 8 '15 at 23:31
• @ilciavo, the quesition is about periodic boundary conditions. – Bill Barth Jul 8 '15 at 23:59

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

• Great and simple solution!! in case someone needs it here the implementations in Python – ilciavo Jul 10 '15 at 11:36
• Great ! I didn't want to be too specific about the method, as I for myself use semi-discretization for instance. But for the equation, you confirm that there is no need to specify a condition on the $x$-derivative ? More precisely, would such a condition lead to an ill-posed problem ? – bela83 Jul 10 '15 at 15:59
• @bela83 You are correct that there is no need to specify anything more than the initial condition. Doing so would result in an over-determined system. All you need to do is be a little careful near the end points of the interval to be sure you properly wrap things around periodically. There are many valid ways to do this. – Doug Lipinski Jul 10 '15 at 16:21

According to this you should impose periodic boundary conditions as:

$$u(0, t) = u(1, t) \\ u_x(0, t) = u_x(1, t)$$

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

$$\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}$$

Solving the system of equations

$$\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}$$

Where $$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}$$

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:

$$u_1 - u_{N} = 0 \\ \frac{u_{2}-u_{0}}{2 \Delta x} - \frac{u_{N+1}-u_{N-1}}{2 \Delta x} = 0$$

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{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}$$

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.

• I don't understand the statement : "add two additional cells $u_{-1}$ and $u_N$" because $u_N$ is already a point in $]0,1[$. I have in mind $x_k=\frac{k}{N+1}$ so that $x_N=\frac{N}{N+1}$. – bela83 Jul 9 '15 at 5:51
• It's just a matter of indexing. You start with $N$ cells (or points) from $u_0$ to $u_{N-1}$ and you add two cells $u_{-1}$ and $u_{N}$. If you have cells going from $u_1$ to $u_N$ then you need to add cells at $u_0$ and $u_{N+1}$ – ilciavo Jul 9 '15 at 5:59
• OK then I don't understand the two more equations ! The first should be (with the notation from the question): $u_0=u_{N+1}$ right ? But I don't understand the second one: why not choose a centered difference approximation ? Finally, that makes $N+1$ unknowns, if first and last values are equal. Please compare to the situation with (homogeneous) Neumann BC in the question. – bela83 Jul 9 '15 at 6:32
• I've change the notation. It depends on the order of approximation for $u_x$. The first comes from $u(0,t)=u(1,t)$ and the second from $u_x(0,t)=u_x(1,t)$ – ilciavo Jul 9 '15 at 10:39
• u1=uN adds an additional restriction(equation u1−uN=0) to your system – ilciavo Jul 9 '15 at 12:33