【游戏流体力学基础及Unity代码（十六）】非线性有限元及牛顿迭代法

Ylisar · 发表于 2021-5-9 21:35

牛顿迭代法

大自然中很多现象都是非线性的。因此来看看怎么用有限元Galerkin方法解决非线性问题吧。首先，所有有限元方法都是把待求解的问题转换成矩阵形式如下。

$F u =b \tag{1}$
u是未知变量。F是已知刚度矩阵。b也已知。如果问题原来是线性的，那么刚度矩阵里面就是常数，求个逆能稳定且准确地算出答案。但是如果F很稀疏，可以用线性方程组的迭代算法来加快速度，比如Jacobi或GuassSeidel或SOR。这三个算法我们挺熟悉的，第十三章解压力泊松方程，与上面一章算逆矩阵也用到了它们。
传统的GaussSeidel和SOR方法虽然速度更快，但由于每次迭代时都要更新很多次结果，所以并不适合并行。因此又衍生出一大堆适合并行的算法。其具体算法可以在各种线性方程数值求解书上看到。此外还有共轭梯度法也可用于解线性方程组，可参考《An Introduction to the Conjugate Gradient Method Without the Agonizing Pain》
比如上一章我的代码就用到了雅可比迭代。
void Jacobi(){
int tmax = 100;
Vector3[] res2 = new Vector3[total];
while (tmax-- > 0)
{
      for(int i = 0; i < total;i++)res2 = res;
      for (int i = 0; i < total; i++)
      {
         Vector3 summ = new Vector3(0.0f,0.0f,0.0f);
         for (int j = 0; j < total; j++)
         {
            if (i == j) continue;
            summ += Kmat[i, j] * res2[j];
         }
         res = (bMat-summ)/ Kmat[i,i];
      }
}
然而如果问题是非线性的，那么就要用非线性方程组(nonlinear system of equation)的迭代算法来算了。常用的算法是牛顿-拉弗森法(Newton-Rapson Method)，也叫牛顿法。公式如下

$x^{k+1} = x^k - \frac{F(x)}{F'(x)} \tag{2}$
这是从一阶泰勒公式推出来的。不过对面最优化数值算法中也有个牛顿迭代法，与上面的牛顿法有一点区别

$x^{k+1} = x^{k} - \frac{F'(x)}{F''(x)} \tag{3}$
这玩意是求最小值的。稍加分析就会发现，非线性方程求根中的牛顿法，只有当F(x) = 0的时候才会停止迭代，此时x就根。而求最小值的牛顿法，只有当F&#39;(x) = 0的时候才会停止迭代，此时斜率为零，x处于极值处。
这种数值最优化方法，隔壁机器学习领域最喜欢用了，因此又发展出来一堆梯度下降，随机梯度下降，动量梯度下降等算法。能见到这些最优化方法的还包括控制理论和凸分析中。数值最优化虽然和非线性方程求根不同，但都是与梯度有关迭代算法，因此各种收敛理论与加速求解技巧也可相互借鉴。
比如求x^4 - x^2 - 2x，它的图形如下，最小值F(-1) = -2，根为0和1.52
写成python代码如下
import numpy as np
tmax = 100
root = np.zeros((tmax))# 求根
minval = np.zeros((tmax))# 求最小值
root[0] = minval[0] = 5 # 初始猜测值
def F(x):# 最小值F(-1) = -2，根为0和1.52
return x*x*x*x - x*x -2*x
for i in range(1,tmax):
f = F(root[i-1])
fx = (-F(root[i-1] + 2) + 8*F(root[i-1]+1) - 8*F(root[i-1]-1) + F(root[i-1] - 2))/12
root = root[i-1] - f/fx

fx = (-F(minval[i-1] + 2) + 8*F(minval[i-1]+1) - 8*F(minval[i-1]-1) + F(minval[i-1] - 2))/12
fxx = (-F(minval[i-1] + 2) + 16*F(minval[i-1]+1) -30*F(minval[i-1]) +  16*F(minval[i-1]-1) + F(minval[i-1] - 2))/12
minval = minval[i-1] - fx/fxx经过几次迭代都能很快收敛。不过上面是仅有一个自变量的方程而已，如果有很多个自变量的话，就得稍加修改。下面两张图来自《数值方法(matlab版)》by John H.Mathews第3.7节。
不过在有限元方法中，我们经常要面对的是这种形式的问题

$\begin{bmatrix} 1 \\ -2\end{bmatrix}\begin{bmatrix} 1 & 2\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix}\begin{bmatrix} 5 & 6\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix} + \begin{bmatrix} x \\ y\end{bmatrix} - \begin{bmatrix} 432 \\ -854\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix} \tag{4}$
上面的矩阵可以有这种变化

