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

Question

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?

Answer 1

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,
   please find it in the link above. *)
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}].