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?

NDSolve currently does not support anti-periodic boundary conditions; only periodic boundary conditions are supported. — user21, Mar 09 '16 at 02: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? — Thomas, Mar 09 '16 at 03: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. — user21, Mar 09 '16 at 04:48

xzczd · Answer 1 · 2019-02-01T15:59:15.523

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"}]

Mathematica graphics

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

Mathematica graphics

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}].

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

1 Answers1

Linked