# Difficulty in getting correct Gaussian curve for diffusion of point source

I want to solve diffusion of a point source numerically and check it against analytical solution. first I define initial profile,

ClearAll["Global*"];
deltaPos = 0;
InitComp = 100;
g[x_] := Piecewise[{{0, x < deltaPos}, {100, x == deltaPos}, {0,
x > deltaPos}}]


and the solve the differential equation,

pde = D[u[x, t], t] - 0.2 D[D[u[x, t], {x, 1}], {x, 1}] == 0;
sol = NDSolve[{pde, u[x, 0] == g[x], u[-10, t] == 0, u[10, t] ==0},
u[x, t], {x, -10, 10}, {t, 0, 20}];


Now, I expect that total mass of the system is essentially conserved, since diffusion is yet to reach boundary points. Initial mass of the system is

NIntegrate[u[x, t] /. sol /. {t -> 0}, {x, -10, 10}]


which outputs $0$, which is acceptable since total mass is concentrated at a single point. Now to check mass after some time(essentially any time t=0.0001,0.1,1,10),

NIntegrate[u[x, t] /. sol /. {t -> 1}, {x, -10, 10}]


outputs to $0.2$, which is confusing. I think integration should produce $100$, which is the total initial mass of the system.

Analytical solution after time t is (with $D = 0.2, C_0 = 100$), \begin{eqnarray} f(x,t) = \frac{C_0}{\sqrt{4\pi Dt}}\exp\left(\frac{-x^2}{4Dt}\right) \end{eqnarray}

Comparing the results with analytical solution,

AnaS[x_] := 100*Exp[-x*x/(4*0.2*10)]/(Sqrt[4*3.14159265*0.2*10])
Plot[AnaS[x], {x, -10, 10}]
Plot[u[x, t] /. sol /. {t -> 10}, {x, -10, 10}]


tells that numerical solution is off by exactly factor of $500$. What am I missing here? Also another question regarding the boundary conditions. How can set boundary conditions like $\frac{\partial u}{\partial x} = 0$ at the boundaries? I tried

D[u[10,t],x]==0


but it did not work.

• I don't think your way of defining a point source is correct. The peak has a finite height but a width of zero, so it basically contains nothing. Apr 29 '15 at 11:20
• @SjoerdC.deVries, I know in c++, this way of defining point source works just fine. I am not sure with Mathematica. Apr 29 '15 at 11:23
• Apr 29 '15 at 11:53
• the boundary conditions for x=-10 and x=10 are completely invalid. A peculiarity of the diffusion equation is that perturbations propagate instantaneously. Another thing is the way, your "delta-function" is defined. NDSolve gives an message that it is looking at 10000 grid points for x. The interval for x is {-10,10} thus it might mistakenly interpret the single point of g[x]=100 as a peak of width around 0.002 and height up to 100 giving a total mass of around 0.2 (500 smaller than the desired 100. If I increase the x range to {-20,20} the solution becomes twice as large. Apr 29 '15 at 15:02
• @alekhine, *I meant any code using finite difference method *. This requires, that a) deltaPos is among the sampling points, b) the effective mass you'll get is dependent on the density of the sampling points. Apr 29 '15 at 16:27

The reason why you get a factor of 500 of I have explained in my comment. Let's replace your g with a better behaved function:

mass = 100;
width = 10^-2;
g[x_] := mass/width HeavisideLambda[x/width]


This is a triangular peak with area 100 and basewidth 0.01.

Now let's impose the desired boundary conditions of $\partial u / \partial x =0$ at the boundaries:

Derivative[1, 0][u][-10, t] == 0
Derivative[1, 0][u][10, t] == 0


pde = D[u[x, t], t] - 0.2 D[D[u[x, t], {x, 1}], {x, 1}] == 0;
sol = NDSolve[{pde, u[x, 0] == g[x], Derivative[1, 0][u][-10, t] == 0,
Derivative[1, 0][u][10, t] == 0}, u[x, t], {x, -10, 10}, {t, 0, 20}]


