【游戏流体力学基础及Unity代码（十四）】舌尖上的有限元Galerkin法
本文主要介绍有限元方法中的Galerkin方法及如何用其解微分方程。今晚吃啥，一直是一个困扰人类的千古难题。各路仁人志士为此耗费了无数的下午的时光，也没能得出统一的结论。而本文为计算出今晚究竟会吃到啥，提出了一种切实可行的方法，取名为“舌尖上的有限元法”。所有代码在https://github.com/clatterrr/FluidSimulationTutorialsUnity/tree/main/python/FiniteElement
今晚吃啥
假设有一条笔直的一维街道，从x = 0处开始，每隔dx = 1米处就有一个卖包子的晚餐店，每份晚餐包含的包子数量u就是自身所在坐标轴数值。如在x = 0处的就叫“零号店”，每顿晚餐有0个包子，在x = 1处的就叫“一号店“，每顿晚餐1个包子，以此类推。
而你的下午的学习生活工作主要集中在x0处。在度过这个美好下午后，你今晚会去哪家店吃，受到与这家店距离的影响。如果x0刚好在整数处，那么你只会去这里的晚餐店吃。如果不在整数处，那么你的左右各有一家店，你晚餐去这两家店的频率之比就是离这两家店距离的反比。比如你在x0 = 2.7处，那么离二号店距离0.7，离三号店距离0.3。那么你去某一晚去二号店的概率就是0.3，去三号店的概率就是0.7。那么平均每顿晚餐你可以吃上的包子数量U，就是0.3x2+ 0.7x3 = 2.7个包子。
写成公式如下。那个圆中间插一竖那个符号叫phi，代表你去某家店的概率，phi1就是去左边那家店的概率。那个反着的3叫做xi，代表你离你左边那家店的距离。xi的取值范围为，因为距离与你超过1的晚餐店你都不会去。第三个式子，x是生活区的坐标，x1是左边店的坐标。
https://www.zhihu.com/equation?tex=%5Cphi_1+%3D+1+-+%5Cxi+%5Cqquad+%5Cphi_2+%3D+%5Cxi+%5Ctag%7B1%7D+%5Cqquad+%5Cxi+%3D+x+-+x_1
剧透一下，这里的phi就是有限元中插值函数，或者说是Galerkin方法中的权函数，或者基底函数等，或者形函数Shape Function。有的教程会把phi写成N，不过它们是同一个东西。
平均每顿晚餐吃到的包子数U = 去左边包子店的概率 乘以 左边店的每顿包子数 + 去右边包子店的概率 乘以 右边店每顿包子数，这句话翻译成数学公式如下。
https://www.zhihu.com/equation?tex=U+%3D+%5Cbegin%7Bbmatrix%7D+1+-+%5Cxi+%26+%5Cxi+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5Cu_2+%5Cend%7Bbmatrix%7D+%3D+%5Cphi%5ET+u+%5Ctag%7B2%7D
但是，如果你开始不知道每个晚餐店提供多少包子的话，就得自己去计算了。接下来都假设你并不知道每个晚餐的店的包子数量。但是你发现，如果你生活的地点x0是2.7，那么可以吃2.7个包子，而如果是2.8的话，那么你每晚可以吃到2.8个包子。换句话，这意味着生活地点只要前进0.1，那么每晚包子数量也能0.1。再换句话，每晚包子数量增加量 除以 生活地点前进的距离再减去一 等于 零。再再换句话：
https://www.zhihu.com/equation?tex=%5Cfrac%7Bd+U%7D%7Bdx%7D+-+1+%3D+0+%5Ctag%7B3%7D
上式是个伟大的定律，在这定律的背后，必定蕴藏着每个包子店晚餐到底提供多少包子的真相。你的任务就是去探索真相。而浮出幕后的第二个线索是而且x = 1处的1号包子店提供1个包子，这也叫初始条件。
而且你发现，这个定律不受与你去哪家店的概率影响，只要前进0.1，那么包子数量就能增强0.1。我们可以利用这个方法来探索另一个伟大的定律。也就是
https://www.zhihu.com/equation?tex=%5Cint_%7Bx1%7D%5E%7Bx2%7D+%5Cphi%28%5Cfrac%7BU%7D%7Bdx%7D+-+1%29dx+%3D+0+%5Ctag%7B4%7D
x1是左边店的坐标，x2是右边店的坐标。
这个方法也叫做Galerkin法。Galerkin法是有限元中的一种加权残差法。除了Galerkin法还有点装配法和子域法等。好了，你可以把这段话给忘了，继续往下看。
上面这个式子4也可以换成如下的式子。记住dx是生活地点前进的距离，dU是多吃到的包子数量，公式3已经说明它们比值为1。
https://www.zhihu.com/equation?tex=%5Cint_%7Bx1%7D%5E%7Bx2%7D+%5Cphi%28dU%29+%3D+%5Cint_%7Bx1%7D%5E%7Bx2%7D+%5Cphi+dx+%5Ctag%7B5%7D
再来看一下公式2并求导。注意不带上标的u是指晚餐包子铺提供的包子数量。
https://www.zhihu.com/equation?tex=dU++%3D+d%28%5Cphi%5ET+u%29++%3D+d%28%5Cphi%5ET%29u+%3D+%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1x%7D+%26+%5Cphi_%7B2x%7D++%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5Cend%7Bbmatrix%7D%5Ctag%7B6%7D
对phi求导得到如下。
https://www.zhihu.com/equation?tex=%5Cphi_%7B1x%7D+%3D+-1+%5Cqquad+%5Cphi_%7B2x%7D+%3D+1+%5Ctag%7B7%7D
上面这个式子上实际上在说，生活地点每前进0.1，那么去左边店的概率就减少0.1，所以phi{1x} = -1。而去右边店的概率就增加0.1，所以phi{2x} = 1。不带上标的u指每个包子店卖的包子。公式5和6翻译成人话就是
距离左边店的距离变化量 乘 phi1x 乘 左边店的包子数量 + 距离右边店的距离变化量 乘 phi2x 乘 右边店的包子数量 等于包子的总变化量，也就是 0.1x(-1)x2 + 0.1x1x3 = 0.1。
将公式5带入公式4，得到如下式子。注意我已将积分范围设置在左右两家店之间，让左边的店坐标为0，右边店坐标为1。那么生活区的坐标dx其实就是距左边店的距离dxi。这也是使用Galerkin方法解一阶偏微分方程的一般形式。
https://www.zhihu.com/equation?tex=%5Cint_%7B0%7D%5E%7B1%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D++%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1x%7D+%26+%5Cphi_%7B2x%7D++%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2++%5Cend%7Bbmatrix%7Dd+%5Cxi+%3D++%5Cint_%7B0%7D%5E%7B1%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D++%5Cend%7Bbmatrix%7Dd%5Cxi+%5Ctag%7B8%7D
上式的未知数就是左边店与右边店的包子数量u1和u2。对于我们这个问题来说，首先将phi{1x}和phi{2x}代入得
https://www.zhihu.com/equation?tex=%5Cint_%7B0%7D%5E%7B1%7D%5Cbegin%7Bbmatrix%7D+-%5Cphi_%7B1%7D+%26+%5Cphi_%7B1%7D%5C%5C++-%5Cphi_%7B2%7D+%26+%5Cphi_%7B2%7D%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2++%5Cend%7Bbmatrix%7Dd+%5Cxi+%3D++%5Cint_%7B0%7D%5E%7B1%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D++%5Cend%7Bbmatrix%7Dd%5Cxi+%5Ctag%7B9%7D
然后积分
https://www.zhihu.com/equation?tex=%5Cint_0%5E1+%5Cphi_1+d%5Cxi+%3D+%5Cint_0%5E1%281-%5Cxi%29d%5Cxi+%3D+%5Cfrac%7B1%7D%7B2%7D+%5Cqquad+%5Cint_0%5E1+%5Cphi_2+d%5Cxi+%3D+%5Cint_0%5E1%5Cxi+d%5Cxi+%3D+%5Cfrac%7B1%7D%7B2%7D+%5Ctag%7B10%7D
将10代入9得
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+-1%2F2+%26+1%2F2%5C%5C++-1%2F2+%26+1%2F2%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2++%5Cend%7Bbmatrix%7D%3D++%5Cbegin%7Bbmatrix%7D+1%2F2+%5C%5C+1%2F2++%5Cend%7Bbmatrix%7D+%5Ctag%7B11%7D
上面虽然只有两个相同的线性方程组，但别忘了第二个线索即初始条件u1 = 1。由此解得u2 = 2。如果也想顺手知道u3和u4的值，那么上式也可写成下面的样子
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+-1%2F2+%26+1%2F2+%26+0+%26+0%5C%5C++-1%2F2+%26+1%2F2+%26+0+%26+0+%5C%5C0+%26+0+%26+0++%26+0+%5C%5C0+%26+0+%26+0++%26+0%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D%3D++%5Cbegin%7Bbmatrix%7D+1%2F2+%5C%5C+1%2F2+%5C%5C+0+%5C%5C+0+%5Cend%7Bbmatrix%7D+%5Ctag%7B12%7D
上面的式子看起来很沙雕也没什么用，实际上确实很沙雕也没什么用。但是如果将其与下面两个式子相加
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+0+%26+0+%26+0++%26+0+%5C%5C0+%26-1%2F2+%26+1%2F2+%26++0%5C%5C+0+%26+-1%2F2+%26+1%2F2+%26++0+%5C%5C0+%26+0+%26+0++%26+0%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D%3D++%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1%2F2+%5C%5C+1%2F2+%5C%5C+0+%5Cend%7Bbmatrix%7D+%5C%5C+%5Cbegin%7Bbmatrix%7D+0+%26+0+%26+0++%26+0+%5C%5C0+%26+0+%26+0++%26+0+%5C%5C0+%26+0+%26+-1%2F2+%26+1%2F2+%5C%5C+0+%26+0+%26+-1%2F2+%26+1%2F2+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D%3D++%5Cbegin%7Bbmatrix%7D+0+%5C%5C+0+%5C%5C+1%2F2+%5C%5C+1%2F2+%5C%5C+%5Cend%7Bbmatrix%7D
就得到了
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+-1%2F2+%26+1%2F2+%26+0++%26+0+%5C%5C-1%2F2+%26+0++%26+1%2F2+%26++0%5C%5C+0+%26+-1%2F2+%26+0+%26++1%2F2+%5C%5C0+%26+0+%26+-1%2F2++%26+1%2F2%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D%3D++%5Cbegin%7Bbmatrix%7D+1%2F2+%5C%5C+1+%5C%5C+1+%5C%5C+1%2F2+%5Cend%7Bbmatrix%7D+%5Ctag%7B13%7D
这就是平衡矩阵，也可以写成下面的形式
https://www.zhihu.com/equation?tex=Ku+%3D+b+%5Ctag%7B14%7D
在固体力学中，上面的方程也叫刚度方程。K叫做刚度矩阵(stiffness matrix)，刚度矩阵相加也叫装配。如果在处理与弹簧有关的有限元问题中，u就是节点位移，b就是节点力。
13式虽然是4x4的矩阵，但只有三个不同的线性方程组。因此还要一个初始条件。为了求出结果，我们假设u1 = 5，那么上式改写成如下
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+1+%26+0+%26+0++%26+0+%5C%5C-1%2F2+%26+0++%26+1%2F2+%26++0%5C%5C+0+%26+-1%2F2+%26+0+%26++1%2F2+%5C%5C0+%26+0+%26+-1%2F2++%26+1%2F2%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D%3D++%5Cbegin%7Bbmatrix%7D+5+%5C%5C+1+%5C%5C+1+%5C%5C+1%2F2+%5Cend%7Bbmatrix%7D+%5Ctag%7B15%7D
初始条件可以是某个点的值，但也可是某处的梯度。如果u1 = u2，那么矩阵应该如下设置
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+1+%26+-1+%26+0++%26+0+%5C%5C-1%2F2+%26+0++%26+1%2F2+%26++0%5C%5C+0+%26+-1%2F2+%26+0+%26++1%2F2+%5C%5C0+%26+0+%26+-1%2F2++%26+1%2F2%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D%3D++%5Cbegin%7Bbmatrix%7D+0+%5C%5C+1+%5C%5C+1+%5C%5C+1%2F2+%5Cend%7Bbmatrix%7D%5Ctag%7B16%7D
求K的逆矩阵时，暂时用高斯消元法即可。由于主对角线上大部分为零，所以将矩阵k和矩阵b整体上移一位就可用普通高斯消元法。由于上面矩阵K大部分元素都是零，因此又被称为稀疏矩阵。如何利用稀疏矩阵的特性快速求它的逆矩阵，则是另外的话题了。完整代码文件为fem1d.py。其中组装矩阵的部分代码如下。注意我的代码文件中u是从0开始编号的。
for i in range(0,nmax-1):
kmat += [[-0.5,0.5],[-0.5,0.5]]
bmat += 今晚吃啥二维函数季
你好不容易算完了，你的朋友们也都已经出去吃晚餐了，你也正打算出去解决晚餐的时候，你突然想到城市现在一般都是二维的，这样城市被化为一块块的正方形，你应该居住在某个正方形之间。正方形的顶点上各有一个晚餐店，因此你周围有四个晚餐店。晚餐提供包子的数量随x轴和y轴变化。
这样，并且假设店的分布位置如下。那么你离一号店和四号店的x轴距离为x，离1号店和二号店的y轴距离为y。下图括号外的数是晚餐店的编号，而括号内的数是实际的配给的包子数量。
那么你去四个店的概率分别为
https://www.zhihu.com/equation?tex=%5Cphi_1+%3D+%281-x%29%281-y%29+%5Cqquad+%5Cphi_2+%3D+x%281-y%29+%5Cqquad%5Cphi_3+%3D+%281-x%29y++%5Cqquad+%5Cphi_4+%3D+xy+%5Ctag%7B17%7D
上面的概率，或者说是权函数加起来刚好为一。假设每个店提供的包子数量为u，、那么你平均每晚吃到的包子数为
https://www.zhihu.com/equation?tex=U+%3D+%5Cphi_1u_1%2B+%5Cphi_2u_2+%2B+%5Cphi_3u_3+%2B+%5Cphi_4u_4++%3D+%5Cphi%5ET+u%5Ctag%7B18%7D
初始条件可能如下，U沿着x轴或y轴都是一次函数，且ab为常数
https://www.zhihu.com/equation?tex=%5Cqquad+%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+x%7D+%2B+a%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+y%7D+%2B+b++%3D+0+%5Ctag%7B19%7D
然后Galerkin法
https://www.zhihu.com/equation?tex=%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+%5Cphi%28%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+x%7D+%2B+a%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+y%7D+%2B+b%29dxdy+%3D+0+%5Ctag%7B20%7D
https://www.zhihu.com/equation?tex=%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+%5Cphi%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+x%7D+dxdy+%2B+a%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+%5Cphi%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+y%7D+dxdy+%3D+-b%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+%5Cphi+dxdy+%5Ctag%7B21%7D
对于上式左边，继续代入并展开
https://www.zhihu.com/equation?tex=%5Cint_0%5E1%5Cint_%7B0%7D%5E%7B1%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D+%5C%5C%5Cphi_%7B3%7D+%5C%5C+%5Cphi_%7B4%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1x%7D+%26+%5Cphi_%7B2x%7D+%26%5Cphi_%7B3x%7D+%26+%5Cphi_%7B4x%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7Ddxdy+%5C%5C%2Ba%5Cint_0%5E1%5Cint_%7B0%7D%5E%7B1%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D+%5C%5C%5Cphi_%7B3%7D+%5C%5C+%5Cphi_%7B4%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1y%7D+%26+%5Cphi_%7B2y%7D+%26%5Cphi_%7B3y%7D+%26+%5Cphi_%7B4y%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7Ddxdy+%3D+-b%5Cint_0%5E1%5Cint_%7B0%7D%5E%7B1%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D+%5C%5C%5Cphi_%7B3%7D+%5C%5C+%5Cphi_%7B4%7D+%5Cend%7Bbmatrix%7Ddxdy+%5Ctag%7B22%7D
算phi的积分与微分。如果你对二重积分不熟悉，可以看看我的一个回答https://www.zhihu.com/question/44875342/answer/821701089，主题仍然是舌尖上的二重积分。
https://www.zhihu.com/equation?tex=%5Cphi_1+%3D+%281-x%29%281-y%29+%5Cqquad+%5Cphi_2+%3D+x%281-y%29+%5Cqquad%5Cphi_3+%3D+%281-x%29y++%5Cqquad+%5Cphi_4+%3D+xy%5C%5C+%5Cfrac%7B%5Cpartial+%5Cphi_1%7D%7B%5Cpartial+x%7D+%3D+y+-+1+%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_2%7D%7B%5Cpartial+x%7D+%3D+1-y%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_3%7D%7B%5Cpartial+x%7D+%3D+-y%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_4%7D%7B%5Cpartial+x%7D+%3D+y+%5C%5C+%5Cfrac%7B%5Cpartial+%5Cphi_1%7D%7B%5Cpartial+y%7D+%3D+x+-+1+%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_2%7D%7B%5Cpartial+y%7D+%3D+-x%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_3%7D%7B%5Cpartial+y%7D+%3D+1-x%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_4%7D%7B%5Cpartial+y%7D+%3D+x%5C%5C+%5Ciint+%5Cphi_1+%3D+%5Cfrac%7B1%7D%7B4%7D+%5Cqquad+%5Ciint+%5Cphi_2+%3D+%5Cfrac%7B1%7D%7B4%7D+%5Cqquad%5Ciint+%5Cphi_3+%3D+%5Cfrac%7B1%7D%7B4%7D+%5Cqquad%5Ciint+%5Cphi_4+%3D+%5Cfrac%7B1%7D%7B4%7D
接下来就是爆肝的算参数时刻了。
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+-1%2F6+%26++1%2F6+%26+-1%2F12+%26+1%2F12+%5C%5C-1%2F6+%26+1%2F6+%26+-1%2F12+%26+1%2F12+%5C%5C+-1%2F12+%26+1%2F12+%26+-1%2F6+%26+1%2F6+%5C%5C+-1%2F12+%26+1%2F12+%26+-1%2F6+%26+1%2F6+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D+%2B+a%5Cbegin%7Bbmatrix%7D+-1%2F6+%26+-1%2F12+%26+1%2F6+%26+1%2F12+%5C%5C-1%2F12+%26+-1%2F6+%26+1%2F12+%26+1%2F6+%5C%5C+-1%2F6+%26+-1%2F12+%26+1%2F6+%26+1%2F12+%5C%5C+-1%2F12+%26+-1%2F6+%26+1%2F12+%26+1%2F6%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D+%3D-+b%5Cbegin%7Bbmatrix%7D+1%2F4+%5C%5C1%2F4+%5C%5C+1%2F4+%5C%5C+1%2F4%5Cend%7Bbmatrix%7D%5Ctag%7B23%7D
稍微解释一下，左边矩阵最最左上角的元素是这样计算的
https://www.zhihu.com/equation?tex=%5Cint_0%5E1+%5Cint_0%5E1%5Cphi_1+%5Cphi_%7B1x%7D+%3D+%5Cint_0%5E1+%5Cint_0%5E1%281-x%29%281-y%29%28y-1%29+%3D+-%5Cfrac%7B1%7D%7B6%7D+%5Ctag%7B24%7D
同样要解上面的矩阵，必须给出一些点的确切值。完整代码文件为fem2d.py。据我测试，似乎是需要三个初始条件才能正确求解出逆矩阵，但我暂时想不通为什么是三个。
如果把式子19稍微改一下，把对y微分改为对时间t微分，那么就成了线性对流方程
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+t%7D+%2B+a%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+x%7D++%3D+0%5Ctag%7B25%7D
源代码稍经修改，结果如下。从左往右为x轴，从上往下为t轴。U为1的地方，在第零秒还在第一个网格，在第三秒就到了第2个网格。完整代码文件为fem2dadvection.py
今晚吃啥二元函数季
你的朋友们都吃完晚餐回来了，他们告诉你晚餐不仅提供包子，也开始提供馒头了。本着严谨认真的态度，你决定继续研究下去。为了化简这个问题，你继续把城市拍扁成一维的。晚餐店提供包子馒头关系式如下
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+x%7D+%2B+a%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+y%7D+%2B+b%5Cfrac%7B%5Cpartial+V%7D%7B%5Cpartial+x%7D+%3D+0+%5Cqquad+%5Cfrac%7B%5Cpartial+V%7D%7B%5Cpartial+y%7D+%3D+0%5Ctag%7B26%7D
其中U是平均每晚能吃到的包子数量，V是平均每晚吃到的馒头数量。
https://www.zhihu.com/equation?tex=U+%3D+%281-x%29%281-y%29u_1+%2B+x%281-y%29+u_2+%2B+%281-x%29yu_3+%2B+xyu_4+%3D+%5Cphi%5ETu%5C%5C+V+%3D+%281-x%29%281-y%29v_1+%2B+x%281-y%29+v_2+%2B+%281-x%29yv_3+%2B+xyv_4+%3D+%5Cphi%5ETv%5C%5C 虽然包子和馒头食材不同，但这和你晚上去哪个晚餐店没什么关系。
https://www.zhihu.com/equation?tex=%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+%5Cphi%28%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+x%7D+%2B+a%5Cfrac%7B%5Cpartial+U%7D%7B%5Cpartial+y%7D+%2B+b%5Cfrac%7B%5Cpartial+V%7D%7B%5Cpartial+x%7D%29dxdy+%3D+0+%5Ctag%7B27%7D
替换
https://www.zhihu.com/equation?tex=%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+%5Cphi%5Cfrac%7Bd%5Cphi_T%7D%7Bdx%7Dudxdy+%2B+a%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+%5Cphi%5Cfrac%7Bd%5Cphi_T%7D%7Bdy%7Dudxdy+%2B%5Cint_%7B0%7D%5E%7B1%7D%5Cint_%7B0%7D%5E%7B1%7D+b%5Cphi%5Cfrac%7Bd%5Cphi_T%7D%7Bdx%7Dvdxdy+%3D+0+%5Ctag%7B28%7D
与之上一节相比，并不需要额外计算什么东西，最后算得如下式子。
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+-1%2F6+%26++1%2F6+%26+-1%2F12+%26+1%2F12+%5C%5C-1%2F6+%26+1%2F6+%26+-1%2F12+%26+1%2F12+%5C%5C+-1%2F12+%26+1%2F12+%26+-1%2F6+%26+1%2F6+%5C%5C+-1%2F12+%26+1%2F12+%26+-1%2F6+%26+1%2F6+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D+%2B+a%5Cbegin%7Bbmatrix%7D+-1%2F6+%26+-1%2F12+%26+1%2F6+%26+1%2F12+%5C%5C-1%2F12+%26+-1%2F6+%26+1%2F12+%26+1%2F6+%5C%5C+-1%2F6+%26+-1%2F12+%26+1%2F6+%26+1%2F12+%5C%5C+-1%2F12+%26+-1%2F6+%26+1%2F12+%26+1%2F6%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5Cu_3+%5C%5Cu_4+%5Cend%7Bbmatrix%7D+%5C%5C%2B+b%5Cbegin%7Bbmatrix%7D+-1%2F6+%26++1%2F6+%26+-1%2F12+%26+1%2F12+%5C%5C-1%2F6+%26+1%2F6+%26+-1%2F12+%26+1%2F12+%5C%5C+-1%2F12+%26+1%2F12+%26+-1%2F6+%26+1%2F6+%5C%5C+-1%2F12+%26+1%2F12+%26+-1%2F6+%26+1%2F6+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+v_1+%5C%5C+v_2+%5C%5Cv_3+%5C%5Cv_4+%5Cend%7Bbmatrix%7D+%3D+0+%5Ctag%7B29%7D
如果剧透原函数为U = x + y，且V = x。那么我们知道的初始条件式子a = 1，b = -2。并且由于dV/dy = 0，那么v1 = v3而且v2 = v4。再加上初始条件后最后算得，u2 = u3 = v2 = v4 = 1，u1 = v1 = v3 = 0，u4 = 2。完美契合原函数。你实在不信，拿那个4x4矩阵第一行检验一下也行，如下
https://www.zhihu.com/equation?tex=%28-%5Cfrac%7B1%7D%7B6%7Du_1+%2B+%5Cfrac%7B1%7D%7B6%7Du_2+-+%5Cfrac%7B1%7D%7B12%7Du_3+%2B+%5Cfrac%7B1%7D%7B12%7Du_4%29+%2B+a%28-%5Cfrac%7B1%7D%7B6%7Du_1+-+%5Cfrac%7B1%7D%7B12%7Du_2+%2B+%5Cfrac%7B1%7D%7B6%7Du_3+%2B+%5Cfrac%7B1%7D%7B12%7Du_4%29%5C%5C+%2B+b%28-%5Cfrac%7B1%7D%7B6%7Dv_1+%2B+%5Cfrac%7B1%7D%7B6%7Dv_2+-+%5Cfrac%7B1%7D%7B12%7Dv_3+%2B+%5Cfrac%7B1%7D%7B12%7Dv_4%29+%3D+0+
如果你无法在一秒钟之内想出这样的验证方法的话，你应该补一下线性代数与矩阵的知识。这部分编程留给读者了。
今晚吃啥二次函数季
终于算出来了每个晚餐店的馒头包子配给，这使你充满了决心。但你突然发现你其实对馒头没什么兴趣，也不愿意去对面街道吃包子。于是继续把包子店拍扁成一维的，而晚餐点仍然只搭配包子这一种小吃。但是吧，包子种类越来越多，比如喷香鲜肉包，多汁榨菜包等，每个晚餐店提供的包子随坐标轴变化呈二次函数关系。为了找出这个二次函数，你需要三个点来确定这个二次函数。也就是说这次你把每三个店首尾相连的分为一组，你每晚可能去这三店中的一个。也就是如果你住在x = 0到x = 2之间，那么就会受到零一二号店的影响。而如果住在x = 2到x = 4之间，就会受到二三四号店的影响，依次类推。
不过二维函数的权函数有两种，第一种类似于二阶精度前向后向差分，计算点为三点中的最左边或最右边。这种权函数的推导方法为，假设每晚吃到包子数量为
https://www.zhihu.com/equation?tex=U+%3D+c_1x%5E2+%2B+c_2x+%2B+c_3+%5Ctag%7B30%7D
那么
https://www.zhihu.com/equation?tex=U%280%29+%3D+c_3+%5Cqquad+U%281%29+%3D+c_1+%2B+c_2+%2B+c_3+%5Cqquad+U%282%29+%3D+4c_1+%2B+2c_2+%2B+c_3+%5Ctag%7B31%7D
令c = Au，那么
https://www.zhihu.com/equation?tex=A+%3D+%5Cbegin%7Bbmatrix%7D+0+%26+0+%26+1+%5C%5C+1+%26+1+%261+%5C%5C+4+%26+2+%26+1%5Cend%7Bbmatrix%7D%5E%7B-1%7D+%3D+%5Cbegin%7Bbmatrix%7D+1%2F2+%26+-1+%26+1%2F2+%5C%5C+-2%2F3+%26+2+%26-1%2F2+%5C%5C+1+%26+0+%26+0%5Cend%7Bbmatrix%7D+%5Ctag%7B32%7D
那么
https://www.zhihu.com/equation?tex=%5Cphi_1+%3D+%5Cfrac%7B1%7D%7B2%7D%5Cxi%5E2+-+%5Cfrac%7B3%7D%7B2%7D%5Cxi+%2B+1+%5Cqquad+%5Cphi_2+%3D+-+%5Cxi%5E2+%2B+2%5Cxi+%5Cqquad+%5Cphi_3+%3D++%5Cfrac%7B1%7D%7B2%7D%5Cxi%5E2-%5Cfrac%7B1%7D%7B2%7D%5Cxi+%5Ctag%7B33%7D
那么U就应该用如下式子表达
https://www.zhihu.com/equation?tex=U+%3D+%5Cphi_1u_1+%2B+%5Cphi_2u_2+%2B+%5Cphi_3+u_3+%3D+%5Cphi%5ETu+%5Ctag%7B34%7D
验证一下式子的正确性。假设U = x^2 + x + 1，那么U0 = 1，U2 = 3，U3为7。也就是
https://www.zhihu.com/equation?tex=%5Cxi%5E2+%2B+%5Cxi+%2B+1+%3D+%28%5Cfrac%7B1%7D%7B2%7D%5Cxi%5E2+-+%5Cfrac%7B3%7D%7B2%7D%5Cxi+%2B+1%29+%2B+3%28-+%5Cxi%5E2+%2B+2%5Cxi%29+%2B+7%28+%5Cfrac%7B1%7D%7B2%7D%5Cxi%5E2-%5Cfrac%7B1%7D%7B2%7D%5Cxi%29+%5Ctag%7B35%7D
如果选择U1 = 3，U2 = 7，U3 = 13，最后算出来的就是(x+1)^2 + (x+1) + 1，也是正确的。毕竟这次最左边的店需要往左移动一个。但是你也可以将31式子改成算U(0)和U(0.5)和U(1)的，翻译成人话就是现在每隔0.5个网格就有一个包子店，因此每个网格包括首尾有三个包子店，这时候你的活动范围就是这一个网格里的三个包子店。最后算出的权函数也有所不同。这种方法的x的积分域是0到2。不过有的教程喜欢用xi来代替x，xi的积分域是0到1，这样x = 2xi，dx = 2dxi。
另一种权函数类似于中心差分的计算方法，计算的位置在中间。针对一维二次函数的权函数为
https://www.zhihu.com/equation?tex=%5Cphi_1+%3D+%5Cfrac%7B1%7D%7B2%7D%5Cxi%28%5Cxi+-+1%29+%5Cqquad+%5Cphi_2+%3D+1+-%5Cxi%5E2+%5Cqquad+%5Cphi_3+%3D+%5Cfrac%7B1%7D%7B2%7D%5Cxi%28%5Cxi%2B1%29+%5Ctag%7B36%7D
再算如果U = 1，U = 3，U = 7，那么上式计算出(x + 1)^2 + (x + 1) + 1，代入x = 0算出的是第二个值3。这种方法的积分域是-1到1。实际算参数的时候，这两种权函数没什么差别。
不过现在我们暂时用第一种形式的权函数。把假设忘了，我们现在只知道下面的公式，U是晚餐店搭配包子的数量
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial%5E2+U%7D%7B%5Cpartial+x%5E2%7D+-+2+%3D+0+%5Ctag%7B37%7D
且零号店搭配一个包子即U(0) = 1，这时候如何反求出U呢？我们仍然可以给上面积分一个权函数
https://www.zhihu.com/equation?tex=%5Cint_%7B0%7D%5E%7B2%7D+%5Cphi%5Cfrac%7BdU%5E2%7D%7Bdx%7D+%3D+2%5Cint_%7B0%7D%5E%7B2%7D+%5Cphi+dx+%5Ctag%7B38%7D
然后分部积分
https://www.zhihu.com/equation?tex=%5B%5Cphi+%5Cfrac%7BdU%7D%7Bd%5Cxi%7D%5D_0%5E2+-+%5Cint_0%5E2+%5Cfrac%7Bd%5Cphi%7D%7Bd%5Cxi%7D%5Cfrac%7BdU%7D%7Bd%5Cxi%7D+d%5Cxi+%3D+2%5Cint_0%5E2+%5Cphi+d+%5Cxi+%5Ctag%7B39%7D
为了验证上面式子正确性，继续手算。先对权函数求导和积分
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+%5Cphi_1%7D%7B%5Cpartial+%5Cxi%7D+%3D+%5Cxi+-+%5Cfrac%7B3%7D%7B2%7D+%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_2%7D%7B%5Cpartial+%5Cxi%7D+%3D+-2%5Cxi+%2B+2+%5Cqquad+%5Cfrac%7B%5Cpartial+%5Cphi_2%7D%7B%5Cpartial+%5Cxi%7D+%3D+%5Cxi+-+%5Cfrac%7B1%7D%7B2%7D%5C%5C+%5Cint_0%5E2+%5Cphi_1d%5Cxi+%3D+%5Cfrac%7B1%7D%7B3%7D+%5Cqquad+%5Cint_0%5E2+%5Cphi_2d%5Cxi+%3D+%5Cfrac%7B4%7D%7B3%7D+%5Cqquad+%5Cint_0%5E2+%5Cphi_3d%5Cxi+%3D+%5Cfrac%7B1%7D%7B3%7D+
再次假设我们提前拿到剧本，知道U = x^2 + x + 1，那么dU = 2x + 1，然后将phi1代入
https://www.zhihu.com/equation?tex=%28%28%5Cfrac%7B1%7D%7B2%7D%5Cxi%5E2+-+%5Cfrac%7B3%7D%7B2%7D%5Cxi+%2B+1%29%282%5Cxi+%2B1%29%29%7C_0%5E2+-+%5Cint_0%5E2%28%5Cxi+-+%5Cfrac%7B3%7D%7B2%7D%29%282%5Cxi%2B1%29d%5Cxi+%3D+2%5Cint_0%5E2+%28%5Cfrac%7B1%7D%7B2%7D%5Cxi%5E2+-+%5Cfrac%7B3%7D%7B2%7D%5Cxi+%2B+1%29d%5Cxi+%5Ctag%7B40%7D
计算可得
https://www.zhihu.com/equation?tex=%28-1%29+-+%28-%5Cfrac%7B5%7D%7B3%7D%29+%3D+2%28%5Cfrac%7B1%7D%7B3%7D%29+%5Ctag%7B41%7D
确实是正确的。代入phi2或phi3一样也不会出错。推导公式时每一步都要保证正确是很重要的。
继续把假设的解析解忘了，现在我们知道式子分部积分后39式肯定是对的，然后再将34式代入39式，得到
https://www.zhihu.com/equation?tex=%5B%5Cphi+%5Cfrac%7Bd%5Cphi%5ET%7D%7Bd%5Cxi%7Du%5D_0%5E2+-+%5Cint_0%5E2+%5Cfrac%7Bd%5Cphi%7D%7Bd%5Cxi%7D%5Cfrac%7Bd%5Cphi%5ET%7D%7Bd%5Cxi%7Du+d%5Cxi+%3D+2%5Cint_0%5E2+%5Cphi+d+%5Cxi+%5Ctag%7B42%7D
然后展开
https://www.zhihu.com/equation?tex=%28%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D+%5C%5C+%5Cphi_%7B3%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1x%7D+%26+%5Cphi_%7B2x%7D+%26+%5Cphi_%7B3x%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5C+u_3+%5Cend%7Bbmatrix%7D%29_0%5E2-+%5C%5C%5Cint_0%5E2+%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1x%7D+%5C%5C+%5Cphi_%7B2x%7D+%5C%5C+%5Cphi_%7B3x%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1x%7D+%26+%5Cphi_%7B2x%7D+%26+%5Cphi_%7B3x%7D+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5C+u_3+%5Cend%7Bbmatrix%7Ddx++%3D+2%5Cint_0%5E2+%5Cbegin%7Bbmatrix%7D+%5Cphi_%7B1%7D+%5C%5C+%5Cphi_%7B2%7D+%5C%5C+%5Cphi_%7B3%7D+%5Cend%7Bbmatrix%7D+dx+%5Ctag%7B43%7D
计算可得。计算这一步不会出现最终的有限元程序中，这些值在编写有限元程序前就需要手动计算好，然后作为参数填入有限元程序。虽然这是用上述第一种权函数算出来的，但如果用式子35的第二种权函数，也可以算得完全一样的值。注意第二种权函数积分范围为-1到1。
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+3%2F2+%26+-2+%26+1%2F2+%5C%5C+0+%26+0+%26+0%5C%5C+1%2F2+%26+-2+%26+3%2F2++%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5C+u_3+%5Cend%7Bbmatrix%7D-+%5Cbegin%7Bbmatrix%7D+7%2F6+%26+-4%2F3+%26+1%2F6+%5C%5C+-4%2F3+%26+8%2F3+%26+-4%2F3%5C%5C+1%2F6+%26-4%2F3+%26+7%2F6++%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5C+u_3+%5Cend%7Bbmatrix%7D++%3D+2%5Cbegin%7Bbmatrix%7D+1%2F3+%5C%5C+4%2F3+%5C%5C+1%2F3+%5Cend%7Bbmatrix%7D+%5Ctag%7B44%7D
合并一下，得到
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+1%2F3+%26+-2%2F3+%26+1%2F3+%5C%5C+4%2F3+%26+-8%2F3+%26+4%2F3%5C%5C+1%2F3+%26-2%2F3+%26+1%2F3+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5C+u_3+%5Cend%7Bbmatrix%7D++%3D+a%5Cbegin%7Bbmatrix%7D+1%2F3+%5C%5C+4%2F3+%5C%5C+1%2F3+%5Cend%7Bbmatrix%7D+%5Ctag%7B45%7D
很好，饿着肚子算半天算出来个这么个穷酸矩阵，表面上看有三个方程组，其实它们都是一样的。一个方程组没法求出三个变量，因此还需要额外两个初始条件。干脆把它们化成如下形式
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+1%2F3+%26+-2%2F3+%26+1%2F3+%5C%5C+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5C+u_3+%5Cend%7Bbmatrix%7D++%3D+a%5Cbegin%7Bbmatrix%7D+1%2F3%5C%5C+0+%5C%5C+0+%5Cend%7Bbmatrix%7D+%5Ctag%7B46%7D
这个方程可以用于解下面的二阶偏导数，其中a是常量，对于我们的问题而言a = 2。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial%5E2+U%7D%7B%5Cpartial+x%5E2%7D+-+a+%3D+0+%5Ctag%7B47%7D
最终用于解一维二次函数的矩阵如下。在解之前还要把k矩阵最后两行替换上初始条件
https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+1%2F3+%26+-2%2F3+%26+1%2F3+%26+0+%26+0%5C%5C+0+%261%2F3+%26+-2%2F3+%26+1%2F3+%26++0%5C%5C+0%26++0++%261%2F3+%26+-2%2F3+%26+1%2F3+%5C%5C0%26+0%26+0+%26+0+%26+0+%5C%5C+0+%26+0+%26+0+%26+0+%26+0+%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D+u_1+%5C%5C+u_2+%5C%5C+u_3+%5C%5Cu_4+%5C%5C+u_5%5Cend%7Bbmatrix%7D++%3D+a%5Cbegin%7Bbmatrix%7D+1%2F3+%5C%5C+1%2F3+%5C%5C+1%2F3+%5C%5C+0+%5C%5C+0+%5Cend%7Bbmatrix%7D+%5Ctag%7B48%7D
上面就是用有限元法Galerkin法解一维二次函数的一般形式了。完整代码文件为fem1dx2.py。二次函数稍微变化一下就成了流体力学中的扩散方程。
突然你意识到你自己晚餐还没吃上，很想抓起身边某应急食品来尽享美味。但用Galerkin法解NavierStokes方程的准备都做好了，怎么能不立刻开始呢？
这篇所有代码都在https://github.com/clatterrr/FluidSimulationTutorialsUnity/tree/main/python/FiniteElement。下一篇用有限元Galerkin法解NavierStokes方程，敬请期待！可压缩流体那几篇预计再跳票几个月。
上一篇：光影帽子：【游戏流体力学基础及Unity代码（十三）】压力泊松方程,python写SIMPLE算法
参考
《计算流体力学有限元方法及其编程详解》作者：毕超。很棒的解释得很详细的一本书。也是从解微分方程开始解释有限元，而不是从杆弹簧模型之类开始。书中附有matlab代码。
https://zhuanlan.zhihu.com/p/264901969
https://www.bilibili.com/video/BV1d4411i7Wr?from=search&seid=7135824331022518280 有关固体力学中的有限元方法慕课，讲的非常详细很棒
《Finite Element Procedures》
剩下的倒没什么推荐的，因为有限元的书里有50%是使用商业软件的，40%是那种满篇公式初学者很难理解的。剩下的不错的基本上都是在说固体力学。我倒是觉得与其一开始就说杆单元弹簧单元什么的，不如从一阶二阶偏微分方程开始说起，多举一些手算的例子，然后把刚度矩阵形函数这些吓人的名词换成比较亲切的比喻，这样对初学者更容易理解一些。
页:
[1]