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["DifferentialEquations`InterpolatingFunctionAnatomy`"]
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.
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. $\endgroup$ – Alex Trounev Feb 7 at 22:00