7
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Sjoerd C. de Vries Apr 29 '15 at 11:20
  • $\begingroup$ @SjoerdC.deVries, I know in c++, this way of defining point source works just fine. I am not sure with Mathematica. $\endgroup$ – alekhine Apr 29 '15 at 11:23
  • 1
    $\begingroup$ Related: mathematica.stackexchange.com/a/71945/1871 $\endgroup$ – xzczd Apr 29 '15 at 11:53
  • 2
    $\begingroup$ 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. $\endgroup$ – LLlAMnYP Apr 29 '15 at 15:02
  • 2
    $\begingroup$ @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. $\endgroup$ – LLlAMnYP Apr 29 '15 at 16:27
8
$\begingroup$

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

Now back to your code:

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]

Initial conditions

Huge peak at the beginning (just a finite approximation of a delta-function).

Plot3D[u[x, t] /. sol, {x, -10, 10}, {t, 0, 20}]

Development

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

Levelling off

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.

Lucky us

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.

$\endgroup$
  • $\begingroup$ thanks for the useful comments and full explanation..! $\endgroup$ – alekhine Apr 29 '15 at 18:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.