# unstable results when solve PDE even with great artificial diffusion term added

I try to solve a group of coupled PDEs, the 1-D continuity equation and the dynamic function:

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

$$\frac {\partial W}{\partial t}=-\frac {1}{m n}\frac {\partial (n k_B T)}{\partial z}$$

It can be seen that, from the second dynamic equation, the velocity $$W$$ would tend to carry the particles from large $$n$$ to small $$n$$, so it would smooth the particles denstity gradient $$\frac {\partial n}{\partial z}$$. Once the denstity gradient disappear, the force drives $$W$$ should die out. So it seems these two equations are self-stabled.

However, when I try to solve these equations

kB = 1.380649*10^(-23);
mp = 1.6726231*10^(-27);
NDSolve[
{D[ni[t, z], t] == -D[W[t, z]*ni[t, z],z],
D[W[t, z], t] == -1/(16 mp*ni[t, z])*(D[ni[t, z]*kB*3000*1.5, z]),
{ni[0, z] == 100*Exp[-z/100], ni[t, 200] == 100*Exp[-200/100],
ni[t, 1000] == 100*Exp[-1000/100], W[0, z] == 0, W[t, 200] == 0,
W[t, 1000] == 0}}, {ni, W}, {t, 0, 100}, {z, 200, 1000}]


The Mathematica returns warning message and a solution stops at t=0.756. The density ni[t,z] at t=0.7 reach a ridiculously large value (~10^15).

I try to stablize the equation by adding really large artificial diffusion term 100000 D[ni[t, z], z, z]:

NDSolve[
{D[ni[t, z], t] - 100000 D[ni[t, z], z, z] == -D[W[t, z]*ni[t, z], z],
D[W[t, z],t] == -1/(16 mp*ni[t, z])*(D[ni[t, z]*kB*3000*1.5, z]),
{ni[0, z] == 100*Exp[-z/100], ni[t, 200] == 100*Exp[-200/100],
ni[t, 1000] == 100*Exp[-1000/100], W[0, z] == 0, W[t, 200] == 0,
W[t, 1000] == 0}}, {ni, W}, {t, 0, 100}, {z, 200, 1000}]


However, this only elongate the solution to t=2.1, ni[t,z] runs away again. What should I do?

• The scaling of your problem seems to be an issue! Oct 30, 2019 at 12:34
• According to Tim's answer below, the i.c.s and b.c.s seem to be inconsistent. For this type of PDE inconsistent i.c. and b.c. is troublesome. Example: mathematica.stackexchange.com/q/145416/1871 Oct 30, 2019 at 13:54
• @xzczd I think that you are probably right about the inconsistency. I updated my answer to zoom in on the boundary areas. As the model stands, it looks like there are some very sharp features to capture. Oct 30, 2019 at 14:50
• @Harry I added a pseudospectral difference order to the NDSolve options. Instabilities can occur, but they don't seem to blow up without bound as they do when the option is not specified. Oct 31, 2019 at 2:45

# Update: Stability Enhancement Using Pseudospectral

Some experimentation shows that a Pseudospectral difference order can prevent the solution instabilities from blowing up and thereby allowing for smaller artificial diffusion.

