# Using NDSolveValue for Solving a parabolic PDE numerically

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: 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.

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

• Thank you very much - your solution is very useful. Jan 9, 2022 at 19:05
• You are welcome! Jan 9, 2022 at 20:49
• @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 Jan 10, 2022 at 4:09
• @xzczd Thank you for the hint. Please show your solution proposal using "TensorProductGrid", would be very helpful. Jan 10, 2022 at 7:07

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] Plot[tst[i, 0.5], {i, -3, 3}, PlotRange -> All] 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];
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!

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