上帝说要有光 于是有了麦克斯韦方程。 于是宇宙开始求解PDE。 写这个教程是PDE还是很有上帝的感觉的
（依赖于时间的瞬态问题） 对常微分方程直接进行数值积分的时间逐步积分法 时间步长只受求解精度限制的隐式算法（Newmark法）
求解时间步长受算法稳定限制的显示算法（中心差分法）
方程的建立
微分方程的等效积分形式和加权余量法
$$!A(u)=\left(\begin{array}{c}A_{1}(u)\\ A_{2}(u)\\ A_{3}(u)\\ A_{4}(u)\\ . \end{array}\right )=0$$ $$!B(u)=\left(\begin{array}{c}B_{1}(u)\\ B_{2}(u)\\ B_{3}(u)\\ B_{4}(u)\\ . \end{array}\right )=0$$ 比如对于热传导方程有 $$!A(\phi)=\frac{\partial}{\partial x}(k\frac{\partial\phi}{\partial x})+\frac{\partial}{\partial y}(k\frac{\partial\phi}{\partial y})+Q=0$$ –Equ .0 在Ω上 以及 $$!B(\phi)=\begin{cases}\phi-\overline{\phi}\\k\frac{\partial\phi}{\partial n}-\overline{q}\end{cases}=0$$在Γ上。 –Equ .1 其中$$\phi$$为温度，k为热传导系数$$\overline{\phi},\overline{q}$$分别是边界上温度和热流的给定值 由于方程组在域Ω上每一点都必须为零。所以有 $$!\int_\Omega v^{T}A(u)d\Omega\equiv\int_{\Omega}(v_{1}A_{1}(u)+v_{2}A_{2}(u)+…)d\Omega\equiv 0$$ –Equ .2 其中，$$!v=\left(\begin{array}{c}v_{1}\\ v_{2}\\ v_{3}\\ v_{4}\\ .\\.\\. \end{array}\right )$$ –Equ .3 为函数向量，是一组和微分方程数量相等的任意函数。 其中Εqu.3和Equ.0是完全等效的积分形式。可以断言： 若Equ.3对于任何的v都可以成立，则Equ.0在域内任何一定得到满足。 同样的，B(u)也可以写成如下的积分形式 $$!\int_{\text{Ω}}\overline{v}^{T}B(u)d\Omega\equiv\int_{\Omega}(\overline{v}_{1}B_{1}(u)+\overline{v}_{2}B_{2}(u)+…)d\Omega\equiv 0$$ –Equ .4 也就是说，积分形式 $$!\int_{\Omega}v^{T}A(u)d\Omega+\int_\Omega \overline{v^{T}}B(u)d\Gamma\equiv 0$$
MATHEMATICA
然后就是大部分东西都可以在MMA的Documation直接搜索FEA得到，博主这边主要是根据自己的一些了解进行扩充，顺便加入一些大学二年级三年级物理专业所涉及的东西（包括四大力学的一些经典问题的，以及流体力学相关的内容），来体会下一个数值解PDE到底可以干什么 另外是这里面的代码实例尽可能和MMA保持视觉上的一致，比如$$\in$$在blog里面就是以Latex输入的。所有代码请见我的Github.io
概括
现谈论下PDE 吧，也就是数理方程处理的问题，本来，作为一个数理方程刚及格的学渣，应该面壁，不过鉴于一直对于计算比较感兴趣，在这里妄谈一下数值计算PDE。 如我们所知，数值计算PDE的方法众多，常用于工程领域的包括： 在CFD（计算流体力学）常用的是有限差分法 或者 有限体积法 在结构力学工程领域常见的是有限元方法。 请参照此篇blog 关于FEA的一些讨论，以后会把FEA的代码实现自己过一遍，用GPU加速的方式。这是一种自我磨砺。 http://blog.stlover.org/?p=239 其他对于比较对称的边界条件还有甚多方法，诸如一堆谱方法（基于傅立叶变换）。 先抛去FEA（有限元）之类的格式不谈，我们先说下我们要关注的是什么
PDE
计算一个PDE，要考虑的问题比较复杂，浩浩也是初入门径 首先是计算的物理对象的尺度和要求的精度，对不同的精度和尺度有不同的方法，举个例子，在等离子体的计算中，就有将流体模拟为例子，仅仅将电磁场作为有限元处理的。 再有的问题是我们要处理的PDE的一些特性，比如是椭圆，抛物还是双曲，很重要的一点是是否时间相关。 比如，如果去求解一个电磁场的稳态，或者无粘流体，或者在外力下的平板弯曲，我们可以认为物体会成为一个稳定状态，这时候我们不需要考虑时间相关的因素。 然而在很多问题，诸如，声场计算，基于N－S方程的卡门涡街，这些都是非稳态的问题，我们则必须加入时间这个自变量。 当然，这就有另一个问题了，方程是否是线性或者非线性的。延伸下来就太多了。 其次的问题是，离散化，我们要求一个领域上的解，不得不考虑将这个解变成一堆的采样点，几乎所有的数值PDE计算都是先讲数据离散化，然后再进行处理，这就是涉及到另一个复杂的问题：
网格划分。
比如我们要计算一个机翼的升力，显而易见的是距离机翼比较进的地方需要计算精细一些，距离机翼比较远的地方，需要计算的疏松一些，这中间就需要很多网格划分的问题。 再有的问题是选择合适的“格式”，就是如何表示这些物体，这就是所谓的有限元，有限差分，有限体积的来历了。 不过在本文中，我们暂时不考虑这些问题，仅仅使用Mathematica提供的有限元方法。未来时间充裕浩浩可能会写一个简单的CUDA的N－S方程解算器，再讲不迟。
入手
首先，你需要准备一些东西，主要包括： 1.请学会正确拼写PDE和知道什么是PDE 2.请至少获得数理方程及格的成绩（呃。。浩浩就是不作不死的刚及格） 3.请获得一台电脑 4.安装Mathematica 10，注意Mathematica官网会自动跳转到中文，中文只有9，请切换到英文下载10版 5.将本blog添加到你的收藏夹。 6.所有的源代码都可以在我的Github.io上找到 7.且行且珍惜，卸载MATLAB入我MMA大法好
入门
我们先从最简单的开始，计算一个稳定状况下的热分布 首先，我们知道，对于一个平面温度场而言，其稳定状况下方程可以写做 $$!\Delta u(x,y) = 0$$ 然后我们可以给定第某类边界条件（博主的学渣本质暴露无疑） $$! U(x,y)=\Gamma(x,y)， (x,y)\in \partial \Omega $$ 这样，我们就可以算出来区域$$\Omega$$的温度分布. 怎么做呢，首先是表示在Mathematica里面啦
Equ = {
Laplacian[u[x, y], {x, y}] == 0,
u[x, 0] == 1,
u[x, 1] == 0,
u[1, y] == 0,
u[0, y] == 0
}这就给定了一个完整的方程，注意， [code lang=”mathematica”]Laplacian[u,{x1,x2,x3}][/code]在MMA里表示 $$!\Delta u_{x1,x2,x3}$$ 记得提前要载入FEM包 [code lang=”mathematica”]Needs[“NDSolve`FEM`”][/code] 然后我们将这个方程拿去计算 [code lang=”mathematica”]uif = NDSolveValue[Equ, u, {x, 0, 1}, {y, 0, 1}][/code] 可以得到一个插值后的uif对象 画图 [code lang=”Mathematica”]Plot3D[uif[x, y], {x, 0, 1}, {y, 0, 1}][/code] 后得到 上面是嵌入网页的Mathematica代码，如果你的电脑配置正确应该可以欣赏到一个动态的可以拿鼠标拉来拉去的热流分布图，托管于Github，请点击信任插件，浩浩不保证IE 360浏览器可以顺畅运行。没有的话请自行安装MMA 10以及Chrome浏览器。 如此一来，我们就仰仗最简单的Mathematica计算出来了一个简单的热流分布。 可能你会问我，你在逗我呢，说好的网格划分呢？其实上面的这个方法是MMA在老版本就内置的PDE求解器，仅仅可以求解一些非复杂情况，如果稍微复杂就跑不动了，不过这不妨碍我们拿之开刀。
区域与边界
此部分需要首先对于布尔运算有一定了解。 按照牛老爷子的证明，上帝做的更多的在于第一推动，所谓第一推动嘛，就是边界条件，可见有时候边界比什么都重要。 Mathematica的PDE功能在前几版都是废的，主要原因就在于前几版没有提供漂亮的区域定义与边界定义，而这一版给出了一个优雅的实现。我们先从解析的部分讲起。 首先我们定义一个区域，可以用一种很简单的方法就是表达式，在式为真的时候为领域。MMA 10 就真的给我们提供了这种方法，感谢大民科Wolfram [code lang=”mathematica”]ImplicitRegion[! (x^2 + y^2 < 5), {{x, -20, 20}, {y, -20, 20}}];[/code] 学过c语言的应该能看出这个感叹号的意思 这个式子的意思是，在$$x\in(－20,20)$$的时候，如果式字$$ !(x^2 + y^2 < 5)$$为真，则是我们的领域，也就是一个方块挖出了一个圆 画出图是 [code lang=”mathematica”]EquArea = ImplicitRegion[! (x^2 + y^2 < 5), {{x, -20, 20}, {y, -20, 20}}]; RegionPlot[EquArea, AspectRatio -> Automatic][/code] 如图的样子RegionPlot函数即是画出一个区域 
然后我们需要定义一下边界的情况。MMA 10同样给了一个漂亮的方式给我们（当然不仅仅限于此） [code]DirichletCondition[beqn,pred];[/code] 这个式字的意思呢就是在后面一个式子成立的时候pred成立的时候，边界条件使得beqn成立，就是后面约定了边界的定义，前面约定了边界的条件 比如 [code]DirichletCondition[u[x, y] == a, x^2 + y^2 == 5];[/code] 就表示， x^2 + y^2 == 5 成立的情况下，边界满足u[x, y] == a。 我们现在做一个简单的问题，一个二维的盒子接地里面装了一个二维的等势体，求电场分布。 如果用电动力学的方法，我们可以写出这个空间的勒让德多项式，然后使得满足边界条件，然后得到解。非常繁琐。在这里面我们定义了以下方程 [code lang=”mathematica”]Equ0 = Laplacian[u[x, y], {x, y}] == 0; a = 10; bdn1 = DirichletCondition[u[x, y] == a, x^2 + y^2 == 5]; bdn2 = DirichletCondition[u[x, y] == 0, x == 20 || x == -20 || y == 20 || y == -20]; Equs = {Equ0, bdn1, bdn2};[/code] Equ0表示了静电方程 bdn1是二维等势体（这里是一个圆柱）的电势，我们定义为a bd2则表示了盒子接地，注意我们这里的布尔运算符和类c语言(c,c++,c#,js etc.)类似，但是 a 或 b 只能用 a||b表示，这是因为在类c语言中，很常见的是单竖线表示位运算，而bool是一个bit，其位运算当然和布尔运算相等，但是mma可不认你这一套 这样，我们有了一堆方程，边界还有一个计算区域，我们就可以直接拿去算了 [code]uif = NDSolveValue[Equs, u, {x, y} $$\in$$ EquArea];[/code] 这里{x, y} $$\in$$ EquArea是什么意思各位理工生一看就能理解，这就是在一个区域内，对Equs进行求解。 注意的是符号$$\in$$在MMA需要摁一下ESC在输入el再摁一下ESC输入，这种两下ESC见输入字符是MMA独有的，当然你也可以用Math输入助手，虽然慢了点，但是比你解数理方程一道题一小时能好点。 我们可以注意到一个事情大括号{}在MMA里面很全能，在Equs的定义里面他表达了一个方程数组的意思，在{x,y}则是一个坐标，熟悉的人知道其实这货还是数组，这就是面对符号的MMA的魅力所在，你可以在数组里面存任何东西，你的方程，你的函数，以至于女神同学的照片，或者博主和某同学半年前吵架那天晚上云后的北斗七星的位置（以后再讲MMA的数据库）。 这个解出来的uif实际上是一个差值函数，我们可以很容易画出来他的电势高度和电场向量 [code lang=”mathematica”]p1 = ContourPlot[uif[x, y], {x, y} $$\in$$ EquArea, ColorFunction -> “Temperature”, AspectRatio -> Automatic, PlotLegends -> Automatic]; vetorField = -Grad[uif[x, y], {x, y}] p2 = StreamPlot[ vetorField, {x, -20, 20}, {y, -20, 20}]; Show[p1, p2][/code] 每一个参数是做什么的建议搜索下文档，需要知道MMA是用箭头给可选参数赋值的 这里用了点MMA的矢量运算[code lang=”mathematica”]-Grad[u[x,y],{x,y}][/code] 表示 $$\Delta_{x,y}u(x,y)$$ 然后用stream画出来，注意因为mma的一点bug直接把vectorField那块换成grad是会挂的。 则出现了下图的样子 
然后我们再研究一个问题，作为对于某冬霜的大一的电磁学大作业嘲讽。 某冬霜同学在大一的电磁学免修考试得了98吧分以后，觉得少了两分不得不抱憾终身，于是一边在复习电动力学免修考试，一边迷上了电磁学大作业，冬霜同学研究这么一个问题
一个以接近光速运动的球体电容器是否与一个拉普拉斯变换坐标后形成的椭球体电容器等价
我们先不用关心什么是他喵的光速情况下的球体电容器，我们可以试试算一下椭球体电容器的电容。在大一那个MMA要做矢量运算需要加载包的年头，这货出了从图书馆找到了一本十数年无人问津的电磁书来解析解以外。我们用了一套蒙特卡洛方法来计算这个问题：
给定一堆电子，约束使得他们只能在椭球上运行，每个电子运动起来只受静电力，且会因为摩擦造成动量损失，让这些百万个电子跑一跑，我们就能得到一个椭球上的电荷分布，也就是其电容等属性
大二学过计算物理以后知道这叫蒙特卡洛。 不过我们可以在今天的MMA用一种更加优雅的方式解决 思路如下， 给定义一个很大的三维空间，给一个椭球电容器参照上面，只要给定电势，就可以得到一个电场分布，根据电动力学的边界条件关系，很容易得到椭球电容器每一点的电荷密度，积分得到电容。
网格划分
上面这个问题按照我们之前的思路，似乎很是简单 [code lang=”mathematica”]Equ0 = Laplacian[u[x, y, z], {x, y, z}] == 0; a = 10; bdn1 = DirichletCondition[u[x, y, z] == a, x^2 + 3 y^2 + z^2 == 5]; bdn2 = DirichletCondition[u[x, y, z] == 0, x == 20 || x == -20 || y == 20 || y == -20 || x == 20 || x == -20]; Equs = {Equ0, bdn1, bdn2}; EquArea = ImplicitRegion[! ($$x^2 + 3 y^2 + z^2$$ < 5), {{x, -20, 20}, {y, -20, 20}, {z, -20, 20}}]; RegionPlot3D[EquArea, AspectRatio -> Automatic, PlotStyle -> Directive[Yellow, Opacity[0.5]]] uif = NDSolveValue[Equs, u, {x, y, z} $$\in$$ EquArea];[/code] 但是解这个方程你会发现 这个图不怎么好看嘛 然后
NDSolveValue::bcnop: No places were found on the boundary where x^2+3 y^2+z^2==5 was True, so DirichletCondition[u==10,x^2+3 y^2+z^2==5] will effectively be ignored. >>
MMA似乎傻了，这也就意味着，对于这个复杂的3D问题，MMA丧失了划分网格的能力。 其实之前我们做的虽然没有显式的制定网格，其实MMA都替我们做了这部分的工作，现在我们还是得一步步来吧。 回到那个2D平面的电容问题 我们可以画出uif解的网格 [code lang=”mathematica”]uif[“ElementMesh”][“Wireframe”][/code] 那么如图 
我们所有的计算，都是在这些网格点上完成的。。是不是十分的粗糙，现在我们就来研究如何获得一个合适的，好看的网格点。 首先我们期待，在中心的圆附近的网格可以更加的稠密，在边界则疏松一些。 其次是要够好看 [code lang=”mathematica”]EquArea = ImplicitRegion[ ! ((x)^2 + (y)^2 < 5), {{x, -20, 20}, {y, -20, 20}}]; bmesh = ToBoundaryMesh[EquArea, “BoundaryMeshGenerator” -> {“RegionPlot”, MaxRecursion -> 10}, “MaxBoundaryCellMeasure” -> .5]; mesh = ToElementMesh[bmesh] mesh[“Wireframe”][/code] 这个就是使用一个区域来生成一个网格，我们先生成边界网络bmesh，再使用边界网络生成mesh。 画出来如图 
在外边框有点不好看，不过可以试试嘛 [code lang=”mathematica”]Equ0 = Laplacian[u[x, y], {x, y}] == 0; a = 10; bdn1 = DirichletCondition[u[x, y] == a, x^2 + y^2 == 5]; bdn2 = DirichletCondition[u[x, y] == 0, x == 20 || x == -20 || y == 20 || y == -20]; Equs = {Equ0, bdn1, bdn2}; uif = NDSolveValue[Equs, u, {x, y} \[Element] mesh]; p1 = ContourPlot[uif[x, y], {x, y} \[Element] EquArea, ColorFunction -> “Temperature”, AspectRatio -> Automatic, PlotLegends -> Automatic]; vetorField = -Grad[uif[x, y], {x, y}] p2 = StreamPlot[ vetorField, {x, -20, 20}, {y, -20, 20}]; Show[p1, p2][/code] 结果如图 
好像是要好看一些啦～～ 不过我们对于外边界的这个耿耿于怀，因为如果MaxBoundaryCellMeasure再取得细致，计算量会增大，但是外边界不需要那么细致，特别对于未来的流体是很明显了，但是要如何优化呢 这时候我翻阅了文档，发现我们得采用一种不是那么优美的方式了… 这里我们以后再转门讲 我们可以试试用这个方法处理椭圆
德莱尼三角
预处理工作，也就是Delaunay三角刨分了。
所谓德莱尼三角刨分，是对平面区域根据点集合的一个优雅刨分，使得平面划分更加“科学合理”，此处的“科学合理”我们无妨认为指的是在我们做曲面积分或者线积分的时候有更小误差，因为我们要算东西必须首先分化之。所以德莱尼三角刨分的应用不可谓不少。
首先通过一些简单的例子了解下什么是德莱尼三角刨分，于是有请Mathematica
- (*this is mathematica code*)
- Needs[“ComputationalGeometry`”](*引入计算包*)
- (*产生一个随机点列，以及与之相应的点坐标规则*)
- points = Table[
- {RandomReal[100], RandomReal[100]}
- , {i, 1, 1000}];
- coorrule = Table[
- i -> points[[i]]
- , {i, 1, Length[points]}
- ];
- (*进行三角刨分，以及将三角刨分的结果转化为mathematica图论组建需要的包*)
- dela = DelaunayTriangulation[points];
- fir = Fold[
- (Join[#1,
- Table[#2[[1]] -> #2[[2]][[j]], {j, 1, Length[ #2[[2]] ]}]]) >,
- {}, dela]
- (*画出结果*)
- GraphPlot[fir, VertexCoordinateRules -> coorrule, DirectedEdges -> False, MultiedgeStyle -> None]
结果如下：
哝～一个漂亮的刨分出现了。这就是我们要求的结果。出现尽量更多“科学合理的”刨分。
那我们无妨像这样一个问题，什么样的刨分是“科学合理”的？
1.何谓漂亮？
有这样一种观点，漂亮的东西一般都是更易于存在的，因为人类的直觉和对美的定义一般是区域合理的，所以正常审美中的美女一般有更强的生育能力，再次我们可以认为那些追逐减肥过度的女生是在反人类。同样，刚才我们也看到的是一个漂亮的刨分，他是合理的。因为对同意的面积区域，总有更小的边长将至围起来，这样在微积分数值运算中，我们可以在更小的边长上作近似，也就会得到更合理的结果。加之于地图绘制，合理的刨分能够从高度采样获得漂亮的地图，这一点我们会再后面演示。那么我们可以下一种定义来寻找一个最漂亮的刨分，然后再讨论如何实现它。
显然，我们畏惧的是小角度角的出现（这意味着小面积和长边破坏了积分的误差），那么无妨（事实上是书上告诉我的）做这样的一个定义，最小角最大，次小角，次次小角。。。都是最大一个三角刨分是最合理的。
用数学化的语言表示，就是刨分A的角度排序$$\beta=(a0,a.,…an)$$为A的角度向量，定义$$\beta_0>\beta_1$$为按照字典序依次比较。
这样我们就可以得到，每个角都是最大情况的一个刨分，现在有一个数学家看起来毫不起眼的问题，如何实现？
COMSOL
折腾了俩小时用COMSOL跑出了一点流体。 其实这玩意用起来很无脑的。比我高中时期用的Fluent用起来简单多了，对付一些不复杂的需求还是可以的，比如经典问题，卡门涡街。 
SOLIDWORKS
四个螺丝给到100N
在接口前的地方是最脆弱的地方，与之前的认识完全相反，一想也是的，位移在2mm碳纤维的情况下有3mm，在2mm铝合金是8mm！顶点在2mm钛合金是4-6mm 于是做了个加强件，优化极其明显，最大应力下降三个数量级，调到了大气压水准，位移掉了六个数量级。。。
拜读楼主文章！