opts = (Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 1000, "DifferenceOrder" -> "Pseudospectral"}});
kB = 1.380649*10^(-23);
mp = 1.6726231*10^(-27);
pfun = ParametricNDSolveValue[{D[ni[t, z], t] +
D[-c1 D[ni[t, z], z] - (-W[t, z]*ni[t, z]), z] ==
0, (16 mp)/(kB*3000*1.5)*ni[t, z] D[W[t, z], t] +
D[-c2 D[W[t, z], z] - (-ni[t, z]), z] == 0,
ni[0, z] == 100*Exp[-z/100], ni[t, 200] == 100*Exp[-200/100],
ni[t, 1000] == 100*Exp[-1000/100], W[0, z] == 0, W[t, 200] == 0,
W[t, 1000] == 0}, ni, {t, 0, 1}, {z, 200, 1000}, {c1, c2}, opts,
MaxStepFraction -> 1/100];
ifun = pfun[1/10, 1/10];
imgs = Plot[Evaluate[ifun[#, z]], {z, 200, 1000},
PlotRange -> {0, All}] & /@ Subdivide[0, 1, 100];
ListAnimate@imgs


The solution really is starting to look like wave. Perhaps there is a way to recast as such and take advantage of the techniques to stabilize wave simulations.

You could try adding artificial diffusion to both $$n_i$$ and $$W$$ and use a finer grid with a higher difference order. For example:

opts = (Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 4000, "DifferenceOrder" -> 10}});
kB = 1.380649*10^(-23);
mp = 1.6726231*10^(-27);
pfun = ParametricNDSolveValue[{D[ni[t, z], t] +
D[-c1 D[ni[t, z], z] - (-W[t, z]*ni[t, z]), z] ==
0, (16 mp)/(kB*3000*1.5)*ni[t, z] D[W[t, z], t] +
D[-c2 D[W[t, z], z] - (-ni[t, z]), z] == 0,
ni[0, z] == 100*Exp[-z/100], ni[t, 200] == 100*Exp[-200/100],
ni[t, 1000] == 100*Exp[-1000/100], W[0, z] == 0, W[t, 200] == 0,
W[t, 1000] == 0}, ni, {t, 0, 100}, {z, 200, 1000}, {c1, c2}, opts,
MaxStepFraction -> 1/100]
ifun = pfun[1, 1];
imgs = Table[
Plot[Evaluate[ifun[t, z]], {z, 200, 1000},
PlotRange -> {0, All}], {t, 0, 100}];
ListAnimate@imgs


You will have to determine if the diffusive errors are acceptable, but it does stabilize the solution.

I agree with xzczd's comment that the BC's and IC's look inconsistent. If they are not, you can see that there are some extremely sharp features to capture in both time and space by zooming in to the solution.

# Adding Growth Factor Refinement to Both Ends of the Domain

As stated previously, there are some fine features to resolve at both ends of the domain. To resolve these features more economically than a uniform mesh, you could consider a non-uniform mesh with geometric growth rate on the ends. Here is one way to accomplish the refinement.

meshGrowth[x0_, xf_, n_, ratio_] := Module[{k, fac, delta},
k = Log[ratio]/(n - 1);
fac = Exp[k];
delta = (xf - x0)/Sum[fac^(i - 1), {i, 1, n - 1}];
N[{x0}~Join~(x0 +
delta Rest@
FoldList[(#1 + #2) &, 0,
PowerRange[fac^0, fac^(n - 3), fac]])~Join~{xf}]
]
unitMeshGrowth[n_, ratio_] := meshGrowth[1, 0, n, ratio]
unitMeshGrowth2Sided [nhalf_,
ratio_] := (1 + Union[-Reverse@#, #])/2 &@
unitMeshGrowth[nhalf, ratio]
grid = 200 + 800*unitMeshGrowth2Sided[500, 100000];
ListPlot@grid
opts = (Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"Coordinates" -> {grid}, "DifferenceOrder" -> 10}});
kB = 1.380649*10^(-23);
mp = 1.6726231*10^(-27);
pfun = ParametricNDSolveValue[{D[ni[t, z], t] +
D[-c1 D[ni[t, z], z] - (-W[t, z]*ni[t, z]), z] ==
0, (16 mp)/(kB*3000*1.5)*ni[t, z] D[W[t, z], t] +
D[-c2 D[W[t, z], z] - (-ni[t, z]), z] == 0,
ni[0, z] == 100*Exp[-z/100], ni[t, 200] == 100*Exp[-200/100],
ni[t, 1000] == 100*Exp[-1000/100], W[0, z] == 0, W[t, 200] == 0,
W[t, 1000] == 0}, ni, {t, 0, 100}, {z, 200, 1000}, {c1, c2}, opts,
MaxStepFraction -> 1/100]
ifun = pfun[1/2, 1/2];
imgs = Plot[Evaluate[ifun[#, z]], {z, 200, 1000},
PlotRange -> {0, All}] & /@ Subdivide[0, 1, 100];
ListAnimate@imgs


By adding the refinement, I was able to reduce the grid by 4x and the artificial diffusion by 2x.

Although the solution my not be valid, stabilizing it through artificial diffusion may give you clues as to where the model breaking down.

• Oh...I understand. Seems that the unequal boundarys ni[t, 200] and ni[t, 1000] are guilty... Oct 31, 2019 at 2:48