PDE: solving Fokker-Planck equation with initial and boundary condition

Question

Here is the problem. We have the following simple PDE: \begin{equation} \frac{\partial p(x,t)}{\partial t}= - a\frac{\partial p(x,t)}{\partial x} + \frac{D}{2} \frac{ \partial^2 p(x,t) }{\partial x^2}, \quad0<x<L, \quad t>0 \end{equation} where $a$ and $D$ are constants. We also have the following initial and boundary conditions: \begin{equation} p(x,0)=0\\ a p(x,t)- \frac{D}{2} \frac{\partial p(x,t)}{\partial x} = f(t) \quad at \quad x=0 \\ \frac{\partial^2 p(x,t)}{\partial x^2}=0 \quad at \quad x=L \end{equation}

The expression, $\frac{\partial^2 p(x,t)}{\partial x^2}=0$, is not a standard well-known boundary condition and makes the PDE difficult to solve. I would be grateful if anyone have any idea for solving this PDE.

Answer 1

The equations in the question can be written as

{D[p[x, t], t] + a D[p[x, t], x] - d/2 D[p[x, t], {x, 2}] == 0, 
 p[x, 0] == 0, 
 a p[0, t] - d/2 (D[p[x, t], x] /. x -> 0) - Sin[t] == 0, 
 D[p[xm, t], t] + a (D[p[x, t], x] /. x -> xm) == 0};

Unfortunately,

NDSolve[%, p, {t, 0, 1}, {x, 0, xm}]

returns unevaluated with the error message,

NDSolveValue::bdord: Boundary condition (p^(0,1))[1,t]+(p^(1,0))[1,t] should have derivatives of order lower than the differential order of the partial differential equation. >>

because the last boundary condition, above, contains a first derivative in time. This situation has much in common with question 78493, even though the equations are different. Following that approach with (for example) f = Sin[t], we decompose the PDE into a set of ODEs at n locations in x.

d = 1; a = 1; xm = 1;
n = 50; h = xm/n;
U[t_] = Table[u[i][t], {i, 0, n}];
eqns = Thread[D[U[t], t] == Join[
   {u[0][t] - Sin[t] - d/2 (u[1][t] - u[0][t])/h}, 
   -a ListCorrelate[{1, 0, -1}/(2 h), U[t]] + d/2 ListCorrelate[{1, -2, 1}/h^2, U[t]], 
   {-a (u[n][t] - u[n - 1][t])/h}]];
eqns = eqns /. eqns[[1, 1]] -> 0;
init = Delete[Thread[U[0] == 0], 1];
lines = NDSolve[{eqns, init}, U[t], {t, 0, 4 Pi}];
ParametricPlot3D[Evaluate[Table[{i h, t, First[u[i][t] /. lines]}, {i, 0, n}]],
  {t, 0, 4 Pi}, PlotRange -> All, AxesLabel -> {"x", "t", "u"}, BoxRatios -> {1, 4, 1}]

Note that the array of equations, eqns, contains the boundary conditions,

eqns[[1]]
eqns[[n + 1]]
(* 0 == -Sin[t] + u[0][t] - 25 (-u[0][t] + u[1][t]) *)
(* Derivative[1][u[50]][t] == -50 (-u[49][t] + u[50][t]) *)

as well as n - 1 equations representing the PDE itself, for instance,

eqns[[n/2]]
(* Derivative[1][u[24]][t] == -25 u[23][t] + 25 u[25][t] + 
   1/2 (2500 u[23][t] - 5000 u[24][t] + 2500 u[25][t]) *)

That such equations are readily solvable numerically suggests (to me, at least) that this capability should be added to NDSolve.