# Boundary conditions (periodic/anti-periodic) for diffusion equation and NDSolve

I am trying to solve a simple diffusion equation in Mathematica. The problem is inspired by spin diffusion, which is why I consider an initial condition that has both positive and negative density rho[x,t].

Here is a simple mathematica code that appears to work correctly

eq1 = D[rho[x, t], t] == Dif* D[rho[x, t], {x, 2}];
a = 3;
iv = {rho[x, 0] == x*Exp[-x^2], rho[-a, t] == 0, rho[a, t] == 0};
Dif = 0.1;
sl1 = NDSolve[{eq1, iv}, {rho[x, t]}, {x, -a, a}, {t, 0, 10}]
Plot3D[rho[x, t] /. sl1, {x, -a, a}, {t, 0, 10}, PlotRange -> {-0.5, 0.5}]


but it gives a warning because the initial condition is only approximately consistent with the boundary condition.

I thought I should do better, and tried to implement periodic boundary conditions. This also works

iv2 = {rho[x, 0] == x*Exp[-x^2], rho[-a, t] == rho[a, t],
Derivative[1, 0][rho][-a, t] == Derivative[1, 0][rho][a, t]};
sl2 = NDSolve[{eq1, iv2}, {rho[x, t]}, {x, -a, a}, {t, 0, 10}]
Plot3D[rho[x, t] /. sl2, {x, -a, a}, {t, 0, 10}, PlotRange -> {-0.5, 0.5}]


but it is still not quite right, because for my problem I should implement anti-periodic boundary conditions for rho[x,t], and periodic boundary conditions for the current D[rho[x,t],x]. So this should be the best solution

iv3 = {rho[x, 0] == x*Exp[-x^2], rho[-a, t] == -rho[a, t],
Derivative[1, 0][rho][-a, t] == Derivative[1, 0][rho][a, t]};
sl3 = NDSolve[{eq1, iv3}, {rho[x, t]}, {x, -a, a}, {t, 0, 10}]


except, it does not work, it produces an error

NDSolve::bcedge: "Boundary condition rho[-3,t]==-rho[3,t] is not
specified on a single edge of the boundary of the computational domain.


which I cannot make sense of. What is the problem here?

• NDSolve currently does not support anti-periodic boundary conditions; only periodic boundary conditions are supported. Mar 9, 2016 at 2:00
• Thank you, I should have guessed that (the error message is a little misleading). Is there a good place where this sort of thing is documented? Mar 9, 2016 at 3:55
• Probably the option you have is the Advanced Numerical Differential Equation Solving in the Wolfram Language. You get there from the NDSolve reference page by selecting the Tutorials icon and selecting the tutorial. I am, however, not sure this is mentioned in the tutorial. Mar 9, 2016 at 4:48

As mentioned by user21 in the comment above, NDSolve can't handle anti-periodic b.c. at the moment, so let's discretize the PDE to a set of ODEs, I'll use pdetoode for the task:

tend = 10;
a = 3;
Dif = 0.1;
eq = D[rho[x, t], t] == Dif*D[rho[x, t], {x, 2}];
ic = rho[x, 0] == x*Exp[-x^2];
bc = {rho[-a, t] == -rho[a, t],
Derivative[1, 0][rho][-a, t] == Derivative[1, 0][rho][a, t]};
points = 25;
grid = Array[# &, points, {-a, a}];
difforder = 4;
(* Definition of pdetoode isn't included in this code piece,
ptoofunc = pdetoode[rho[x, t], t, grid, difforder];
del = #[[3 ;;]] &;
{ode, odeic} = del@ptoofunc@# & /@ {eq, ic};
odebc = Map[ptoofunc, bc, {2}];
sollst = NDSolveValue[{ode, odeic, odebc}, rho /@ grid, {t, 0, tend}];
sol = rebuild[sollst, grid, -1];

Plot[{sol[a, t], sol[-a, t]}, {t, 0, tend}, PlotRange -> All,
PlotLegends -> {"Right boundary", "Left boundary"}]


Plot3D[sol[x, t], {x, -a, a}, {t, 0, tend}, PlotRange -> All]


Notice pdetoode currently can't discretize periodic/anti-periodic b.c. e.g. u[-1, t] == u[1, t] directly, so I circumvent the issue with Map[ptoofunc, bc, {2}].