# How to preserve normalization in NDSolve?

I have a probability density function: $$P_{init}(x)=\exp(-(x-x0)^2)/\sqrt{\pi}$$.

I am trying to use it as the initial condition for the following partial differential equation:

 Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]
V[x] = (-(x/5)^4)/Cosh[x/5];
F[x] = -D[V[x], x];
x0=5;
Pinit[x_] := Exp[-(x - x0)^2]/(Sqrt[Pi]);
T = 100;
BoundaryCondition = 250
uval = NDSolveValue[{D[u[x, t], t] + D[F[x]*u[x, t], x] -
D[u[x, t], x, x] == 0, u[x, 0] == Pinit[x],
u[-BoundaryCondition, t] == 0, u[BoundaryCondition, t] == 0},
u, {x, -BoundaryCondition, BoundaryCondition}, {t, 0, T}]


The above is a Fokker-Planck equation, which shows how the probability density expands in time.

The initial distribution is normalized, namely $$\int_{-\infty}^\infty {P_{init}(x)}dx=1$$, as it should.

However, it seems that no matter what T I choose, uval[x,T] never remains normalized.

Importantly: I get that uval[x,0] is different than Pinit(x), which is a contradiction.

How do I force Mathematica to solve the Fokker-Planck equation, whilst maintaining normalization?

Note that the reason that the integration boundaries are big, is since I would like to estimate the distribution at a long time, where the function might be much wider than the initial condition. This means that if I take boundaries which are too closely apart, I introduce mistakes because I force the function to be zero at a place and time where it shouldn't.

• People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful Commented Feb 7, 2019 at 19:34
• Commented Feb 7, 2019 at 19:40
• With such parameters T=100, -250<=x<=250 and initial data, the numerical solution cannot be sufficiently accurate. It is necessary to limit the area of integration within reasonable limits. Commented Feb 7, 2019 at 22:00
• Commented Feb 8, 2019 at 4:56
• @user1611107 I have just seen your edit to the question. I guess you are the same person as ForgotMyUserDetail? Please log in to your account. Then you can write comments to you own post and it would be much less confusing for other users. Moreover, you will be able to upvote helpful answers (this is what drives the community) and to mark the best answer as accepted. Commented Feb 8, 2019 at 9:33

OK, let me summarize my comments to an answer. This problem is related to this one. The i.c. is awry in uval because the peak in the i.c. is so narrow compared to the domain of definition that the default spatial grid is too coarse to capture it:

Clear[V, F]
V[x_] = (-(x/5)^4)/Cosh[x/5];
F[x_] = -D[V[x], x];
x0 = 5;
BoundaryCondition = 250;
Pinit[x_] = Exp[-(x - x0)^2]/(Sqrt[Pi]);
T = 100;

uval = NDSolveValue[{D[u[x, t], t] + D[F[x]*u[x, t], x] - D[u[x, t], x, x] == 0,
u[x, 0] == Pinit[x], u[-BoundaryCondition, t] == 0, u[BoundaryCondition, t] == 0},
u, {x, -BoundaryCondition, BoundaryCondition}, {t, 0, T}];

coordx = uval["Coordinates"][[1]]
(*
{-250.,    -229.167, -208.333, -187.5,   -166.667,
-145.833, -125.,    -104.167, -83.3333, -62.5,
-41.6667, -20.8333,  0.,       20.8333,  41.6667,
62.5,     83.3333,  104.167,  125.,     145.833,
166.667,  187.5,    208.333,  229.167,  250.}
*)

ptsx = Point[{#, 0} & /@ coordx]

Plot[{Pinit[x]}, {x, -BoundaryCondition, BoundaryCondition}, PlotRange -> All,
Epilog -> {Red, ptsx}]


Then the solution is simple: make the spatial grid dense enough to capture the peak and approximate it in an accurate enough way:

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

uvalfixed =
NDSolveValue[{D[u[x, t], t] + D[F[x]*u[x, t], x] - D[u[x, t], x, x] == 0,
u[x, 0] == Pinit[x], u[-BoundaryCondition, t] == 0, u[BoundaryCondition, t] == 0},
u, {x, -BoundaryCondition, BoundaryCondition}, {t, 0, T}, Method -> mol[2000, 4]];

Plot[{Pinit[x], uvalfixed[x, 0]}, {x, -BoundaryCondition, BoundaryCondition},
PlotRange -> All, PlotStyle -> {Automatic, {Thick, Dashed}}]


NIntegrate[uvalfixed[x, T], {x, -BoundaryCondition, BoundaryCondition}]

(* 1. *)

• @ xzczd: Great answer man, as far as I can see now, this really seems to solve the problem. Thanks! Commented Feb 10, 2019 at 11:49
• @ xzczd: Quick followup question, please? In your method, what would you do if the problem had a V(x) which diverges as x approaches +0, so that you want to place a reflecting boundary at some small positive x and solve numerically only in the half positive plane? Thanks! Commented Feb 10, 2019 at 20:51
• @user1611107 I'm sorry, but I don't understand what you mean… Commented Feb 11, 2019 at 4:10

As correctly pointed out in the comment the domain of integration is very large.

BoundaryCondition = 20;
Plot[{Pinit[x], uval[x, 0]}, {x, -BoundaryCondition, BoundaryCondition}, PlotRange -> All,
PlotStyle -> {Red, {Dashed, Green}}]


BoundaryCondition = 100;


BoundaryCondition = 150;


Let's try with MethodOfLines while keeping the original domain of interest,

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf : False | True, sf_: Automatic] := {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}

pts = 150;

uval = NDSolveValue[{D[u[x, t], t] + D[F[x]*u[x, t], x] -
D[u[x, t], x, x] == 0, u[x, 0] == Pinit[x],
u[-BoundaryCondition, t] == 0, u[BoundaryCondition, t] == 0},
u, {x, -BoundaryCondition, BoundaryCondition}, {t, 0, T},
Method -> Union[mol[pts, 6], mol[True, 100]]]

Plot[{Pinit[x], uval[x, 0]}, {x, -BoundaryCondition,
BoundaryCondition}, PlotRange -> All,
PlotStyle -> {Red, {Dashed, Green}}]


pts = 500;


It is still evident that the two doesn't match perfectly.

• It's not necessary to use a high difference order in this case, just use a dense enough grid, say, 500. Commented Feb 8, 2019 at 4:52
• @xzczd I agree. I tried many different combinations.
– zhk
Commented Feb 8, 2019 at 4:53
• Actually the most commonly used difference order is 2 or 4, and nowadays any average PC can bear with thousands of grid points. Commented Feb 8, 2019 at 5:15
• How would you solve it if D[u[x, t], t] + D[F[x]*u[x, t], x] - D[x*u[x, t], x, x] == 0? All I did is changed D[u[x, t], x, x] to D[x*u[x, t], x, x].
– sra
Commented Jan 31, 2023 at 21:38