# My code for diffusion does not conserve mass

I am trying to numerically integrate a diffusion equation with ad advective flux around x=0, with zero flux at the boundaries. The problem is that although the equation should conserve the total mass, this is not the case in the numerical solution of NDSolve.

How can I get the conserved mass to work?

The problem is the following: I write the local flux as

J=-D[y[t,x],x]-G[x]y[t,x].


y[t,x] is the density at time t in point x and G[x] is G[x]=x/(Abs[x]^3+1).

I write the Mathematica code as follows:

G[x_] := x/(Abs[x]^3 + 1);
D[ D[y[t, x], x] + G[x] y[t, x], x]

y[t,x]/(1+Abs[x]^3)-(3 x Abs[x]^2 y[t,x] (Abs^\[Prime])[x])/(1+Abs[x]^3)^2+
(x Derivative[0,1][y][t,x])/(1+Abs[x]^3)+Derivative[0,2][y][t,x]


Here I substitute Abs'[x] with Sign[x] otherwise NDSolve does not work.

Lastly, I compute:

sol = NDSolve[{D[y[t, x], t] ==
y[t, x]/(1 +
Abs[x]^3) - (3 x Abs[x]^2 y[t, x] Sign[x])/(1 +
Abs[x]^3)^2 + (x Derivative[0, 1][y][t, x])/(1 + Abs[x]^3) +
Derivative[0, 2][y][t, x], y[0, x] == 1,
Derivative[0, 1][y][t, -1] + G[-1] y[t, -1] == 0,
Derivative[0, 1][y][t, 1] + G[1] y[t, 1] == 0},
y, {t, 0, 10}, {x, -1, 1}, MaxStepSize -> 0.1, MaxSteps -> 100000]


If I now plot the solution at different times, though, the total amount of mass in the closed interval [-1,1] grows a lot, so clearly I can't trust this solution.

-
I'm not sure what y^(0,1) means, but your function for G is not peaked at zero. Furthermore, you define it two different ways: x/(Abs[x]^2+1) and x/(Abs[x]^3+1). Perhaps a Lorentzian like 1/(Abs[x]^2+1) would do. –  Daniel Flatin Oct 9 at 12:18
You're right, thank you. The description is wrong, but the problem is with the G[x] that now appears in the above text, I've corrected it. Cheers! –  Andrea Oct 9 at 13:06
y^(0,1) is the derivative in space, y^(1,0) is the derivative in time, I'm not sure why the copy-paste pasted it this way –  Andrea Oct 9 at 13:19
I think you should use Derivative[0,1][y][t,-1] and similar. –  b.gatessucks Oct 9 at 13:25
Most numerical solvers are not conservative. For example, plain 4th order Runge-Kutta is dissipative. It takes care to generate a conservative method, and, unfortunately, I never learned what was necessary to generate them. I believe Crank-Nicholson is unitary, i.e. length preserving, but it has been a while. –  rcollyer Oct 9 at 13:27