$\begin{bmatrix} 1 \\ 2\end{bmatrix}\begin{bmatrix} 1 & 2\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix}\begin{bmatrix} 5 & 6\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix} = \begin{bmatrix} 1 \\ 2\end{bmatrix}\begin{bmatrix} 1 & 2\end{bmatrix}\begin{bmatrix} 5 & 6\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix} \\= \begin{bmatrix} 1 \\ 2\end{bmatrix}\begin{bmatrix} 5 & 6 & 10 & 12\end{bmatrix}\begin{bmatrix} x^2 \\ xy \\ yx \\y^2\end{bmatrix} \tag{6}$
上面这些都可以写成下面这种形式

$f_1(x，y) = (x + 2y)(5x + 6y) + x - 432 = 0\\ f_2(x,y) = -2(x + 2y)(5x + 6y) + y + 854 = 0\tag{7}$
求解的话，第一步，先求雅可比矩阵

$J= \begin{bmatrix} df_1/dx & df_1/dy \\ df_2/dx & df_2/dy\end{bmatrix} = \begin{bmatrix} 10x + 16y + 1 & 16x + 24y \\ -20x - 32y & -32x - 48y + 1\end{bmatrix} \tag{8}$
第二步，初始值瞎猜一个x = 0，y = 0。代入原方程和雅可比矩阵

$F(0,0) = \begin{bmatrix} -432 \\ 854\end{bmatrix} \qquad J(0,0) =\begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix}\tag{9}$
第三步，然后算dx和dy，注意等号右边有个负号

$\begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y\end{bmatrix} = -\begin{bmatrix} -432 \\ 854\end{bmatrix} \tag{10}$
第四步，算出dx = 432，dy = 862。然后算出得到xy的下一步迭代值

$\begin{bmatrix} x \\ y\end{bmatrix} = \begin{bmatrix} 0 \\ 0\end{bmatrix} + \begin{bmatrix} \Delta x \\ \Delta y\end{bmatrix} = \begin{bmatrix} 432 \\ -854\end{bmatrix} \tag{11}$
解得x = 432，y = -854。完成这次迭代。然后继续回到第二步，把x = 432和y = -854作为初始值，继续下一轮迭代。继续这么下去，对于这个方程来说第十二次迭代就能得到正确的结果x = 3以及y = 4。写成python代码如下
import numpy as np
def fun(n,x,y):
if n == 0:
      return (x + 2*y)*(5*x + 6*y) + x - 432
else:
      return -2*(x+2*y)*(5*x+6*y)+y + 854
