本文主要参考【1】。
我大学时候大物和工数都挂过,如果你觉得我说的有不对的地方,那多半是我真的不对,欢迎真正的流体力学大佬指正。
若有方程组形式如下:
![equation?tex=%5Cbegin%7Bcases%7D+ax%2Bby%3Dc%5C%5C+dx%2Bey%3Df+%5Cend%7Bcases%7D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/637da12cfed59544fb4e29a4e604f9ca.png)
那么我们可以尝试把它写成矩阵形式:
![equation?tex=%5Cleft%5B++%5Cbegin%7Barray%7D%7Bccc%7D+a+%26+b%5C%5C+d+%26+e%5C%5C+%5Cend%7Barray%7D+%5Cright%5D+%5Cleft%5B++%5Cbegin%7Barray%7D%7Bccc%7D+x%5C%5C+y%5C%5C+%5Cend%7Barray%7D+%5Cright%5D%3D+%5Cleft%5B++%5Cbegin%7Barray%7D%7Bccc%7D+c%5C%5C+f%5C%5C+%5Cend%7Barray%7D+%5Cright%5D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/781635c7804eaca5198718c50316eef4.png)
或者说
![equation?tex=+%5Cbold+A%5Cbold+x%3D%5Cbold+b](https://i-blog.csdnimg.cn/blog_migrate/7454333fb4951b70ffba5f39cd9f1455.png)
![equation?tex=+%5Cbold+x%3D%5Cbold+A%5E%7B-1%7D%5Cbold+b](https://i-blog.csdnimg.cn/blog_migrate/84d58568fb8d98cecaa9215ad2b8418f.png)
雅可比法
在数值线性代数中,雅可比法(Jacobi Method)是一种解对角元素几乎都是各行和各列的绝对值最大的值的线性方程组的算法。求解出每个对角元素并插入近似值。不断迭代直至收敛。这个算法是雅可比矩阵的精简版。方法的名字来源于德国数学家卡尔·雅可比。
给定一个n×n的线性方程组
![equation?tex=%7B%5Cdisplaystyle+A%5Cmathbf+%7Bx%7D+%3D%5Cmathbf+%7Bb%7D+%7D](https://i-blog.csdnimg.cn/blog_migrate/edc8e9e36a2dc403ef49fbb1fe41f89a.png)
![equation?tex=%7B%5Cdisplaystyle+A%3D%7B%5Cbegin%7Bbmatrix%7Da_%7B11%7D%26a_%7B12%7D%26%5Ccdots+%26a_%7B1n%7D%5C%5Ca_%7B21%7D%26a_%7B22%7D%26%5Ccdots+%26a_%7B2n%7D%5C%5C%5Cvdots+%26%5Cvdots+%26%5Cddots+%26%5Cvdots+%5C%5Ca_%7Bn1%7D%26a_%7Bn2%7D%26%5Ccdots+%26a_%7Bnn%7D%5Cend%7Bbmatrix%7D%7D%2C%5Cqquad+%5Cmathbf+%7Bx%7D+%3D%7B%5Cbegin%7Bbmatrix%7Dx_%7B1%7D%5C%5Cx_%7B2%7D%5C%5C%5Cvdots+%5C%5Cx_%7Bn%7D%5Cend%7Bbmatrix%7D%7D%2C%5Cqquad+%5Cmathbf+%7Bb%7D+%3D%7B%5Cbegin%7Bbmatrix%7Db_%7B1%7D%5C%5Cb_%7B2%7D%5C%5C%5Cvdots+%5C%5Cb_%7Bn%7D%5Cend%7Bbmatrix%7D%7D.%7D](https://i-blog.csdnimg.cn/blog_migrate/a459438f98485af0aa2453e7e0937aa3.png)
![equation?tex=%5Cbold+A](https://i-blog.csdnimg.cn/blog_migrate/9a5a03eb7d2a1b053d8b2fff4f4e5d18.png)
![equation?tex=%5Cbold+D](https://i-blog.csdnimg.cn/blog_migrate/9be3fafb8a03574c175f37f0e54c31eb.png)
![equation?tex=%5Cbold+R](https://i-blog.csdnimg.cn/blog_migrate/176cf33d3955875da86af54c9b4967f1.png)
![equation?tex=%7B%5Cdisplaystyle+%5Cbold+A%3D%5Cbold+D%2B%5Cbold+R%5Cqquad+%7B%5Ctext%7B%E5%85%B6%E4%B8%AD%7D%7D%5Cqquad+%5Cbold+D%3D%7B%5Cbegin%7Bbmatrix%7Da_%7B11%7D%260%26%5Ccdots+%260%5C%5C0%26a_%7B22%7D%26%5Ccdots+%260%5C%5C%5Cvdots+%26%5Cvdots+%26%5Cddots+%26%5Cvdots+%5C%5C0%260%26%5Ccdots+%26a_%7Bnn%7D%5Cend%7Bbmatrix%7D%7D%2C%5Cqquad+%5Cbold+R%3D%7B%5Cbegin%7Bbmatrix%7D0%26a_%7B12%7D%26%5Ccdots+%26a_%7B1n%7D%5C%5Ca_%7B21%7D%260%26%5Ccdots+%26a_%7B2n%7D%5C%5C%5Cvdots+%26%5Cvdots+%26%5Cddots+%26%5Cvdots+%5C%5Ca_%7Bn1%7D%26a_%7Bn2%7D%26%5Ccdots+%260%5Cend%7Bbmatrix%7D%7D%7D](https://i-blog.csdnimg.cn/blog_migrate/0f4d7dd2540e138da6294014f804c0c1.png)
线性方程组可以重写为:
![equation?tex=%7B%5Cdisplaystyle+%5Cbold++D%5Cmathbf+%7Bx%7D+%3D%5Cmathbf+%7Bb%7D+-+%5Cbold+R%5Cmathbf+%7Bx%7D+%7D+%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/45d9b97cb6b7b9391ecb6ff863675dcf.png)
雅可比法是一种迭代方法。在每一次迭代中,上一次算出的x被用在右侧,用来算出左侧的新的x。这个过程可以如下表示:
![equation?tex=%7B%5Cdisplaystyle+%5Cmathbf+%7Bx%7D+%5E%7B%28k%2B1%29%7D%3D%5Cbold+D%5E%7B-1%7D%28%5Cmathbf+%7Bb%7D+-%5Cbold+R%5Cmathbf+%7Bx%7D+%5E%7B%28k%29%7D%29%7D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/8d53e51af99ba0e0ac360d6359bb1c0c.png)
对每个元素可以用以下公式:
![equation?tex=%7B%5Cdisplaystyle+x_%7Bi%7D%5E%7B%28k%2B1%29%7D%3D%7B%5Cfrac+%7B1%7D%7Ba_%7Bii%7D%7D%7D%5Cleft%28b_%7Bi%7D-%5Csum+_%7Bj%5Cneq+i%7Da_%7Bij%7Dx_%7Bj%7D%5E%7B%28k%29%7D%5Cright%29%2C%5Cquad+i%3D1%2C2%2C%5Cldots+%2Cn.%7D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/1400df266586281f8cdb8bf2180e12b0.png)
注意计算
![equation?tex=%5Cbold+x%5E%7B%28k%2B1%29%7D](https://i-blog.csdnimg.cn/blog_migrate/a6b90bbd55d8123d2b8e9469f29f3065.png)
![equation?tex=%5Cbold+x%5E%7Bk%7D](https://i-blog.csdnimg.cn/blog_migrate/bb2bbb14df6287dd90188e47476b0c2e.png)
![equation?tex=x_i%5E%7Bk%2B1%7D](https://i-blog.csdnimg.cn/blog_migrate/3ece40eaeb058577d03e394a3c118047.png)
![equation?tex=x_i%5Ek](https://i-blog.csdnimg.cn/blog_migrate/ddea4d23b5cfba6969d94512ae7c8b0d.png)
一般来说,雅可比法可以通过不断的迭代提升解的精度。但是并非绝对,有时候雅可比法反而会产生一个发散的迭代序列。保证收敛的条件是矩阵
![equation?tex=%5Cbold+A](https://i-blog.csdnimg.cn/blog_migrate/9a5a03eb7d2a1b053d8b2fff4f4e5d18.png)
![equation?tex=%7B%5Cdisplaystyle+%5Cleft%7Ca_%7Bii%7D%5Cright%7C%3E%5Csum+_%7Bj%5Cneq+i%7D%7B%5Cleft%7Ca_%7Bij%7D%5Cright%7C%7D.%7D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/a1a5db052aded2fef34ccb2c40a2a34a.png)
有时即使不满足此条件,雅可比法仍可收敛。
高斯-赛德尔法
高斯-赛德尔法是数值线性代数中的一个迭代法,可用来求出线性方程组解的近似值。该方法以卡尔·弗里德里希·高斯和路德维希·赛德尔命名。
跟雅可比法类似,但不同的地方在于,求解过程中,所有已解出的k+1步元素都被直接使用,对于每个元素可以使用如下公式:
![equation?tex=%7B%5Cdisplaystyle+x_%7Bm%7D%5E%7Bk%2B1%7D%3D%7B%5Cfrac+%7B1%7D%7Ba_%7Bmm%7D%7D%7D%5Cleft%28b_%7Bm%7D-%5Csum+_%7Bj%3D1%7D%5E%7Bm-1%7Da_%7Bmj%7D%5Ccdot+x_%7Bj%7D%5E%7Bk%2B1%7D-%5Csum+_%7Bj%3Dm%2B1%7D%5E%7Bn%7Da_%7Bmj%7D%5Ccdot+x_%7Bj%7D%5E%7Bk%7D%5Cright%29%2C%5Cquad+1%5Cleq+m%5Cleq+n.%7D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/4de45569d9c1c7e6984a788afe85d489.png)
为了保证该方法可以进行,主对角线元素
![equation?tex=a_%7Bmm%7D](https://i-blog.csdnimg.cn/blog_migrate/7414e53281b73e817329ea14154eb187.png)
论文内容
首先,Jos提出速度的N-S方程,以及在该速度场下的密度方程:
![equation?tex=%5Cfrac+%7B%5Cpartial+%5Cbold+u%7D%7B%5Cpartial+t%7D+%3D+-%28%5Cbold+u%5Ccdot+%E2%88%87%29%5Cbold+u+%2B+%5Cnu+%E2%88%87%5E2%5Cbold+u+%2B+%5Cbold+f%5C%5C+%5Cfrac+%7B%5Cpartial+%5Crho%7D%7B%5Cpartial+t%7D+%3D+-%28%5Cbold+u%5Ccdot+%E2%88%87%29%5Crho+%2B+%5Ckappa+%E2%88%87%5E2%5Crho+%2B+S%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/42c7af71ac406765ff43557dfa34a23f.png)
首先,将求解密度方程,这允许读者熟悉求解器的不同组成部分,稍后会把求解密度方程的想法转移到相对更难的速度方程上。
![11de547d154a28de3518581b8b555901.png](https://i-blog.csdnimg.cn/blog_migrate/6d0894f551b01692961c8c91b64c73d4.png)
static u[size], v[size], u_prev[size], v_prev[size];
static dens[size], dens_prev[size];
//速度和密度的一维数组,通过如下形式使用
#define IX(i,j) ((i)+(N+2)*(j))
论文设计的求解器原型中,用户添加用鼠标添加密度源。假设相邻快照间的时间步长为dt。整个密度求解器的步骤如下:
![a7287cd352ce295d7e6c81cbb9889765.png](https://i-blog.csdnimg.cn/blog_migrate/de0495c6540c38116210c5ba50446f53.png)
密度求解器
先描述密度场在不随时间变化的固定速度场中移动的求解器。密度方程式表明,单个时间步长上的密度变化由三个原因引起。这些原因是方程中等号右侧的三项:第一项表明密度变化应跟随速度场;第二项表明密度以一定速度扩散;第三项表明密度由于密度源的存在而增大。求解器将按相反顺序解析这些项,我们可以称之为添加密度源、扩散、移动。
密度源
第一项即密度源,很容易实现,一般会通过游戏引擎的某些功能填充,本例中由数组s[]记录用户鼠标移动来填充。
void add_source ( int N, float * x, float * s, float dt)
{
int i, size=(N+2)*(N+2);
for ( i=0 ; i<size ; i++ )
x[i] += dt*s[i];
}
扩散
![f60ba2afa43302db5943ab3f3cd4d95e.png](https://i-blog.csdnimg.cn/blog_migrate/591143e47761ff6f9b9717372fe60a45.png)
第二项考虑到了可能的以
![equation?tex=diff](https://i-blog.csdnimg.cn/blog_migrate/79121949c1f63692d72b739ce639a4ef.png)
![equation?tex=diff%3E0](https://i-blog.csdnimg.cn/blog_migrate/79121949c1f63692d72b739ce639a4ef.png%3E0)
x0[IX(i-1,j)] + x0[IX(i+1,j)] + x0[IX(i,j-1)] + x0[IX(i,j+1)] - 4*x0[IX(i,j)]
扩散的简单求解方法,就是在每个格子里单独计算这些交换并更新现有值:
void diffuse_bad ( int N, int b, float * x, float * x0, float diff, float dt ) {
int i, j;
float a= dt * diff * N * N;
for ( i=1 ; i<=N ; i++ ) {
for ( j=1 ; j<=N ; j++ ) {
x[IX(i,j)] = x0[IX(i,j)] + a*(x0[IX(i-1,j)] + x0[IX(i+1,j)] +
x0[IX(i,j-1)] + x0[IX(i,j+1)] - 4*x0[IX(i,j)]);
}
}
set_bnd ( N, b, x );
}
set_bnd()正是稍后讨论的边界格子。但不幸的是,如果针对一个较大的负数值(指
![equation?tex=diff](https://i-blog.csdnimg.cn/blog_migrate/79121949c1f63692d72b739ce639a4ef.png)
x0[IX(i,j)] = x[IX(i,j)] - a * (x[IX(i-1,j)] + x[IX(i+1,j)] + x[IX(i,j-1)] + x[IX(i,j+1)] - 4 * x[IX(i,j)]);
这是由未知数x[IX(i,j)]构成的线性系统,我们可以为该系统建立矩阵然后调用标准矩阵求逆方法。但是由于矩阵过于稀疏,可以使用更简单的迭代技术来反转,比如高斯-赛德尔方法:
void diffuse ( int N, int b, float * x, float * x0, float diff, float dt )
{
int i, j, k;
float a = dt * diff * N * N;
for ( k=0 ; k<20 ; k++ ) {
for ( i=1 ; i<=N ; i++ ) {
for ( j=1 ; j<=N ; j++ ) {
x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+
x[IX(i,j-1)] + x[IX(i,j+1)]))/(1+4*a);
}
}
set_bnd ( N, b, x );
}
}
移动
与扩散步骤类似,可以建立一个线性系统,然后用迭代方法搞定。但是这次生成的线性方程将与速度有关,因此求解很棘手。
考虑将密度建模成一组粒子可能会更加有效。这种情况下只需要通过速度场去推动粒子。假定每个网格单元的中心是一个粒子,使用速度场对其进行推动,然后再将粒子转化为网格的值,如下图(b)所示。为了正确取得网格中心的密度值,并不需要对网格中心的粒子进行推动,而是找到推动之后恰好处于网格中心的粒子,它携带的密度值可以通过四个邻居的插值进行计算,如下图(c)所示。所以这是要维护两个网格数组的缘故:
![473f161c08635e7fd379f6ce23d1b197.png](https://i-blog.csdnimg.cn/blog_migrate/09508775f5b1ec510ef7c735b9fc3861.png)
![a513797e9ee8471e9439a729b9c6b406.png](https://i-blog.csdnimg.cn/blog_migrate/160d99b4f591b3b4e6368648e365b71c.png)
![b6626236a551dab6cf8af62bab66144d.png](https://i-blog.csdnimg.cn/blog_migrate/4542f418c16d09480a70138fe1108227.png)
假定最开始的密度源包含在x0数组中,则一个时间步长内的变化为:
![6b781e298ad68db49101fe96c995dac1.png](https://i-blog.csdnimg.cn/blog_migrate/e2cb8c8f14ecab01107f7a4adcc83528.png)
SWAP为交换数组指针的宏。
一个小小的实践
【1】中举了一个小小的例子作为两种迭代法的试验场。
假定有一根棍子,其左端被恒定为单位温度0,右端被恒定为单位温度1。最终棍子上的温度分布将趋于稳定,这可以被视作一维下的扩散。
将棍子均匀划分为8个点,每个点的温度被记为
![equation?tex=x_i](https://i-blog.csdnimg.cn/blog_migrate/55abd6ee2c07065bfd8875b4e79a1cf9.png)
![equation?tex=x_1%3D0%2Cx_8%3D1](https://i-blog.csdnimg.cn/blog_migrate/35102e3a502d9be731b4616040d9c981.png)
![equation?tex=x_i%27%3Dx_%7Bi-1%7D%2Bx_%7Bi%2B1%7D-2x_i%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/55abd6ee2c07065bfd8875b4e79a1cf9.png%27%3Dx_%7Bi-1%7D%2Bx_%7Bi%2B1%7D-2x_i%5C%5C)
![equation?tex=x%27_i](https://i-blog.csdnimg.cn/blog_migrate/d2e2a5143af317c386310f9c39f50ee4.png)
![equation?tex=%5Cbegin%7Bcases%7Dx_2%27%3Dx_%7B1%7D%2Bx_%7B3%7D-2x_2%3D0%2Bx_%7B3%7D-2x_2%5C%5C+x_3%27%3Dx_%7B2%7D%2Bx_%7B4%7D-2x_3%5C%5C+x_4%27%3Dx_%7B3%7D%2Bx_%7B5%7D-2x_4%5C%5C+x_5%27%3Dx_%7B4%7D%2Bx_%7B6%7D-2x_5%5C%5C+x_6%27%3Dx_%7B5%7D%2Bx_%7B7%7D-2x_6%5C%5C+x_7%27%3Dx_%7B6%7D%2Bx_%7B8%7D-2x_7%3Dx_%7B6%7D%2B1-2x_7%5C%5C+%5Cend%7Bcases%7D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/c63f3bcc8a1ba11b289a6e28df7a5534.png)
写成矩阵就是:
![equation?tex=%5Cleft%5B%5Cbegin%7Barray%7D%7Bccc%7D+-2+%26+1+%26+0+%26+0+%26+0+%26+0%5C%5C+1+%26+-2+%26+1+%26+0+%26+0+%26+0%5C%5C+0+%26+1+%26+-2+%26+1+%26+0+%26+0%5C%5C+0+%26+0+%26+1+%26+-2+%26+1+%26+0%5C%5C+0+%26+0+%26+0+%26+1+%26+-2+%26+1%5C%5C+0+%26+0+%26+0+%26+0+%26+1+%26+-2%5C%5C+%5Cend%7Barray%7D%5Cright%5D+%5Cleft%5B%5Cbegin%7Barray%7D%7Bccc%7D+x_2%5C%5Cx_3%5C%5Cx_4%5C%5Cx_5%5C%5Cx_6%5C%5Cx_7+%5Cend%7Barray%7D%5Cright%5D%3D+%5Cleft%5B%5Cbegin%7Barray%7D%7Bccc%7D+0%5C%5C0%5C%5C0%5C%5C0%5C%5C0%5C%5C-1%5C%5C+%5Cend%7Barray%7D%5Cright%5D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/c9a01527153d6c30e7266c12043a2978.png)
可见这是一个规律且稀疏的矩阵,我们可以将其拆成更加简便的三个矩阵:
![equation?tex=%5Cleft%5B%5Cbegin%7Barray%7D%7Bccc%7D+-2+%26+1+%26+0+%26+0+%26+0+%26+0%5C%5C+1+%26+-2+%26+1+%26+0+%26+0+%26+0%5C%5C+0+%26+1+%26+-2+%26+1+%26+0+%26+0%5C%5C+0+%26+0+%26+1+%26+-2+%26+1+%26+0%5C%5C+0+%26+0+%26+0+%26+1+%26+-2+%26+1%5C%5C+0+%26+0+%26+0+%26+0+%26+1+%26+-2%5C%5C+%5Cend%7Barray%7D%5Cright%5D%3D+%5Cleft%5B%5Cbegin%7Barray%7D%7Bccc%7D+0+%26+1+%26+0+%26+0+%26+0+%26+0%5C%5C+0%26+0+%26+1+%26+0+%26+0+%26+0%5C%5C+0+%26+0+%26+0+%26+1+%26+0+%26+0%5C%5C+0+%26+0+%26+0+%26+0+%26+1+%26+0%5C%5C+0+%26+0+%26+0+%26+0+%26+0+%26+1%5C%5C+0+%26+0+%26+0+%26+0+%26+0%26+0%5C%5C+%5Cend%7Barray%7D%5Cright%5D%2B+%5Cleft%5B%5Cbegin%7Barray%7D%7Bccc%7D+-2+%26+0+%26+0+%26+0+%26+0+%26+0%5C%5C+0+%26+-2+%26+0+%26+0+%26+0+%26+0%5C%5C+0+%26+0+%26+-2+%26+0+%26+0+%26+0%5C%5C+0+%26+0+%26+0+%26+-2+%26+0+%26+0%5C%5C+0+%26+0+%26+0+%26+0+%26+-2+%26+0%5C%5C+0+%26+0+%26+0+%26+0+%26+0+%26+-2%5C%5C+%5Cend%7Barray%7D%5Cright%5D%2B+%5Cleft%5B%5Cbegin%7Barray%7D%7Bccc%7D+0+%26+0+%26+0+%26+0+%26+0+%26+0%5C%5C+1+%26+0+%26+0+%26+0+%26+0+%26+0%5C%5C+0+%26+1+%26+0+%26+0+%26+0+%26+0%5C%5C+0+%26+0+%26+1+%26+0+%26+0+%26+0%5C%5C+0+%26+0+%26+0+%26+1+%26+0+%26+0%5C%5C+0+%26+0+%26+0+%26+0+%26+1+%26+0%5C%5C+%5Cend%7Barray%7D%5Cright%5D%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/90d294796835c97cfdc7de169e871e90.png)
可以凝练出一个符号表达:
![equation?tex=%5Cbold+A%3D%5Cbold+L%2B%5Cbold+D%2B%5Cbold+U](https://i-blog.csdnimg.cn/blog_migrate/9a5a03eb7d2a1b053d8b2fff4f4e5d18.png%3D%5Cbold+L%2B%5Cbold+D%2B%5Cbold+U)
雅可比迭代法
将变形过程表示为:
![equation?tex=%5Cbold+D%5Cbold+x%3D-%28%5Cbold+L%2B%5Cbold+U%29%5Cbold+x%2B%5Cbold+b%5C%5C+%5Cbold+x%3D-%5Cbold+D%5E%7B-1%7D%28%5Cbold+L%2B%5Cbold+U%29%5Cbold+x%2B%5Cbold+D%5E%7B-1%7D%5Cbold+b](https://i-blog.csdnimg.cn/blog_migrate/9be3fafb8a03574c175f37f0e54c31eb.png%5Cbold+x%3D-%28%5Cbold+L%2B%5Cbold+U%29%5Cbold+x%2B%5Cbold+b%5C%5C+%5Cbold+x%3D-%5Cbold+D%5E%7B-1%7D%28%5Cbold+L%2B%5Cbold+U%29%5Cbold+x%2B%5Cbold+D%5E%7B-1%7D%5Cbold+b)
于是对于
![equation?tex=%5Cbold+x](https://i-blog.csdnimg.cn/blog_migrate/1ff7b56bcb4c40192eb47906a05033fb.png)
![equation?tex=x_%7Bi%7D%3D-%5Cfrac%7B1%7D%7Ba_%7Bii%7D%7D%28b_i-%5Csum_%7Bi%21%3Dj%7D%5E%7B%7D%7Ba_%7Bij%7Dx%27_j%7D%29](https://i-blog.csdnimg.cn/blog_migrate/f151754f32ab62c214435c1e20df372f.png)
高斯-赛德尔法
雅可比迭代法为什么适用于并行计算呢?区别在于是否使用
![equation?tex=x_i%5E%7Bk%2B1%7D](https://i-blog.csdnimg.cn/blog_migrate/3ece40eaeb058577d03e394a3c118047.png)
![equation?tex=x_i%5Ek](https://i-blog.csdnimg.cn/blog_migrate/ddea4d23b5cfba6969d94512ae7c8b0d.png)
我们假定在初始时棍子温度均为0,相当于设定了一个初始解
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class Experiment : MonoBehaviour
{
float[] array_Jacobi, array_Jacobi_prev;
float[] array_GaussSeidel;
float[] rhs;
float[] matrix;
int times;
private void Start()
{
rhs = new float[] { 0, 0, 0, 0, 0, -1 };
matrix = new float[]
{
-2, 1, 0, 0, 0, 0,
1, -2, 1, 0, 0, 0,
0, 1, -2, 1, 0, 0,
0, 0, 1, -2, 1, 0,
0, 0, 0, 1, -2, 1,
0, 0, 0, 0, 1, -2,
};
array_Jacobi = new float[6];
array_Jacobi_prev = new float[6];
array_GaussSeidel = new float[6];
times = 0;
//
for(int i = 0; i < 6; i++)
{
array_Jacobi_prev[i] = 0;
array_GaussSeidel[i] = 0;
}
}
float sum(int i, float[] a, float[] x0)
{
float result = 0;
for (int j = 0; j < 6; j++)
{
if (j != i)
result += a[i * 6 + j] * x0[j];
}
return result;
}
void Update()
{
times++;
for (int i = 0; i < 6; i++)
{
array_Jacobi[i] = (rhs[i] - sum(i, matrix, array_Jacobi_prev)) / matrix[i * 6 + i];
array_GaussSeidel[i] = (rhs[i] - sum(i, matrix, array_GaussSeidel)) / matrix[i * 6 + i];
}
var t = array_Jacobi_prev;
array_Jacobi_prev = array_Jacobi;
array_Jacobi = t;
Debug.LogErrorFormat("Times:{0}", times);
Debug.LogFormat("{0} {1} {2} {3} {4} {5}", array_Jacobi[0], array_Jacobi[1], array_Jacobi[2], array_Jacobi[3], array_Jacobi[4], array_Jacobi[5]);
Debug.LogWarningFormat("{0} {1} {2} {3} {4} {5}", array_GaussSeidel[0], array_GaussSeidel[1], array_GaussSeidel[2], array_GaussSeidel[3], array_GaussSeidel[4], array_GaussSeidel[5]);
}
}
![dcd23918fa9adf6723254d6615d7169d.png](https://i-blog.csdnimg.cn/blog_migrate/3c1e02a9887fbdf4fced0b872c7df2f1.jpeg)
精度达到第1、2、3、4位有效数字,高斯-赛德尔迭代用了8次、18次、37次、43次;而雅可比迭代用了16次、36次、74次、86次。
实际上,迭代的过程,就是一维情况下密度源(温度分布)扩散的过程。
速度求解器
![equation?tex=%5Cfrac+%7B%5Cpartial+%5Cbold+u%7D%7B%5Cpartial+t%7D+%3D+-%28%5Cbold+u%5Ccdot+%E2%88%87%29%5Cbold+u+%2B+%5Cnu+%E2%88%87%5E2%5Cbold+u+%2B+%5Cbold+f%5C%5C+%5Cfrac+%7B%5Cpartial+%5Crho%7D%7B%5Cpartial+t%7D+%3D+-%28%5Cbold+u%5Ccdot+%E2%88%87%29%5Crho+%2B+%5Ckappa+%E2%88%87%5E2%5Crho+%2B+S%5C%5C](https://i-blog.csdnimg.cn/blog_migrate/42c7af71ac406765ff43557dfa34a23f.png)
之前的文章中,探讨了密度分布的变化,并提出了密度求解器的工作步骤。接下来,就轮到速度了。我们将速度方程解读为,速度随着时间步长的变化可归结于三个原因:力的增加、黏度扩散和自对流。自对流不太明显,可以简单的解释为速度场沿着自身移动方向移动。
正如同密度求解器那样,速度求解器亦可以如此:
![b5571cbfa68979ff07bb2f7f93f344a3.png](https://i-blog.csdnimg.cn/blog_migrate/b89e0adac4575bdae020273c9b6ec4b2.png)
在上述代码中,有一个名为project()的函数跟密度求解器不一致。该函数强迫速度保证质量守恒,这是被驱动的真实流体所具有的属性,如下图所示,也如同该系列第一篇专栏《【前置技能】流体力学的一些基础知识》所说的那样,质量守恒则带来无散场,无散场中拥有大量的涡流状图案,其视觉效果就是流体将具备丰富的涡流。
![21347cb9a3eec55cfd3f364b5ddaf615.png](https://i-blog.csdnimg.cn/blog_migrate/f79b95d50d829be9717ae42c8df8e971.jpeg)
以下是将速度投影到一个无散场的project()的代码:
![bd9c8e2809202394d690945d0d2be0aa.png](https://i-blog.csdnimg.cn/blog_migrate/8d9a2f1a7e489a6703d3b7e69434827a.png)
在流程中两次调用了project(),为了让advect()表现的更为精准。
边界条件
边界条件由set_bnd()函数处理,假定流体的边界是固体墙壁,意味着垂直于边界的速度水平分量为零。对于代码中考虑的密度和其他场,假定其是连续的:
![8b7049bbbe6813d8905e4d9c1536c6c7.png](https://i-blog.csdnimg.cn/blog_migrate/10ea476dcc8f93b68c5085f60a7f9e7d.png)
![3f5117c486c017aec84292e02e3ea112.png](https://i-blog.csdnimg.cn/blog_migrate/dd9d4210c610204e38ab57837567fc2e.png)
以下展示了各个部分的代码是如何组合并运作的:
![6fcb574fc5a6ec3f628fbe07e77ca3d1.png](https://i-blog.csdnimg.cn/blog_migrate/f15637a2cbd0bbd8900095fa6ff23891.png)
效果演示
Gif效果如下,想要学习请移步【6】:
![d29a10bb56d15a02a2f05c1ddb3554f6.gif](https://i-blog.csdnimg.cn/blog_migrate/62bb987db5128800a4e2945355feb9ef.gif)
![8cb35ecf1483c4acab8ac4c15e54b015.gif](https://i-blog.csdnimg.cn/blog_migrate/7ecde244934be541c106271f4297c005.gif)
总结
本文介绍了Jos的论文,并以一个简单的实验为例,说明了论文中关于密度求解器部分的内容,并用代码稍作演示;此外介绍了日本keijiro大神的一个实现。
参考文献
【1】Real-Time Fluid Dynamics for Games
【2】雅克比法 - 维基百科
【3】高斯-赛德尔法 - 维基百科
【4】https://blog.csdn.net/baimafujinji/article/details/50582462
【5】https://zhuanlan.zhihu.com/p/30965284
【6】keijiro/StableFluids