I use NDSolve to solve a 1-D continuity equation:

$$\frac{\partial n}{\partial t}+\frac{\partial(W n)}{\partial z}=0$$

When Integrated over the all height (z), it leads to

$$\frac \partial {\partial t} \int n dz =-W(z_\text{max})n(z_\text{max}) + W(z_\text{min})n(z_\text{min})$$

If $W(z_\text{max})=W(z_\text{min})=0$.

Then $\int n dz = \text{const}$, so it is called "continuity equation". Here I chose $W(z)=1×\sin\left(\frac{2×6 \pi(z - 200)}{1000 - 200}\right)$ and solve the eqution from $z=200$ to $z=1000$. It can be seen that $W$ equals $0$ at boundary so the total "particle number" $n$ (ni) should be conserved. Noted that this W is a simplified example. The real $W$ I have to handle is much more complicated and evolved with time $t$.

fun = NDSolve[{D[ni[t, z], t] == -D[1*Sin[2*6 Pi (z - 200)/(1000 - 200)]*ni[t, z], z], {ni[
  0, z] == 100, ni[t, 200] == 100, ni[t, 1000] == 100}}, {ni}, {t,
0, 100}, {z, 200, 1000}]

However, after solving this equation by NDSolve, I found that

((ni /. First@fun)[0, #] & /@ Range[200, 1000]) // Total

gives 80100. and

((ni /. First@fun)[20, #] & /@ Range[200, 1000]) // Total

gives 78698. and

((ni /. First@fun)[50, #] & /@ Range[200, 1000]) // Total

gives 19582.5.

It seems $\int n dz$ decreases as time evolved. The solution I get can't be correct. How could I get the right solution?


2 Answers 2


Making the spatial grid denser helps:

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

fun = NDSolve[{D[ni[t, z], 
     t] == -D[1*Sin[2*6 Pi (z - 200)/(1000 - 200)]*ni[t, z], z], {ni[0, z] == 100, 
    ni[t, 200] == 100}}, {ni}, {t, 0, 100}, {z, 200, 1000}, 
  Method -> mol[10^5, 2]]
NIntegrate[(ni /. First@fun)[50, z], {z, 200, 1000}, 
 Method -> {Automatic, SymbolicProcessing -> 0}, MaxRecursion -> 20]
(* 80001. *)

The solution is rather sharp for large t, that's the reason why a extremely dense grid is needed:

ListLinePlot[(ni /. fun[[1]])["ValuesOnGrid"] // Last, PlotRange -> All, 
 DataRange -> {200, 1000}]

enter image description here

  • $\begingroup$ Thanks a lot! This solved my problem. And a little question: shouldn't the peaks in the final solution share same value? Since each peak concentrates the paticles 100m around it. $\endgroup$
    – Harry
    Oct 25, 2019 at 12:29
  • $\begingroup$ @Harry Seems that I'm too tired… In this case it's better to use ListPlot for visualizing, and the peak height at $t=100$ for 10^5 is already (almost) the same. See my update. $\endgroup$
    – xzczd
    Oct 25, 2019 at 13:10


This update is a cleaner implementation of adding artificial diffusion on a finer mesh.

First, for consistency, it is good to pose the PDE in coefficient form as described here.

$$m\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

Not only will this provide one with a more consistent workflow within the Wolfram Language, but you will be able to more easily map the problem to other solvers.

In the original answer, I applied artificial diffusion in a quick and dirty way to a system that had much numerical diffusion due to its coarse mesh. Now, I will apply artificial diffusion using coefficient form on a much finer mesh with a high difference order. This should allow us to use a much smaller artificial diffusion, $c$. The following solves and plots in a few seconds on my system.

opts = (Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 2000, "DifferenceOrder" -> 10}});
pfun = ParametricNDSolveValue[{D[ni[t, z], t] + 
     D[-c D[ni[t, z], 
         z] - (-Sin[2*6 Pi (z - 200)/(1000 - 200)] ni[t, z]), z] == 
    0, {ni[0, z] == 100, ni[t, 200] == 100, ni[t, 1000] == 100}}, 
  ni, {t, 0, 100}, {z, 200, 1000}, {c}, opts];
Plot[Evaluate[pfun[1/50][100, z]], {z, 200, 1000}, PlotRange -> All] 

New Artificial Diffusion

The solution looks stable throughout the entire time range.

imgs = Table[
   Plot[Evaluate[pfun[1/50][t, z]], {z, 200, 1000}, 
    PlotRange -> {0, 5500}], {t, 0, 100}];

Peaks Animation

Using the conservative coefficient form of the PDE, it is easy to map to another solver, such as COMSOL, to show that the solutions are similar. COMSOL comparison

Now, this answer is not as resolved nor accurate as @xcczd's, but it can help stabilize unstable systems. Studying the solutions of stabilized systems may give insights to future refinements. One can consume a lot time fighting instability.

Original Answer

If you can tolerate adding some artificial diffusion as described in the Stabilization of Convection-Dominated Equations Section of FEM Best Practices to stabilize your solution, that could speed your computations long term.

We can see that your initial approach is numerically unstable.

ifun = ni /. 
  First@NDSolve[{D[ni[t, z], 
       t] == -D[1*Sin[2*6 Pi (z - 200)/(1000 - 200)]*ni[t, z], 
        z], {ni[0, z] == 100, ni[t, 200] == 100, 
      ni[t, 1000] == 100}}, {ni}, {t, 0, 100}, {z, 200, 1000}]
Plot3D[ifun[t, z], {t, 0, 100}, {z, 200, 1000}, PlotRange -> All]

Unstable solution

We can use ParametricNDSolveValue to find the minimum artificial diffusion needed to stabilize the solution. An artificial diffusion coefficient of 200 seems to work pretty well.

pfun = ParametricNDSolveValue[{D[ni[t, z], t] - 
     dart D[ni[t, z], z, z] == -D[
      1*Sin[2*6 Pi (z - 200)/(1000 - 200)]*ni[t, z], z], {ni[0, z] == 
     100, ni[t, 200] == 100, ni[t, 1000] == 100}}, 
  ni, {t, 0, 100}, {z, 200, 1000}, {dart}, StartingStepSize -> 0.1]
 Plot3D[pfun[d][t, z], {t, 0, 100}, {z, 200, 1000}, PlotPoints -> 50, 
  PlotRange -> All], {d, 0, 400, Appearance -> "Labeled"}]

Artificial Diffusion

We can see that particles are conserved to about 0.2% over the time range.

pts = Table[
   NIntegrate[pfun[200][t, z], {z, 200, 1000}, 
    Method -> {Automatic, SymbolicProcessing -> 0}], {t, 0, 100}];

Conservation Plot

  • $\begingroup$ Emmm...Thanks for your enlightening anwser. I wonder that once diffusion term involved, the particles are no longer expected to be conserved, since they can diffuse in or out through boundary...And It seems that the diffusion term changes the solution signifantly when compared with the answer posted by xzczd. $\endgroup$
    – Harry
    Oct 25, 2019 at 15:26
  • $\begingroup$ The initial answer applied artificial diffusion to a coarse mesh, which inherently has a lot of numerical diffusion. I used a finer discretization so that much less artificial diffusion can be used. It depends on what your needs are, but if you find yourself spending a lot of time fighting stability issues, you can use this approach to stabilize the numerics understanding that you introduced error. The stabilized solution can give you insights for further refinements. $\endgroup$
    – Tim Laska
    Oct 26, 2019 at 13:02

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.