tmax = 100
val = np.zeros((2,tmax))
val[0,0] = 0
val[1,0] = 0
for t in range(0,tmax-1):
#二阶精度的中心差分求雅可比矩阵，对于本方程来说够用了
a = (fun(0,val[0,t]+1,val[1,t]) - fun(0,val[0,t]-1,val[1,t]))/2
b = (fun(0,val[0,t],val[1,t]+1) - fun(0,val[0,t],val[1,t]-1))/2
c = (fun(1,val[0,t]+1,val[1,t]) - fun(1,val[0,t]-1,val[1,t]))/2
d = (fun(1,val[0,t],val[1,t]+1) - fun(1,val[0,t],val[1,t]-1))/2
&#39;&#39;&#39;
\ a b \ \ x\ = \ e \
\ c d \ \ y\ = \ f \
&#39;&#39;&#39;
e = fun(0,val[0,t],val[1,t])
f = fun(1,val[0,t],val[1,t])
div = 1/(a*d - b*c)#二阶矩阵求逆公式
val[0,t+1] = val[0,t] - (d*e - b*f)*div
val[1,t+1] = val[1,t] - (-c*e + a*f)*div扔到unity中可视化如下，红色曲面是第一个方程，蓝色曲面是第二个方程，在红色曲面高度与蓝色曲面高度都为零的地方即为根。为了更好地可视化，我将牛顿迭代法的速度放慢了。
这份可视化的代码在clatterrr/FluidSimulationTutorialsUnity。有了牛顿迭代法这件趁手的兵器，我们就可以上路去讨伐非线性偏微分方程了。不过请注意，由于我没有找到我能。有了牛顿迭代法这件趁手的兵器，我们就可以上路去讨伐非线性偏微分方程了。不过请注意，由于我没有找到我能看懂的代码，所以以下代码都是我基于自己的理解写的，可能存在错误。
第一关

现在需要从最简单的非线性方程算起。NavierStokes方程太复杂，甚至无粘性Burgers方程也不是最简单的。真正最基础的非线性方程如下，其中u是待求解的变量，x是坐标轴。d1d2是已知常数。

$u\frac{\partial U}{\partial x} = d_1x + d_2 \tag{12}$
假设u = c1x + cx，那么上式其实就是

$(c_1x + c_2)\frac{\partial (c_1x + c_2)}{\partial x} = c_1^2x + c_1c_2 = d_1x + d_2 \tag{13}$
权函数就用x和1-x。用Galerkin法展开

$\int \phi U\frac{\partial U}{\partial x} = d_1\int \phi x + d_2\int\phi \tag{14}$
等式右侧好算，直接数值积分一波算完
def simpon(p0,p1,p2,h):
return h/3*(p0 + 4*p1 + p2)
for i in range(0,num):
rhs[0,i] = (phi(i,0) + phi(i,1))/2 #梯形公式算数值积分 , 权函数本身的积分
for xr in range(0,3):
      x0 = xr / 2
      inte[xr] = x0*phi(i,x0) # 权函数乘上x的积分
rhs[1,i] = simpon(inte[0], inte[1], inte[2], 0.5)等式左侧，化成下面这个样子

$\int \begin{bmatrix} \phi_1 \\ \phi_2 \\ ..\end{bmatrix}\begin{bmatrix} \phi_1 & \phi_2 & ..\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ ..\end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} & ..\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ ..\end{bmatrix} \\ = \int \begin{bmatrix} \phi_1 \\ \phi_2 \\ ..\end{bmatrix}\begin{bmatrix} \phi_1 & \phi_2 & ..\end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} & ..\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ ..\end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \\ ..\end{bmatrix} \tag{15}$
鉴于现在只有两个权函数，所以最后得到的形式很简单

