嗯……吧里关心微分方程的人不多，不过这个月新帖也少，姑且发个帖说说近期出现在Stackexchange的这个我觉得还比较有意思的帖子。我们知道，对于不规则区域上的偏微分方程求解，目前有限元法（finite element method, FEM）几乎可以说是一家独大（我甚至见过某些对微分方程求解不甚了解的童鞋竟模模糊糊地觉得“解偏微分方程=用一下有限元”）；目前，在Mathematica的NDSolve函数的内置方法里，它也是唯一一种可以在不规则域上求解偏微分方程的方法。但是，可以用于不规则域的偏微分方程数值求解方法其实远不止有限元一种，今天我要讲的就是一种较为小众（？）的无网格偏微分方程求解方法——广义有限差分法（generalized finite difference method, GFDM），也称无网格有限差分法（meshless finite difference method, MFDM）。
我们知道传统的有限差分法（finite difference method, FDM）所用的差分公式是基于单变量泰勒展开的，这使得其通常只适用于规则区域的求解（当然了，基于外差法或是坐标变换将FDM用于不规则域是一个可行的思路，但是非常繁琐，一个例子是SE帖子《Poisson PDE over a irregular region with FDM》（编号：64585）），GFDM的改进，在于使用多元泰勒展开（不是内置的那个，而是“点附近的展开”，可以参看《Multivariable Taylor expansion does not work as expected》（编号：15023））代替单变量泰勒展开，也就是说，对于选定的中心结点，我们选定它附近的n个点，分别作n次多元m阶泰勒展开（m和n不必相同），获得n个方程，再使用最小二乘法求解这n个方程，即可得到该点处的1至m阶导数——显然，GFDM的这套思想还是非常简单自然有吸引力的。
好了，重点来了。近期，在SE帖子《How can I construct the derivative matrix for an irregular domain?》（编号：131705）下面，Ulrich Neumann 终于对GFDM作了尝试，这一尝试最终被整理成了一个基于GFDM法的微分矩阵计算程序包，代码如下：
BeginPackage["GFDM`"];
taylorterms; matGFDM;
Begin["`private`"];
taylorterms[dim_, order_] :=
Block[{t, f, vars, difforder, x, xlst}, xlst = x /@ Range@dim;
With[{expr =
Normal@Series[f @@ (t xlst), {t, 0, order}] /. t -> 1}, {vars,
difforder} =
Cases[expr, a : Derivative[o__][f][__] :> {a, {o}}, Infinity]\[Transpose];
{difforder,
Compile @@ {{#, _Real, 1} & /@
xlst, (CoefficientArrays[expr, vars] // Last // Normal)}}]]
w = Function[{d, dm}, #] &[
With[{x = d/dm}, If[d <= dm, 1 - 6 x^2 + 8 x^3 - 3 x^4, 0]] //
PiecewiseExpand // Simplify`PWToUnitStep];
dmfunc[sf_, pcenter_, pinf_] := sf Sqrt[# . # &@(pcenter - pinf)];
matGFDM[dim_, order_, plst_, neighborsize_,
sf_ : False] := (meshsize = Length[plst];
nf = Nearest[plst -> "Index"];
neighbors = nf[plst, neighborsize + 1];
{difforderlst, taylortermfunc} = taylorterms[dim, order];
With[{lst = difforderlst}, take = First@FirstPosition[lst, #] &];
diffmat = Function[indexlst, pcenter = plst[[indexlst[[1]]]];
If[NumericQ@sf, dm = dmfunc[sf, pcenter, plst[[indexlst[[-1]]]]];
wlst =
w[Sqrt[((pcenter - Transpose@plst[[indexlst // Rest]])^2 //
Total)], dm], wlst = 1.];
sparse =
Module[{spa =
SparseArray[{{Range@neighborsize,
indexlst // Rest}\[Transpose] ->
ConstantArray[1., neighborsize]}, {neighborsize, meshsize}]},
spa[[All, First@indexlst]] = -1.; wlst spa];
pinv =
PseudoInverse[
wlst Transpose[
taylortermfunc @@ (plst[[indexlst // Rest]]\[Transpose] -
pcenter)]];
pinv . sparse] /@ neighbors;
With[{mat = diffmat\[Transpose], take = take}, mat[[take@#]] &]);
End[];
EndPackage[];
程序包的使用方式如下：
(* 在给定的方形区域 mesh 上计算微分矩阵 *)
<< NDSolve`FEM`
mesh = ToElementMesh[Rectangle[],
"MeshElementType" -> "TriangleElement", "MeshOrder" -> 1];
plst = mesh["Coordinates"];
dim = 2; order = 8; size = 29;
getmat = matGFDM[dim, order, plst, size]; // AbsoluteTiming
(*{0.194138,Null}*)
(* 提取混合导数 D[..., x, y] 的微分矩阵 *)
dxymat = getmat@{1, 1};
(* 验算 *)
Clear[x,y]
f = Function[{x, y}, Sin[2 Pi (x - Pi/2)] Cos[Pi y]];
funcvalues = f @@ (plst\[Transpose]);
dfunc = Function[{x, y}, D[#, x, y] &@f[x, y] // Evaluate];
ref = dfunc @@@ plst;
dxylst = dxymat . funcvalues;
dgfdmfunc = ElementMeshInterpolation[{mesh}, dxylst];
Manipulate[
Plot[{dgfdmfunc[x, y], dfunc[x, y]}, {x, 0, 1}, PlotRange -> 20,
PlotStyle -> {Automatic, Dashed}], {y, 0, 1}]
如大家所见，这个精度并不是特别理想，尽管数值基本能对上，但是误差肉眼可见（要知道这个计算足足用了29个点作最小二乘），并且这个精度似乎也远远不及GFDM相关论文中所声称的精度，至于是不是我对论文的理解有问题我就不太清楚了。
部分论文在计算时使用了权函数，但就我个人的测试结果而言，效果似乎不明显。
大家觉得是哪里出了问题？
另，在上述 131705 号帖子中还有径向基函数-有限差分法（radial basis function-finite difference, RBF-FD，也是一种以FDM后继者自居的方法）的编程尝试，效果也不是很好，有兴趣的童鞋可以看看。
我们知道传统的有限差分法（finite difference method, FDM）所用的差分公式是基于单变量泰勒展开的，这使得其通常只适用于规则区域的求解（当然了，基于外差法或是坐标变换将FDM用于不规则域是一个可行的思路，但是非常繁琐，一个例子是SE帖子《Poisson PDE over a irregular region with FDM》（编号：64585）），GFDM的改进，在于使用多元泰勒展开（不是内置的那个，而是“点附近的展开”，可以参看《Multivariable Taylor expansion does not work as expected》（编号：15023））代替单变量泰勒展开，也就是说，对于选定的中心结点，我们选定它附近的n个点，分别作n次多元m阶泰勒展开（m和n不必相同），获得n个方程，再使用最小二乘法求解这n个方程，即可得到该点处的1至m阶导数——显然，GFDM的这套思想还是非常简单自然有吸引力的。
好了，重点来了。近期，在SE帖子《How can I construct the derivative matrix for an irregular domain?》（编号：131705）下面，Ulrich Neumann 终于对GFDM作了尝试，这一尝试最终被整理成了一个基于GFDM法的微分矩阵计算程序包，代码如下：
BeginPackage["GFDM`"];
taylorterms; matGFDM;
Begin["`private`"];
taylorterms[dim_, order_] :=
Block[{t, f, vars, difforder, x, xlst}, xlst = x /@ Range@dim;
With[{expr =
Normal@Series[f @@ (t xlst), {t, 0, order}] /. t -> 1}, {vars,
difforder} =
Cases[expr, a : Derivative[o__][f][__] :> {a, {o}}, Infinity]\[Transpose];
{difforder,
Compile @@ {{#, _Real, 1} & /@
xlst, (CoefficientArrays[expr, vars] // Last // Normal)}}]]
w = Function[{d, dm}, #] &[
With[{x = d/dm}, If[d <= dm, 1 - 6 x^2 + 8 x^3 - 3 x^4, 0]] //
PiecewiseExpand // Simplify`PWToUnitStep];
dmfunc[sf_, pcenter_, pinf_] := sf Sqrt[# . # &@(pcenter - pinf)];
matGFDM[dim_, order_, plst_, neighborsize_,
sf_ : False] := (meshsize = Length[plst];
nf = Nearest[plst -> "Index"];
neighbors = nf[plst, neighborsize + 1];
{difforderlst, taylortermfunc} = taylorterms[dim, order];
With[{lst = difforderlst}, take = First@FirstPosition[lst, #] &];
diffmat = Function[indexlst, pcenter = plst[[indexlst[[1]]]];
If[NumericQ@sf, dm = dmfunc[sf, pcenter, plst[[indexlst[[-1]]]]];
wlst =
w[Sqrt[((pcenter - Transpose@plst[[indexlst // Rest]])^2 //
Total)], dm], wlst = 1.];
sparse =
Module[{spa =
SparseArray[{{Range@neighborsize,
indexlst // Rest}\[Transpose] ->
ConstantArray[1., neighborsize]}, {neighborsize, meshsize}]},
spa[[All, First@indexlst]] = -1.; wlst spa];
pinv =
PseudoInverse[
wlst Transpose[
taylortermfunc @@ (plst[[indexlst // Rest]]\[Transpose] -
pcenter)]];
pinv . sparse] /@ neighbors;
With[{mat = diffmat\[Transpose], take = take}, mat[[take@#]] &]);
End[];
EndPackage[];
程序包的使用方式如下：
(* 在给定的方形区域 mesh 上计算微分矩阵 *)
<< NDSolve`FEM`
mesh = ToElementMesh[Rectangle[],
"MeshElementType" -> "TriangleElement", "MeshOrder" -> 1];
plst = mesh["Coordinates"];
dim = 2; order = 8; size = 29;
getmat = matGFDM[dim, order, plst, size]; // AbsoluteTiming
(*{0.194138,Null}*)
(* 提取混合导数 D[..., x, y] 的微分矩阵 *)
dxymat = getmat@{1, 1};
(* 验算 *)
Clear[x,y]
f = Function[{x, y}, Sin[2 Pi (x - Pi/2)] Cos[Pi y]];
funcvalues = f @@ (plst\[Transpose]);
dfunc = Function[{x, y}, D[#, x, y] &@f[x, y] // Evaluate];
ref = dfunc @@@ plst;
dxylst = dxymat . funcvalues;
dgfdmfunc = ElementMeshInterpolation[{mesh}, dxylst];
Manipulate[
Plot[{dgfdmfunc[x, y], dfunc[x, y]}, {x, 0, 1}, PlotRange -> 20,
PlotStyle -> {Automatic, Dashed}], {y, 0, 1}]
如大家所见，这个精度并不是特别理想，尽管数值基本能对上，但是误差肉眼可见（要知道这个计算足足用了29个点作最小二乘），并且这个精度似乎也远远不及GFDM相关论文中所声称的精度，至于是不是我对论文的理解有问题我就不太清楚了。
部分论文在计算时使用了权函数，但就我个人的测试结果而言，效果似乎不明显。
大家觉得是哪里出了问题？
另，在上述 131705 号帖子中还有径向基函数-有限差分法（radial basis function-finite difference, RBF-FD，也是一种以FDM后继者自居的方法）的编程尝试，效果也不是很好，有兴趣的童鞋可以看看。