This will of course complain due to the absence of a continuous derivative (much less, a continuous second derivative) of g[x], but it will come up with something.

Plot3D[u[x, t] /. sol, {x, -5, 5}, {t, 0, 20}, PlotRange -> Full] Huge peak at the beginning (just a finite approximation of a delta-function).

Plot3D[u[x, t] /. sol, {x, -10, 10}, {t, 0, 20}] Something saner, as diffusion takes place.

{-10, 10} is quite a large box. It's hard to see, that the process of diffusion is really contained within it. Let's confine the region a bit and give more time to develop.

sol = NDSolve[{pde, u[x, 0] == g[x], Derivative[1, 0][u][-3, t] == 0,
Derivative[1, 0][u][3, t] == 0}, u[x, t], {x, -3, 3}, {t, 0, 50}]
Plot3D[u[x, t] /. sol, {x, -3, 3}, {t, 0, 30}, PlotRange -> {0, 30}] Now the peak fills the box up to about 15.

NIntegrate[u[x, t] /. (First@sol) /. t -> 0, {x, -0.1, 0.1}]


100.08

I've chosen this range, because at t=0 the peak is a bit narrow and NIntegrate can miss it.

Let's check, how everything evolves with time:

Do[Print[NIntegrate[u[x, t] /. (First@sol) /. t -> tt, {x, -3, 3}]],
{tt, {1, 5, 10, 20, 30, 40}}]


100.08
100.08
100.08
100.08
100.08
100.08

Very consistent. But I don't like the 0.08 error. Let's try to make NDSolve's job a bit easier by specifying easier initial conditions (broader peak, so numerical precision isn't so much of an issue).

mass = 100;
width = 1;
g[x_] := mass /width HeavisideLambda[x/width]
sol = NDSolve[{pde, u[x, 0] == g[x], Derivative[1, 0][u][-3, t] == 0,
Derivative[1, 0][u][3, t] == 0}, u[x, t], {x, -3, 3}, {t, 0, 50}]
Do[Print[NIntegrate[u[x, t] /. (First@sol) /. t -> tt, {x, -3, 3}]],
{tt, {1, 5, 10, 20, 30, 40}}]


100.
100.
100.
100.
100.
100.

EDIT follow-up for ill-defined peak functions :-)

Clear[g]; deltaPos = 0;
g[x_] := Piecewise[{{0, x < deltaPos}, {100, x == deltaPos}, {0, x > deltaPos}}]
pde = D[u[x, t], t] - 0.2 D[D[u[x, t], {x, 1}], {x, 1}] == 0;
sol = NDSolve[{pde, u[x, 0] == g[x], Derivative[1, 0][u][-10, t] == 0,
Derivative[1, 0][u][10, t] == 0}, u[x, t], {x, -10, 10}, {t, 0, 50}]
Plot3D[u[x, t] /. sol, {x, -3, 3}, {t, 0.1, 3}, PlotRange -> Full]


produces this (below), as x=0 is included in the initial sampling points. Naturally, NDSolve complains about reaching the maximum of allowed grid points, because it detects a peak in g[x], but can't get a handle on how narrow it is:

NDSolve::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x.

Clear[g]; deltaPos = 2;
g[x_] := Piecewise[{{0, x < deltaPos}, {100, x == deltaPos}, {0, x > deltaPos}}]


With a definition like above NDSolve shuts up and returns u[x, t] that is zero across the domain. The initial sampling grid seems to be rather coarse. It simply doesn't realize, that there's an out-of-place peak in your starting function. deltaPos = 5 is, however, among the sampling points. That returns a much saner result, similar to what we get for deltaPos = 0`. The lesson to learn is "define your point sources properly".

Thanks for asking the question, it brought a great deal about solving ODEs with Mathematica to my attention, especially in the linked question in the comments.

• thanks for the useful comments and full explanation..! Apr 29 '15 at 18:02