$\int \begin{bmatrix} \phi_1 \\ \phi_2 \end{bmatrix}\begin{bmatrix} \phi_1 & \phi_2 \end{bmatrix}\begin{bmatrix} \phi_{1x} & \phi_{2x} \end{bmatrix}\begin{bmatrix} u_1 \\ u_2 \end{bmatrix}\begin{bmatrix} u_1 \\ u_2\end{bmatrix} \\ = \int \begin{bmatrix} \phi_1 \\ \phi_2 \end{bmatrix}\begin{bmatrix} \phi_1\phi_{1x} & \phi_1\phi_{2x} & \phi_2\phi_{1x} & \phi_2\phi_{2x}\end{bmatrix}\begin{bmatrix} u_1^2 \\ u_1u_2 \\ u_2u_1 \\ u_2^2 \end{bmatrix} \tag{16}$
上面那个横着的一行四列的矩阵，乘上后面一列四行的矩阵，最终仅仅是等于一个数而已。不过那个一列四行的矩阵里面都是未知量，所以先把那些已知的权函数积分了。写成代码如下
for k in range(0,num):# num是权函的数量，是2
for i in range(0,num):
      for j in range(0,num):
         for xr in range(0,3):
            x0 = xr / 2
            inte[xr] = phi(k,x0)*phi(i,x0)*(phi(j,x0+1) - phi(j,x0-1))/2
         Fcoeff[k,i,j] = simpon(inte[0], inte[1], inte[2], 0.5)Fcoeff变量有三个维度。不同的k代表方程组中不同的方程。而i和j维度，代表的就是那个一行四列的矩阵。组装矩阵，牛顿迭代法中的F矩阵计算如下
for p in range(0,num):
for i in range(0,num):
      for j in range(0,num):
         k0 = int(idxk[p])# 求解的位置
         i0 = int(idxk)# u1的索引
         j0 = int(idxk[j])# u2的索引
         Fmat[k0] = Fmat[k0] + Fcoeff[p,i,j]*res[t,i0]*res[t,j0]Fcoeff这么写还有个好处，就是算雅可比矩阵很简单。

$J = \begin{bmatrix} \partial F_1 / \partial u_1 & \partial F_1 / \partial u_2 \\\partial F_2 / \partial u_1 & \partial F_2 / \partial u_2 \end{bmatrix} \tag{17}$
雅可比矩阵第一行第一列展开如下

$F_1 = \int\phi_1(\phi_1\phi_{1x}u_1^2 + \phi_1\phi_{2x}u_1u_2 + \phi_2\phi_{1x}u_1u_2 + \phi_2\phi_{2x}u_2^2) \\ \partial F_1 / \partial u_1 = \int\phi_1(2\phi_1\phi_{1x}u_1 + \phi_1\phi_{2x}u_2 + \phi_2\phi_{1x}u_2) \tag{18}$
当下面代码i = 0且 j = 0的时候，算的就是雅可比矩阵第一行第一列，也就是上面这个式子。
for i in range(0,num):
for j in range(0,num):
      i0 = int(idxk)
      j0 = int(idxk[j])
      for p in range(0,num):
         for q in range(0,num):
            if p == j:
                  Jmat[i0,j0] += Fcoeff[i,p,q]*res[t,int(idxk[q])]
            if q == j:
                  Jmat[i0,j0] += Fcoeff[i,p,q]*res[t,int(idxk[p])]在代码中，雅可比矩阵第一列第一行还需要继续展开成p行q列的矩阵。

$F_1 = \begin{bmatrix} \int\phi_1 \phi_1\phi_{1x}u_1^2 & \int\phi_1 \phi_1\phi_{2x}u_1u_2 \\ \int\phi_1 \phi_2\phi_{1x}u_1u_2 & \int\phi_1 \phi_2\phi_{2x}u_2^2 \end{bmatrix} \tag{19}$
Fcoeff中k = 0时就是上面这个权函积分结果
把F1求导得到

