unknown cause of scaling factors in NDSolve solutions to a partial differential equation

Question

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?

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]

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

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.