I am trying to solve a coupled diffusion partial differential equation with a Gaussian profile at t = 0. I would like to see how the two functions, y and z, evolve in time and in one dimension. This can be solved analytically using Fourier transform and I have been comparing my numerical solution to the analytic solution.

parms = {Ds -> 150, m -> 1, a -> -1, e -> 1, b -> 0.5};
xMax = 30000.;
tMax = 40000.;
nsol = NDSolve[{Ds D[y[x, t], x, x] - a/e b m D[z[x, t], 
     x] == D[y[x, t], t] /. parms, 
Ds D[z[x, t], x, x] - a/e b m D[y[x, t], x] == 
  D[z[x, t], t] /. parms, 
y[x, 0] == PDF[NormalDistribution[0, 1.], x], y[xMax, t] == 0., 
y[-xMax, t] == 0., z[x, 0] == 0, z[xMax, t] == 0., 
z[-xMax, t] == 0.}, {y[x, t], z[x, t]}, {x, -xMax, xMax}, {t, 0, 
tMax}, Method -> {"MethodOfLines", "TemporalVariable" -> t, 
 "SpatialDiscretization" -> {"TensorProductGrid", 
   "MaxPoints" -> 625, "MinPoints" -> 625}}];

Implementing the above gives good agreement with analytic result to within a scaling factor. I've tried several other "Methods" with varying degrees of success. I have the following questions.

  1. Why does the numerical solution change scale as I change xMax?

  2. MaxPoints and MinPoints are lucky choices that give a solid result. However if I change those numbers to 626, each solution goes to zero. How do I know what to choose for Max/Min Points?

  3. When I solve analytically I don't need to specify boundary conditions at xMax or tMax - I only enforce some condition at t = 0 and that suffices. Is there a way to do this here with NDSolve where I can reduce my number of boundary conditions? In the same vein, would a another "Method" be better for my equations?

  • 1
    $\begingroup$ If you did not explicitly introduce boundary conditions in your analytical solution, you may have implicitly assumed periodicity in x. If so, you may wish to run your code with periodic boundary conditions for a more rigorous comparison. $\endgroup$
    – bbgodfrey
    Aug 7, 2017 at 23:25
  • $\begingroup$ The code produces the same results for "MaxPoints" -> 625, "MinPoints" -> 625 and for 623 and 627. However, nearby even values yield essentially zero as a solution. Likewise a periodic solution is essentially the same as the zero boundary condition case, except at late time, when the propagating waves in the solution collide in the periodic solution $\endgroup$
    – bbgodfrey
    Aug 8, 2017 at 0:44
  • 1
    $\begingroup$ @bbgodfrey As to the implicit b.c. introduced when Fourier transform is used for finding analytic solution, you may have a look at this post. $\endgroup$
    – xzczd
    Aug 8, 2017 at 10:46

1 Answer 1


First, I'd like to point out, your statement

When I solve analytically (using Fourier transform) I don't need to specify boundary conditions at xMax or tMax - I only enforce some condition at t = 0 and that suffices…

is not true. When finding analytic solution with Fourier transform, you've actually impose b.c. implicitly. See this post for more information.

Currently methods implmemented in NDSolve for solving PDEs only handle problem defined in a finite domain (well, actually I'm not aware of any numerical methods that can handle PDEs defined in infinite domain directly), so, when dealing with problem defined in an infinite domain, artifical b.c. is necessary to approximate b.c. at infinity, which can be really troublesome in some cases. (See here for example. )

However, AFAIK, when dealing with heat equation and its close relatives, b.c.s (Dirichlet b.c. or Neumann b.c.) at far enough places are usually good enough approximation for b.c.s at infinity i.e. you should have obtain a numeric solution that is in not bad agreement with the analytic one while you didn't, why?

Because your xMax is way too large.

Just plot the i.c. and you'll see how unnecessarily large it is:

xMax = 30000.;
Plot[PDF[NormalDistribution[0, 1.], x] // Evaluate, {x, -xMax, xMax}, Axes -> None, 
 PlotPoints -> 101]

Mathematica graphics

Making the domain unnecessarily large is a waste of computing resource, yet it's not the biggest problem. The biggest problem is, to describe the i.c. in such a domain with a uniform grid accurately, a very dense uniform spatial grid is needed, which generally overwhelms the RAM of an average computer. (One may argue that we can set a non-uniform grid with the "Coordinate" option, but again AFAIK, "Coordinate" is a subtle option, when one sets it arbitrarily, NDSolve may fail easily so I'd like not to touch it in this case.)

The grid point number 625 is undoubtedly too sparse, luckily it's an odd number, so one of the grid point hits the peak of the i.c., that's why you obtain a result that "gives good agreement with analytic result to within a scaling factor", but when the grid point number is an even number, every grid point hits a position where the i.c. almost equals 0, that's the reason why you got 0 as solution when changing spatial grid points to 626.

Finally let me fix your code. As mentioned above, we just need to make xMax a bit smaller. (BTW I've modified the b.c. to Neumann b.c. because it turns out to be a bit better):

parms = {Ds -> 150, m -> 1, a -> -1, e -> 1, b -> 0.5};

xMax = 30000./1000;
tMax = 40000.;

eq = {Ds D[y[x, t], x, x] - a/e b m D[z[x, t], x] == D[y[x, t], t], 
   Ds D[z[x, t], x, x] - a/e b m D[y[x, t], x] == D[z[x, t], t]};
ic = {y[x, 0] == PDF[NormalDistribution[0, 1], x], z[x, 0] == 0};

{nysol, nzsol} = 
  NDSolveValue[{eq, ic, 
     D[#[x, t] == 0, x] & /@ {y, z} /. {{x -> -xMax}, {x -> xMax}}} /. parms, {y, 
    z}, {x, -xMax, xMax}, {t, 0, tMax}];

Compare it to the analytic solution found by Fourier transform (Definition of ft can be found here):

(*Definition of ft isn't included in this post,
  please find it in the link above.  *)
tset = ft[{eq, ic}, x, w] /. HoldPattern@FourierTransform[a_, __] :> a

tsol = DSolve[tset, {y[x, t], z[x, t]}, t]

{ysol[x_, t_], zsol[x_, t_]} = 
 InverseFourierTransform[tsol[[1, All, -1]], w, x] // Simplify

Manipulate[Plot[{nysol[x, t], ysol[x, t]} /. parms // TrigToExp // Evaluate, 
   {x, -xMax, xMax}, PlotRange -> {0, 0.5}, PlotStyle -> {{Dashed, Thick}, Red}, 
  PlotLegends -> {"Numeric", "Analytic"}], {t, 0, 1}]

enter image description here

As one can see, the numerical solution deviates a bit when t gets larger, you can make xMax a bit larger as a countermeasure, but once again, remember not to make it too large.

  • 1
    $\begingroup$ Thanks, this is all very helpful. About the range of xMax...look at Ds = 0.1 where the drift part dominates the diffusion part (I should have posted with this value instead of 150). You see packets now traveling left and right and they take long to diffuse. That being said I can still reduce xMax depending on which parameters I'm choosing. $\endgroup$
    – BeauGeste
    Aug 8, 2017 at 16:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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