$\partial F_1 /\partial u_1= \begin{bmatrix} 2\int\phi_1 \phi_1\phi_{1x}u_1 & \int\phi_1 \phi_1\phi_{2x}u_2 \\ \int\phi_1 \phi_2\phi_{1x}u_2 & 0 \end{bmatrix} \tag{20}$
我的代码中就是用p和q去遍历这个矩阵。而且我用的if p == j和if q == j，而不是else if，就是为了方便算左上角的x^2 = x + x。
然后值得一写的就是非线性有限元如何组装矩阵。对于一开始求解偏微分方程的左边项，直接相加就行了。但对于右边项有点特殊。比如最终结果是2x + 1，组装矩阵时那么，对于0处和1处的值应该是1和3，这没错。但对于1和2处，结果应该是3和5，但此时1又恰巧是第一个，这么一算难道不应该是2 *0 + 1 = 1么？注意非线性有限元的刚度矩阵不需要初始条件也能求解。所以不能通过给初始条件解决这个问题。
不过解法也很简单，对于0和1处的它的结果应该是2x+1，而对于1和2处的结果应该是2(x+1) + 1，这样1处作为组装矩阵中的第一个，也就是将0代入2(x+1) + 1得到3就能准确无误了。
for i in range(0,num):
idx = int(idxk)
tempd1 = d1
tempd2 = d2 + d1 * idxk[0]
Fmat[idx] = Fmat[idx]  - tempd2*rhs[0,i] -  tempd1*rhs[1,i]最后，牛顿迭代法要求初始随便猜一个值。但最好别全部猜零，否则雅可比矩阵全为零是无法求逆的。最后结果如下，横轴是x轴，竖轴是牛顿迭代法的迭代次数，这时候牛顿法还是很快很稳定很准确的。这里完整代码为nonlinear01.py。
当然了，初学者写这些迭代算法的时候，写着写着就就发散了。不过也有一些特殊的方法来判断自己程序哪里出错，就像上篇说的那样，就把自己知道的最终稳定的结果当成初始条件代进去，那么正确的牛顿迭代法肯定算得之后的变化量都为零。借此可确定自己程序究竟哪里有问题。
第二关

闯关一次函数的关卡了，现在来挑战二次函数。设U的解形式如下

$U = c_1x^2 + c_2x + c_3 \tag{21}$
那么

$U\frac{\partial U}{\partial x} = (2c_1^2) x^3 + (3c_1c_2)x^2 + (2c_1c_3 + c_2^2) x + c_2c_3 \tag{22}$
其大体步骤与之前的差不多。主要变化由于二次函数需要用到三个权函数，为组装矩阵时一次组装三点
for k in range(0,nmax - 2):
idxk[0] = k
idxk[1] = k + 1
idxk[2] = k + 2
for i in range(0,num):
      idx = int(idxk)
      td1 = d1
      td2 = d1*3*k + d2
      td3 = d1*3*k*k + d2*2*k + d3
      td4 = d1*k*k*k + d2*k*k + d3*k + d4
      Fmat[idx] = Fmat[idx]  - td1*rhs[0,i] - td2*rhs[1,i] - td3*rhs[2,i] - td4*rhs[3,i]上面这么写的依据是

$d_1(x+k)^3 + d_2(x+k)^2 + d_3(x+k) + d_4 \\ = d_1x^3 + (3d_1k + d_2)x^2 + (3d_1k^2 + 2d_2k + d_3)x + (d_1k^3 + d_2k^2 + d_3k + d_4) \tag{23}$
另外这里如果不设置初始条件的话，最终可能算出错误的结果。看看之前的那个方程，里面的未知量组成都是两个未知量相乘，这说明肯定有每个未知量肯定有两个不同解，牛顿迭代法就会选择一个容易得到的解返回。

