关于（遵循Fokker-Planck的）波包演化的问题

1. 你为什么觉得这个结果是有问题的？你所期待的结果是啥？
2. 你确定方程本身正确吗？可否给一下出处？
3. 默认情况下 eerr 警告的误差显然大了些，而尝试加密空间网格后时间步进的速度开始变得很小。（在我写这答案这段时间里还没算完。）你确定时间区间的设置合理吗？
4. 你是不是有一层楼被吞/仅自己可见了？是的话最好申诉下。

5. 你的两个边界是否是对无穷远的近似？

……折腾了半天发现还是网格不够密。添加选项
Method -> mol[2000, 4]
即可，其中mol的定义是
mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
2 GHz的电脑，耗时约 60 秒；边界其实可以用两个为 0 的第一类边界代替，求解速度略快，50 秒。（不过我只测了一次。）

……坑了个爹的，NDSolve好像又出现了性能退化。同样的代码（mol[3000, 4]），版本9画出来是这样的：
Plot[pdeSoln /. {{rb -> rbs}, {rb -> rbl}} // Evaluate, {t, 0, tend}]
pdeSoln[[0]]["Coordinates"][[1]] // ListPlot

版本12.2是这样的：

12.2明显不对劲。

嗬，这个问题居然会用到这个我是真没想到。
仔细试验后发现，本问题若想高效求解，必须使用特殊的方式进行空间离散，此问题在Stackexchange帖子《How to solve the tsunami model and animate the shallow water wave?》（编号107988）里有所论及，简单地说就是，对于某些基于守恒率的偏微分方程，如果空间离散不将守恒率以某种方式纳入考量，数值求解就会出问题。（求解效率低下，甚至无法得出正确解。）以往见到的例子，守恒率内部的项貌似都具有某种非线性，可LZ的这个方程只是变系数，居然也有类似问题，还真是让我开了眼了。
总之上代码。离散将使用 pdetoode 函数（mathematica.stackexchange.com/a/127997/1871）进行辅助。没改的参数我就不重复贴了：
{eq, ic, bc} = {D[ρ[t, rb], t] == 500 TEb D[mid[t, rb], rb], ρ[0, rb] ==
1/(0.1 Sqrt[π]) E^(-((rb - rbs)^2/0.01)), ρ[t, rb] ==
0 /. {{rb -> r1}, {rb -> r2}}};
eqmid = mid[t, rb] == E^(-GLb/TEb) D[E^(GLb/TEb) ρ[t, rb], rb];
(*eqmid2=mid2[t,rb]==E^(GLb/TEb) ρ[t,rb];*)
tend = 2 10^4;
domain = {r1, r2};
points = 3000; difforder = 2;
grid = Array[# &, points, domain];
ptoofunc = pdetoode[{ρ, mid}[t, rb], t, grid, difforder];
del = #[[2 ;; -2]] &;
ode = Block[{mid}, Evaluate@ptoofunc[eqmid[[1]]] = ptoofunc[eqmid[[2]]];
ptoofunc[eq] // del];
odeic = ptoofunc[ic];
odebc = With[{sf = 1}, Map[sf # + D[#, t] &, bc, {2}] // ptoofunc];
sollst = NDSolveValue[{ode, odeic, odebc}, ρ /@ grid, {t, 0, tend}, jianshi[t](*,
SolveDelayed\[Rule]True*)]; // AbsoluteTiming
sol = rebuild[sollst, grid]
计算只需约15 s。结果：
Plot[sol[t, rb] /. {{rb -> rbs}, {rb -> rbl}} // Evaluate, {t, 0, tend}]
Plot3D[sol[t, rb], {t, 0, tend}, {rb, r1, r2}, PlotRange -> All,
PlotPoints -> 50]

可算像话了嗯。

……我傻了。仔细想想，有限元法的标准形式就是基于守恒率的啊：
bc = \[Rho][t, rb] == 0 /. {{rb -> r1}, {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 - rbs)^2/0.01)), bc}, \[Rho][t, rb], {t, 0,
tend}, {rb, r1, r2},
Method -> {MethodOfLines,
SpatialDiscretization -> {FiniteElement,
MeshOptions -> MaxCellMeasure -> (r2 - r1)/3000}}]; // AbsoluteTiming
耗时只要1秒不到。结果和13楼是一样的。