1
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • $\begingroup$ The scaling of your problem seems to be an issue! $\endgroup$ Commented Oct 30, 2019 at 12:34
  • 1
    $\begingroup$ 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 $\endgroup$ Commented Oct 30, 2019 at 13:54
  • $\begingroup$ @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. $\endgroup$ Commented Oct 30, 2019 at 14:50
  • $\begingroup$ @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. $\endgroup$ Commented Oct 31, 2019 at 2:45

1 Answer 1

Reset to default
4
$\begingroup$

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

Pseudospectral solution

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.

Original Answer

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

Coupled problem with diffusion on both terms

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

Additional Observations

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.

Zoom of Time from 0 - 1 s.

Time Zoom

Space and Time Zoom

Space and Time Zoom

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

Mesh refined solution

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.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Oh...I understand. Seems that the unequal boundarys ni[t, 200] and ni[t, 1000] are guilty... $\endgroup$ Commented Oct 31, 2019 at 2:48

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.