c++用类实现高斯消元法求解线性方程组的解_【前置技能】论文和其他人的实现...

本文主要参考【1】。

我大学时候大物和工数都挂过，如果你觉得我说的有不对的地方，那多半是我真的不对，欢迎真正的流体力学大佬指正。

---

若有方程组形式如下：

那么我们可以尝试把它写成矩阵形式：

或者说

。那么反过来说，

，也可以理解为原方程组的解的矩阵形式。在书中介绍了两种方法，分别是

雅可比法和 高斯-赛德尔法。

雅可比法

在数值线性代数中，雅可比法（Jacobi Method）是一种解对角元素几乎都是各行和各列的绝对值最大的值的线性方程组的算法。求解出每个对角元素并插入近似值。不断迭代直至收敛。这个算法是雅可比矩阵的精简版。方法的名字来源于德国数学家卡尔·雅可比。

给定一个n×n的线性方程组

。其中：

可以分解成对角部分

和剩余部分

:

线性方程组可以重写为：

雅可比法是一种迭代方法。在每一次迭代中，上一次算出的x被用在右侧，用来算出左侧的新的x。这个过程可以如下表示：

对每个元素可以用以下公式：

注意计算

需要

中除了自己之外的每个元素。我们不能立刻用

覆盖

，因为其他元素的计算也需要这些值。正是如此它比高斯-赛德尔法更加适合并行计算，更加适合GPU使用。

一般来说，雅可比法可以通过不断的迭代提升解的精度。但是并非绝对，有时候雅可比法反而会产生一个发散的迭代序列。保证收敛的条件是矩阵

为严格或不可约的

对角占优矩阵。严格的行对角占优矩阵指对于每一行，对角线上的元素之绝对值大于其余元素绝对值的和，即

有时即使不满足此条件，雅可比法仍可收敛。

高斯-赛德尔法

高斯-赛德尔法是数值线性代数中的一个迭代法，可用来求出线性方程组解的近似值。该方法以卡尔·弗里德里希·高斯和路德维希·赛德尔命名。

跟雅可比法类似，但不同的地方在于，求解过程中，所有已解出的k+1步元素都被直接使用，对于每个元素可以使用如下公式：

为了保证该方法可以进行，主对角线元素

需非零。

论文内容

首先，Jos提出速度的N-S方程，以及在该速度场下的密度方程：

首先，将求解密度方程，这允许读者熟悉求解器的不同组成部分，稍后会把求解密度方程的想法转移到相对更难的速度方程上。

论文中设计的网格。密度和速度指的均是格子的中心。网格包含额外的一圈以设定边界条件

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。整个密度求解器的步骤如下：

密度求解器的基本结构。每个时间步长都要经过三个阶段，求解密度方程的右侧三个项。

密度求解器

先描述密度场在不随时间变化的固定速度场中移动的求解器。密度方程式表明，单个时间步长上的密度变化由三个原因引起。这些原因是方程中等号右侧的三项：第一项表明密度变化应跟随速度场；第二项表明密度以一定速度扩散；第三项表明密度由于密度源的存在而增大。求解器将按相反顺序解析这些项，我们可以称之为添加密度源、扩散、移动。

密度源

第一项即密度源，很容易实现，一般会通过游戏引擎的某些功能填充，本例中由数组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]; 
}

扩散

通过扩散，格子与其邻居交换密度

第二项考虑到了可能的以

速率传播的扩散，当

时，密度变化最终将传播到整个格子体系内。本例假定格子仅与其直接相邻的4个格子进行密度交换。格子密度会因为流出而降低，同时也会因为邻居的流入而增加，总的来说，其净差为：

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()正是稍后讨论的边界格子。但不幸的是，如果针对一个较大的负数值（指

），密度有可能会产生震荡最后发散，从而模拟失败。因此，反过来，我们去寻找那些扩散后得到当前密度情况的密度情况。

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)所示。所以这是要维护两个网格数组的缘故：

代表密度的粒子在速度场中运动；(b)网格中心的粒子经过一个时间步长后到达下一位置；(c)反之，寻找一个时间步长后到达网格中心的粒子。

假定最开始的密度源包含在x0数组中，则一个时间步长内的变化为：

SWAP为交换数组指针的宏。

一个小小的实践

【1】中举了一个小小的例子作为两种迭代法的试验场。

假定有一根棍子，其左端被恒定为单位温度0，右端被恒定为单位温度1。最终棍子上的温度分布将趋于稳定，这可以被视作一维下的扩散。

将棍子均匀划分为8个点，每个点的温度被记为

，我们可以知道的是

，中间的六个值如何计算呢？很明显是线性插值，但是我们先提出一个方程来表示单个时间步长上的温度扩散情况：

是该点在时间步长前的温度。那么我们可以写出如下方程组：

写成矩阵就是：

可见这是一个规律且稀疏的矩阵，我们可以将其拆成更加简便的三个矩阵：

可以凝练出一个符号表达：

，即由对角矩阵和严格上、下三角矩阵组成。

雅可比迭代法

将变形过程表示为：

于是对于

的每个分量，我们可以写成：

高斯-赛德尔法

雅可比迭代法为什么适用于并行计算呢？区别在于是否使用

覆盖

。而高斯-赛德尔法选择覆盖，相当于雅可比法迭代一次才用最新的已知量修正一次；而后者每计算一行，都使用最新的已知量，相当于已经渐进地修正了。

我们假定在初始时棍子温度均为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]);
    }
}

6/7≈0.85714；白色是雅可比法；黄色是高斯-赛德尔法

精度达到第1、2、3、4位有效数字，高斯-赛德尔迭代用了8次、18次、37次、43次；而雅可比迭代用了16次、36次、74次、86次。

实际上，迭代的过程，就是一维情况下密度源（温度分布）扩散的过程。

速度求解器

之前的文章中，探讨了密度分布的变化，并提出了密度求解器的工作步骤。接下来，就轮到速度了。我们将速度方程解读为，速度随着时间步长的变化可归结于三个原因：力的增加、黏度扩散和自对流。自对流不太明显，可以简单的解释为速度场沿着自身移动方向移动。

正如同密度求解器那样，速度求解器亦可以如此：

在上述代码中，有一个名为project()的函数跟密度求解器不一致。该函数强迫速度保证质量守恒，这是被驱动的真实流体所具有的属性，如下图所示，也如同该系列第一篇专栏《【前置技能】流体力学的一些基础知识》所说的那样，质量守恒则带来无散场，无散场中拥有大量的涡流状图案，其视觉效果就是流体将具备丰富的涡流。

每个速度场都是一个不可压缩场（即无散场）和一个梯度场的和。

以下是将速度投影到一个无散场的project()的代码：

在流程中两次调用了project()，为了让advect()表现的更为精准。

边界条件

边界条件由set_bnd()函数处理，假定流体的边界是固体墙壁，意味着垂直于边界的速度水平分量为零。对于代码中考虑的密度和其他场，假定其是连续的：

以下展示了各个部分的代码是如何组合并运作的：

效果演示

Gif效果如下，想要学习请移步【6】：

总结

本文介绍了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