3
$\begingroup$

I am trying to solve numerically the diffusion equation with absorbing boundary conditions on a finite domain. The initial condition is an approximation of a delta function centered in x0 with standard deviation sigma.

I find that changing the size of the domain slightly has a dramatic (non-physical) effect on the solution evaluated at a finite time. Why is this? Why don't I get a warning from Mathematica? This code reproduces the problem:

sigma = 1/32;
T = 10;
x0 = 1;
L = 24;
sol1 = NDSolve[{D[pn[x, t], {t, 1}] == 1/8 D[pn[x, t], {x, 2}], 
pn[0, t] == 0, pn[L, t] == 0, 
pn[x, 0] == 1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]}, 
pn[x, t], {x, 0, L}, {t, 0, T}];
L = 26;
sol2 = NDSolve[{D[pn[x, t], {t, 1}] == 1/8 D[pn[x, t], {x, 2}], 
pn[0, t] == 0, pn[L, t] == 0, 
pn[x, 0] == 1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]}, 
pn[x, t], {x, 0, L}, {t, 0, T}];
t = T;
GraphicsRow[{Plot[Evaluate[pn[x, t] /. sol1], {x, 0, 24}, 
PlotRange -> All], 
Plot[Evaluate[pn[x, t] /. sol2], {x, 0, 26}, PlotRange -> All]}]

On the left I plot the solution with domain size L=24, on the right I plot the solution with L=26. As you can see, the change is dramatic (notice the y axis values) and Mathematica produces no warning.

How can I solve this issue and be able to solve the equation numerically for L larger than 24?

enter image description here

$\endgroup$
  • 2
    $\begingroup$ Tip; Setting the variable of integration t to a numeric value with t = T messes things up when you (or I) recompute the solutions. It's very inconvenient. You can avoid needing such a workaround by solving for pn instead of pn[x, t] in your NDSolve calls. Alternatively, you could use the code pn[x, t] /. sol /. t -> T. $\endgroup$ – Michael E2 May 16 '17 at 21:00
3
$\begingroup$

You have control the quality of the grid/mesh of the spatial variable. Either with "MaxStepSize":

NDSolve[{D[ppn[x, t], {t, 1}] == 1/8 D[ppn[x, t], {x, 2}], 
  ppn[0, t] == 0, ppn[L, t] == 0, 
  ppn[x, 0] == 
   1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]}, ppn, {x, 0, 
  L}, {t, 0, T},
 Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid", 
     "MaxStepSize" -> 1}}]

or with "MinPoints":

NDSolve[{D[ppn[x, t], {t, 1}] == 1/8 D[ppn[x, t], {x, 2}], 
  ppn[0, t] == 0, ppn[L, t] == 0, 
  ppn[x, 0] == 
   1/Sqrt[2 Pi sigma^2] Exp[-(x - x0)^2/(2 sigma^2)]}, ppn, {x, 0, 
  L}, {t, 0, T},
 Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid", 
     "MinPoints" -> Ceiling@L}}]
$\endgroup$

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.