关于NDSolve在运行时所使用的方法

Question

想请教两个问题：

1.  Mathematica 9（以及这之前的版本）中NDSolve解微分方程时，它用到的方法是否全都是有限差分法？

2.  Mathematica 10 中用NDSolve解微分方程时，在用户不主动加载有限元法的程序包（并且也没有在任何选项中指定）时，是否在某些时候也会自动用到有限元法？

 

补充：问题1中用词有纰漏。这里我只是说到“有限差分法”，其实本意也包括和差分法原理相同或相似的方法，诸如ODE中那些单步法、多部法等（也不知它们是否能称为差分法）。都是沿着（某个）自变量的方向一步步向后递推计算的。

Comments

看这里：http://mathematica.stackexchange.com/q/145/1871 ，顺便虽然这帖的回答获得的投票都很多并且也产生了被采纳的答案，但是实际上问题并没有这么简单；再顺便，你的原始提问倒还没啥，可补充部分却句句是糟点，现在没时间写长答案，只简单提一下：1. 仔细读读“有限差分法”的维基百科；2. 查查“隐式差分格式”。

---

其实我认为隐式差分格式和我补充里提到的那些也是类似的原理，只是有解方程的需求。也许原理一词不准确？

比如微分几何里有局域（local）和全局（global）的区分；类似的我觉得有限差分法就对应这种“局域”的概念，而有限元法就是对应global的概念的。

不过我这些只算是个人感觉，不曾仔细研究。

---

我那里的表述也确实很不严密。

刚看了那个链接，受用。

---

所以说哪里“类似”了啊？类似在它们把连续的方程给离散化了？要是这样的话我好像还真不记得有哪个微分方程数值解法在根本上用的不是“将连续方程用离散近似”的思想，那全体数值解法都可以改名叫差分法了；微分几何我不懂，但是把有限差分对应到局域、有限元对应到全局也是显然不对的，因为至少就我所见的例子来看，有限元法的网格若是规则的，那它就可以等效于隐式差分格式。有限元法和有限差分的区别，就我所见，应该在于网格是否规则。（尽管它们在构建基本理论时的出发点不太相同。）关于这部分内容可以参看Olver所写的 Introduction to Partial Differential Equations的相关章节。（这本书原本在作者主页上有预印本下载，不过现在本书已经出版，所以下载链在作者主页上已经没了，但是，你懂的。）

---

哦你是这样看的，你说的有道理。那本书谢谢推荐！

说到“全局”一词，一个闭合曲面的面积就是这个曲面的一个全局量。有限元法的基本思想是用变分法将解微分方程的问题转换为一族允许曲面的能量积分取极值的问题，这个能量积分也正是一个全局概念。我就是这个思路。

Answer 1

ndsolveMethod[expr___] := 
 DeleteDuplicates /@ {DeleteCases[
    Cases[Select[
      Flatten@Trace[NDSolve[expr], TraceInternal -> True], ! 
        FreeQ[#, NDSolve`MethodData] &], 
     x_Symbol :> 
      If[Context[x] == "NDSolve`" && StringFreeQ[SymbolName[x], "$"], 
       SymbolName@x, Nothing], Infinity, Heads -> True], 
    "MethodData"], 
   Cases[Flatten@Trace[NDSolve[expr], TraceInternal -> True], 
    HoldForm[Method -> ___], Infinity]}

可以用这个查看解特定方程的时候使用的方法

更新：mma9运行题主给的code：

Comments

-----------------很棒！---------------

谢谢！

Answer 2

如图如图如图如图如图如图如图如图如图

Comments

嗯，感觉明晰不少。谢谢！

另外不知这个出处是哪儿？我还想对照一下版本。

---

http://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.html#3826

---

哦，好的。我搜了搜发现里面没有提有限元法，我想这个应该是对应版本9的。

Answer 3

再次赞一下happyfish给出的ndsolveMethod代码。我在没有加载有限元程序包的情况下，用ndsolveMethod对MMA文档中的一个微分方程测试过，其结果显示，的确使用了有限元法（我用的是MMA10.3）。也就是说，在版本10中，有些情况下NDSolve会自动调用有限元法来解方程。

下面是测试用的完整代码：

uif=NDSolve[{Laplacian[u[x,y],{x,y}]==0,
u[x,0]==0,u[x,1]==1,
u[0,y]==0,u[1,y]==0},
u,{x,0,1},{y,0,1}]

Plot3D[u[x,y]/.uif[[1]],{x,0,1},{y,0,1},
ColorFunction->"TemperatureMap"]

ndsolveMethod[expr___]:=DeleteDuplicates/@{DeleteCases[Cases[
Select[Flatten@Trace[NDSolve[expr],TraceInternal->True],
!FreeQ[#,NDSolve`MethodData]&],x_Symbol:>If[
Context[x]=="NDSolve`"&&StringFreeQ[SymbolName[x],"$"],
SymbolName@x,Nothing],Infinity,Heads->True],
"MethodData"],Cases[Flatten@Trace[NDSolve[expr],
TraceInternal->True],HoldForm[Method->___],Infinity]};

ndsolveMethod[{Laplacian[u[x,y],{x,y}]==0,
u[x,0]==0,u[x,1]==1,
u[0,y]==0,u[1,y]==0},
u,{x,0,1},{y,0,1}]

以上代码在MMA10.3中运行后相关结果为：

{Method->{FiniteElement,IntegrationOrder->Automatic,
InterpolationOrder->Automatic,PrecomputeGeometryData->True}}

但不知在本版9中会怎样。不知谁电脑上还装有Mathematica 9？可否帮忙运行一下以上程序，把结果贴出来，看一看都用到什么方法？

 

（顺便提一下，我的原问题后来有所补充修正。）

Comments

我运行了一下，出了一堆错误，貌似9版本解不了那个方程。截图更新在我的答案里了

---

原来如此啊！谢谢啦。

那个方程的边界条件看起来也挺规整的，竟然非有限元不能解。。。也不知为什么不能解。