在初始时刻将一个波包放在rbs处考察波包随时间的演化
r1 = 0.05;
r2 = 5;
xx1 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r1;
xx2 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r2;
GLb = (1 + 45 rb^8 + 24 PEb rb^10 - 64 rb^9 TEb)/(
6 rb^4); rbs = 0.6804971882808477`;
rbl = 2.780560197330779`; PEb = 0.2379247924446251`; TEb = 0.5`;
pdeSoln =
NDSolveValue[{D[\[Rho][t, rb], t] ==
500 TEb D[E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb],
rb], \[Rho][0, rb] ==
1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbs)^2/0.01)), xx1, xx2}, \[Rho][t,
rb], {t, 0, 20000}, {rb, r1, r2}]
Plot3D[pdeSoln, {t, 0, 20000}, {rb, r1, r2}, PlotRange -> {0, 6},
AspectRatio -> 10/13, PlotPoints -> 150, MaxRecursion -> 2,
ColorFunction -> "Rainbow"]
给出了一个有问题的结果

GLb图如下所示，rbs和rbl是对应两个GLb极小值

如果把波包放在rbl处，即\[Rho][0, rb] == 1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbl)^2/0.01))
就会得到好结果，请问问题出在哪

方程是Fokker planck方程，ρ是演化波包，波包有一个振幅，我选的振幅是1/（0.1√π），即小于6，但结果演化到后面发现波包振幅随着时间的增长都越来越大

1.我所期待的结果是3楼的结果两个峰交换位置，演化时间可能不一样。
2.方程是Fokker Planck方程，两个边界方程，一个初始波包。GLb是势函数。
3.时间区间我不确定。
4.3楼结果已显示

3楼代码是
r1 = 0.05;
r2 = 5;
xx1 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r1;
xx2 = (D[\[Rho][t, rb], rb]) == -(1/TEb) (D[GLb, rb]) \[Rho][t,
rb] /. rb -> r2;
pdeSoln =
NDSolveValue[{D[\[Rho][t, rb], t] ==
500*TEb D[E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb],
rb], \[Rho][0, rb] ==
1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbl)^2/0.01)), xx1, xx2}, \[Rho][t,
rb], {t, 0, 20000}, {rb, r1, r2}]
Plot3D[pdeSoln, {t, 0, 20000}, {rb, r1, r2}, PlotRange -> {0, 6},
AspectRatio -> 10/13, PlotPoints -> 150, MaxRecursion -> 2,
ColorFunction -> "Rainbow"]

回复 xzcyr :两个波包演化到最后都应该达到各自的稳定值，但无论是2000格还是3000格都有一个波包达不到

网格2000和3000的图像，最后都没有达到平稳值

@xzcyr 之所以不令边界处波包等于零，是因为GLb函数一旦变化的话，边界处相对来讲可能并不趋于0。另外此处波包指的是概率密度，直接令概率密度等于0，意味着吸收，我的边界条件是概率流密度为0，即反射。

GLb = 3/(20 rb^7) + 18 rb^7 + (49 PEb rb^9)/4 - (147 rb^8 TEb)/
5; rbs = 0.7790043927824004`;
rbl = 1.8725304891620458`; PEb = 0.4145939392739674`; TEb = 0.65;
pdeSoln =
NDSolveValue[{D[\[Rho][t, rb], t] ==
TEb D[E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb], rb] +
NeumannValue[-TEb D[
E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb], rb],
rb == 0.1] +
NeumannValue[-TEb D[
E^(-GLb/TEb) D[E^(GLb/TEb) \[Rho][t, rb], rb], rb],
rb == 2], \[Rho][0, rb] ==
1/(0.1 Sqrt[\[Pi]]) E^(-((rb - rbs)^2/0.001))}, \[Rho][t, rb], {t,
0, 50000}, {rb, 0.1, 2},
Method -> {MethodOfLines,
SpatialDiscretization -> {FiniteElement,
MeshOptions -> MaxCellMeasure -> (2 - 0.1)/3000}}]
发现自己不会用NeumannValue,求助 @xzcyr