$F_1 = \begin{bmatrix} \int\phi_1 \phi_1\phi_{1x}u_1^2 & \int\phi_1 \phi_1\phi_{2x}u_1u_2 \\ \int\phi_1 \phi_2\phi_{1x}u_1u_2 & \int\phi_1 \phi_2\phi_{2x}u_2^2 \end{bmatrix} \tag{24}$
然而事实上我们要解的偏微分方程只有一个解，如果得到另一个错误的解，并不是牛顿迭代法的锅。因此我们一开始就规定至少一个初始条件，就可最终得到正确的解。代码写起来如下。
for i in range(0,nmax): # 循环里是初始猜值，随便猜，只要不全是零就行
res[0,i] = i*i*i + 5
res[:,0] = c1*(-1)*(-1) + c2*(-1) + c3 # 初始条件，因为权函数是从-1积分起的
... # 下面是牛顿迭代中
# J du = F，我们希望du_0 = 0，也就是u_0一直不变，因为这是在前面设置好的初始条件
Jmat[0,:] = 0
Jmat[0,0] = 1
Fmat[0] = Fmat2[0] = 0最后求解过程如下，很快很稳定地求了出来。设c1 = 3,c2 = 5,c3 = 2。注意横轴是x轴，是从-1开始的
如果不设定至少一个初始条件的话，牛顿迭代法也能收敛。但除非初始猜测的值和实际的值很接近，否则就会返回错误的结果。下面这样的结果是错误的，虽然它离正确答案很接近，但到目前为止我们没有用任何近似算法，所以也不会允许任何误差，所以迭代成这样的程序肯定是错的。完整代码为nonlinear02.py。
第三关

来给我们的偏微分方程加上时间变量。如下

$U = c_1x + c_2t + c_3 \tag{25}$
那么就可以写出一个类似无粘性Burgers方程的偏微分方程。

$\frac{\partial U}{\partial t} + U\frac{\partial U}{\partial x} = c_0^2x + c_0c_1t + c_0c_2 + c_1 \tag{26}$
真正的非粘性Burgers方程等式右边是零，此时U的解析解非常复杂。于是我只好先把它化成一个比较简单的问题。与之前的算法第一处不同是组装矩阵时

$d_0(x + x_0) + d_1(t + t_0) + d_2 = (d_0)x + (d_1)t + (d_0x_0 + d_1t_0 + d_2) \tag{27}$
与之前不同的就是U有了两个维度，一个x轴，一个t轴。从一个维度到两个维度，最主要就是要注意网格编号。我的是习惯是先遍历x轴再遍历t轴，因此一维数据中，第一个数据t=0,x=0，第二个t=0,x=1，第三个t=0,x=2....第nmax-1个，t = 0,x = nmax - 1，第nmax个，t = 1,x = 0，第nmax + 1个t = 1,x = 0如此循环下去。从我的权函数顺序和组装矩阵时的顺序也可以看出来。
def phi(n,t,x):
if n == 0:
      return (1 - t)*(1 - x)
elif n == 1:
      return (1 - t)*x
elif n == 2:
      return t*(1 - x)
elif n == 3:
      return t*x
...
for k in range(0,tmax*nmax):
      colu = int(k / nmax) # t 轴
      row = k % nmax # x 轴
      idxk[0] = colu * nmax + row
      idxk[1] = idxk[0] + 1
      idxk[2] = idxk[0] + nmax
      idxk[3] = idxk[2] + 1
      if(colu >= tmax - 1)|(row >= nmax - 1):
         continue与仅有x轴不同，这次U还要求对时间的偏导。组装矩阵时别忘了。雅可比矩阵对未知量求导之后，就只剩下权函数积分的结果了。完整代码为nonlinear03.py。
for i in range(0,num):
for j in range(0,num):
      i0 = int(idxk)
      j0 = int(idxk[j])
      Fmat3[i0] = Fmat3[i0] + Fdt[i,j]*res[t,j0]
      Jmat[i0,j0] += Fdt[i,j]第四关

本章最后，假设U具有如下形式

$U= c_0x^2 + c_1x + c_2t^2 + c_3t + c_4 \tag{28}$
那么下面的偏微分方程解如下

$\frac{\partial U}{\partial t} + U\frac{\partial U}{\partial x} = 2c_0^2x^3 + 3c_0c_1x^2 + (2c_0x_4 + c_1^2)x\\ + c_1c_2t^2 + (c_1c_3 + 2c_2)t + + 2c_0c_2xt^2 + 2c_0c_3xt + (c_1c_4 + c_4) \tag{29}$
或者是加上粘性项，仅仅是多了个常数

