4
$\begingroup$

Inspired by this question I am trying to solve the following PDE numerically on $x \in [-3, 3]$ and $t \in [0, 0.5]$ using NDSolveValue:

$$ \frac{\partial p}{\partial t} = (12x^2-4) p + \left[4x(x^2-1)+0.1\right] \frac{\partial p}{\partial x} + \frac{1}{\beta} \frac{\partial^2 p}{\partial x^2} \\ $$

The initial condition is $p(x,0)=\frac{1}{\sqrt{2\pi}} \exp(-\frac{x^2}{2})$ and the zero boundary condition is $p(-3,t)=p(3,t)=0$.

As a first orientation I looked at the approach described by this post entitled "Solving 2D+1 PDE with Pseudospectral in one direction with periodic boundary condition?" and used:

b = 200;
pval = NDSolveValue[{D[p[x,t],t] == (12x^2-4)p[x,t]+(4x(x^2-1)+0.1) *
   D[p[x,t],x] + D[p[x,t]/b, {x,2}], p[x,0] == Exp[-x^2/2]/Sqrt[2Pi],
       p[-3,t] == 0, p[3,t] == 0}, p, {x,-3,3}, {t,0,0.5}];

ContourPlot[pval[i,j], {i,-3,3}, {j,0,0.5}]

After choosing $b=-200$, $b=-1$, $b=1$ and $b=200$, it generated the following Contour Plots:

enter image description here

My problem/question:

Unfortunatelly the zero boundary condition doesn't seem to hold, since we should have $p(−3,t)=p(3,t)=0$ for any $t$. Something I am doing wrong and it would be great if anyone could help me.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
8
$\begingroup$

MethodOfLines gives a solution to your problem:

b = 200;
P = NDSolveValue[{D[p[x, t],t] == (12 x^2 - 4) p[x, t] + (4 x (x^2 - 1) + 0.1)*D[p[x, t], x] + D[p[x, t]/b, {x, 2}], 
p[x, 0] == Exp[-x^2/2]/Sqrt[2 Pi], p[-3, t] == 0, p[3, t] == 0}, 
p, {x, -3, 3}, {t, 0, 0.5}, 
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement","MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}];

Plot3D[P[x, t], {x, -3, 3}, {t, 0, .5}]

enter image description here

The boundary conditions aren't fullfilled exactly, because boundary and initial conditions are inconsistent.

inconsistence might be removed if you change the initial conditions to p[x, 0] == Which[0 < RealAbs[x] <= 3, Min[(3 - RealAbs[x])/.01, 1], True, 0] Exp[-x^2/2]/Sqrt[2 Pi]

Hope it helps!

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you very much - your solution is very useful. $\endgroup$
    – Eldar Sultanow
    Jan 9, 2022 at 19:05
  • $\begingroup$ You are welcome! $\endgroup$
    – Ulrich Neumann
    Jan 9, 2022 at 20:49
  • 1
    $\begingroup$ @EldarSultanow and Ulrich: The critical part is not MethodOfLines (Eldar's code also uses it under the hood, and this is the only available method built in NDSolve for time dependent problem at the moment, but the sub-option is TensorProduct by default in this case, related: mathematica.stackexchange.com/q/140773/1871 ), but FiniteElement. Related: mathematica.stackexchange.com/q/230429/1871 $\endgroup$
    – xzczd
    Jan 10, 2022 at 4:09
  • $\begingroup$ @xzczd Thank you for the hint. Please show your solution proposal using "TensorProductGrid", would be very helpful. $\endgroup$
    – Ulrich Neumann
    Jan 10, 2022 at 7:07
8
$\begingroup$

OK, let me extend my comment to an answer. First of all, I'd like to point out that OP's observation

the zero boundary condition doesn't seem to hold

isn't really a problem here. Though ibcinc warning does pop up, and the zero boundary condition is not exactly hold, the error is rather small in this case. For more info you may refer to the this post. The real problem here, is the eerr warning pops up and the corresponding error is so large (it's 280235.11872106727! ) that we cannot ignore it.

Then why does the estimated error so large? In short, this problem is quite similar with this one i.e. naive spatial discretization won't work for the PDE and we need to respect the conservation law in some way when discretizing. The simplest solution is to, as shown by Ulrich in his answer, turn to FiniteElement method, but still, we can use TensorProductGrid method. We just need to transform the PDE a little and set a low DifferenceOrder:

With[{p = p[x, t], mid = mid[x, t]},
 eq = D[p, t] == D[mid, x];
 eqadd = mid == (4 x (x^2 - 1) + 0.1) p + D[p/b, x];]
ic = p[x, 0] == Exp[-x^2/2]/Sqrt[2 Pi];
icadd = mid[x, 0] == eqadd[[-1]] /. p -> Function[{x, t}, Evaluate@ic[[-1]]];
bc = {p[-3, t] == 0, p[3, t] == 0};

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

tst = NDSolveValue[{eq, eqadd, ic, icadd, bc}, p, {x, -3, 3}, {t, 0, 0.5}, 
    Method -> mol[300, 2]]; // AbsoluteTiming

The eerr warning still pops up, but the corresponding error (815.9 ) isn't too large. Let's check how it looks like:

Plot3D[tst[i, j], {i, -3, 3}, {j, 0, 0.5}, PlotRange -> All, PlotPoints -> 200]

Mathematica graphics

Plot[tst[i, 0.5], {i, -3, 3}, PlotRange -> All]

Mathematica graphics

Looks reasonable, and it's consistent with that of FEM. Further check shows, the solution doesn't change significantly even if we make the grid denser, so it's probably reliable.

Update

Aha, the magic of Pseudospectral seems to be effective for this problem:

b = 200;
pval = NDSolveValue[{D[p[x, t], 
      t] == (12 x^2 - 4) p[x, t] + (4 x (x^2 - 1) + 0.1) D[p[x, t], x] + 
      D[p[x, t]/b, {x, 2}], p[x, 0] == Exp[-x^2/2]/Sqrt[2 Pi], p[-3, t] == 0, 
    p[3, t] == 0}, p, {x, -3, 3}, {t, 0, 0.5}, Method -> mol[401]];
Plot[pval[i, 0.5], {i, -3, 3}, PlotRange -> All]

It's not necessary to modify the form of PDE even if TensorProductGrid is chosen!

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you very much for your extended comment. PseudoSpectral has "saved" the method TensorProductGrid( ;-) ). But FiniteElement seems to be significantly faster (factor 6) $\endgroup$
    – Ulrich Neumann
    Jan 10, 2022 at 15:47
  • $\begingroup$ @UlrichNeumann Which version are you in? In v12.3.1, win10, the timing of mol[401] and "MaxCellMeasure" -> 0.01 is 0.736 seconds v.s. 5.86 seconds, and precision of FEM is worse near the boundary. $\endgroup$
    – xzczd
    Jan 11, 2022 at 1:02
  • $\begingroup$ My version is v12.2 Win10. $\endgroup$
    – Ulrich Neumann
    Jan 11, 2022 at 7:00

Your Answer

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