$\frac{\partial U}{\partial t} + U\frac{\partial U}{\partial x} - \mu \frac{\partial^2 U}{\partial x^2}= 2c_0^2x^3 + 3c_0c_1x^2 + (2c_0x_4 + c_1^2)x\\ + c_1c_2t^2 + (c_1c_3 + 2c_2)t + + 2c_0c_2xt^2 + 2c_0c_3xt + (c_1c_4 + c_4 - 2\mu c_0) \tag{30}$
代码文件分别为nonlinear04.py和nonlinear05.py。真正的Burgers等式右边应该是零，然而这意味着U的解析解非常复杂，不大可能用二次函数去拟合。因此这算是一个非常简化的偏微分方程了。
而且在我的代码中，要保证牛顿迭代法最终收敛，要么一开始猜的值足够接近真实值，要么给出的限制条件足够多。我认为这是牛顿迭代法有某些特殊的性质我没发现导致的。并且还要在每步迭代时，给du乘上一个非常小的值比如0.01，控制U的变化量，防止迭代法在解的附近迈的步子太大回不来。而且我用的网格非常少，比如nonlinear04.py只有8x8个网格。
我暂时写不出更复杂的代码了，但我也不打算在此纠结，第一是即使这样，我也能用非线性有限元问题也能解决很多问题了，很多结构力学，电磁学中的非线性偏微分方程，与上面这个一维粘性Burgers方程是很相似的。第二是我只用了牛顿法这一种方法，等之后我会尝试拟牛顿法，非精确牛顿法等，借此或许能更多地认识牛顿法本身。第三是我倒是很怕自己陷入某种思维定势，认为只有这一种方法解非线性有限元了。不如从现在起把它忘了，再过一年再来看肯定有很多新认识的。第四是除了非线性有限元外，还有很多有趣的东西。
有关非线性方程组数值解法的书有很多很不错的书，不过我参考的是《迭代方法和预处理技术》谷同祥，安恒斌等编著。好了，有限元之旅暂时到此。下一篇是频谱法模拟水波。敬请期待。
上一篇，光影帽子：【游戏流体力学基础及Unity代码（十五）】线性有限元及弹性物体模拟
搜索代码的技巧

除了在搜索引擎还有github上外，还有一些不为人知的源代码搜索技巧。某些Springer出版的书籍会附一些代码以及其它一些资料，可以在https://extras.springer.com/这个网站免费下到。按照提示找即可。比如《MATLAB Codes for Finite Element Analysis》上有一些代码，但照着敲太麻烦且意义不大，于是可以先去Springer上找到本书的信息，然后确定它的年份，以及它的实体书(Softcover或者Hardcover)ISBN。这本书是2009年的，那么就可以打开https://extras.springer.com/2009/，然后在一大堆ISBN号中找到它，点击下载就行了。
以及我一般如果觉得某些资料不错的话，会去搜索这本书作者的个人主页，一般他都会是某学校的教授，又正在讲课或曾经讲过课，于是会有一些讲义附在个人主页上可免费下载，一般Teaching一栏上。某些讲义质量非常高，但又很难通过搜索引擎搜到。
如果运气好的话，这位德高望重的教授还会附一些简洁易懂的教学源代码在他的个人主页上，而且是独家的未上传过github的代码。偶尔还能找到一些已发表论文的开源代码。
而且既然这位教授技术如此高超且非常慷慨，那么想必他的同事也...嘿嘿
同理，我在github上看到不错的代码的话，star之余也肯定会点开作者的github主页，看看他还有没有其它的不错代码。这三种方法经常能搜到一些不错的代码，不过这些代码能不能运行起来又是另一回事了。不过有总比没有好。

BlaXuan · 发表于 2021-5-9 21:39

老弟真的厉害！可以可以！

		自动登录	找回密码
密码			立即注册

[笔记] 【游戏流体力学基础及Unity代码（十六）】非线性有限元及牛顿迭代法

本帖子中包含更多资源

浏览过的版块