Conjugate Gradient

Posted on 2020-07-08

共轭梯度法（Conjugate Gradient）及其相关算法

1. 引言

一般来说，解$Ax=b$这种方程，我们的数值计算书[4]中讲得方法有：直接法（高斯消元，行、列主元，LU分解，LDLT分解），迭代法（Jacobi，Gauss-Seidel，松弛法）。分解法固然是最直接的方法，但是对于大型的稀疏矩阵，用分解法分解出来的矩阵，非零项很多，很吃内存。书中所讲的几种迭代法，迭代次数会比较多，而且迭代次数越多，迭代的效率越低。共轭梯度法是很好地解决了这一问题的算法。

2. 梯度下降法

首先考虑这样一个函数的极值问题，$f(\mathbf{x})=\frac{1}{2}\mathbf{x}^TA\mathbf{x}-\mathbf{x}^Tb+c$。与极值密切相关的是函数的导数，对$f(\mathbf{x})$关于$\mathbf{x}$求导，得到，$\frac{df(\mathbf{x})}{d\mathbf{x}}=\frac{1}{2}(A+A^T)\mathbf{x}-b$，当$A=A^T$时，即$A$为对称阵时，$\frac{df(\mathbf{x})}{d\mathbf{x}}=A\mathbf{x}-b$。这样，求解方程$A\mathbf{x}=b$的问题被转化为一个函数求极值的问题。而且形式上要求$A$为对称阵。也就是说共轭梯度法并不是对所有形式的矩阵A都适用，需要A满足对称正定的条件。

梯度下降法的思想是，给某个解$\mathbf{x_n}$，残差$\mathbf{r_n}=b-A\mathbf{x_n}=-f’(\mathbf{x_n})$。按照牛顿法的思想，让$\mathbf{x_n}$沿着负梯度$\mathbf{r_n}$的方向移动到$\mathbf{x_{n+1}}$。这里要求$\mathbf{x_{n+1}}$是$f(\mathbf{x})$在$\mathbf{x_n}+k\mathbf{r_n}$上的最小值。因此，有$\frac{df(\mathbf{x_n}+k\mathbf{r_n})}{dk}=\mathbf{r_n}^T\mathbf{r_{n+1}}=0$。所以，$\mathbf{r_n}^T(b-A(\mathbf{x_n}+k\mathbf{r_n}))=\mathbf{r_n}^T\mathbf{r_n}-k\mathbf{r_n}^TA\mathbf{r_n}=0$，所以，$k=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{r_n}^TA\mathbf{r_n}}$。这样，算法的流程就很清楚了。

$\alpha_n=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{r_n}^TA\mathbf{r_n}}$

$\mathbf{x_{n+1}}=\mathbf{x_n}+\alpha_n\mathbf{r_n}$

$\mathbf{r_{n+1}}=\mathbf{r_n}-\alpha_nA\mathbf{r_n}$

但是，存在一个问题，即便$\mathbf{r_j}^T\mathbf{r_n}=0,j<n$，$\mathbf{r_j}^T\mathbf{r_{n+1}}=\mathbf{r_j}^T\mathbf{r_n}-k\mathbf{r_j}^TA\mathbf{r_n}=0$也不一定满足，$\mathbf{r_{n+1}}$也不能保证与前面$n$个$\mathbf{r_i}$都垂直。即后面前进的方向会有部分与前面的某步重合，这使得梯度下降法的计算效率比较低。注意到，式子中有一项$\mathbf{r_j}^TA\mathbf{r_n}$干扰了正交性，所以需要在这上面进行一些调整。（这可能是共轭梯度法设计出来的早先思路）

3. 共轭梯度法

基于此，定义了一个$\mathbf{d_j}$与$\mathbf{d_k}$关于$A$正交，即$\mathbf{d_j}^TA\mathbf{d_k}=0$。

现在，我们尝试构建下CG的计算方法，从最简单的情形考虑起。设$\mathbf{r_0}=b-A\mathbf{x_0}$，$\mathbf{d_0}=\mathbf{r_0}$。此时，$\alpha_0=\frac{\mathbf{r_0}^T\mathbf{r_0}}{\mathbf{d_0}^TA\mathbf{d_0}}$，$\mathbf{x_1}=\mathbf{x_0}+\alpha_0\mathbf{d_0}$，$\mathbf{r_1}=\mathbf{r_0}-\alpha_0A\mathbf{d_0}$，此时，可以看出，$\mathbf{d_0}^T\mathbf{r_0}=\mathbf{r_0}^T\mathbf{r_0}$，$\mathbf{r_0}^T\mathbf{r_1}=\mathbf{r_0}^T\mathbf{r_0}-\alpha_0\mathbf{r_0}^TA\mathbf{d_0}=0$。可以知道$\mathbf{r_1}$与$\mathbf{r_0}$正交，因此由$\mathbf{r_1}$构建出一个新的更新方向就是(Schmit)，$\mathbf{d_1}=\mathbf{r_1}+\beta_0\mathbf{d_0}$，因为要满足关于$A$正交，所以，$\mathbf{d_0}^TA\mathbf{d_1}=\mathbf{d_0}^TA(\mathbf{r_1}+\beta_0\mathbf{d_0})=0$。此时，利用$A\mathbf{d_0}=\frac{\mathbf{r_0}-\mathbf{r_1}}{\alpha_0}$，得$\frac{\mathbf{r_0}^T-\mathbf{r_1}^T}{\alpha_0}\mathbf{r_1}=-\beta_0\mathbf{d_0}^TA\mathbf{d_0}$，所以，$\beta_0=\frac{\mathbf{r_1}^T\mathbf{r_1}}{\alpha_0\mathbf{d_0}^TA\mathbf{d_0}}=\frac{\mathbf{r_1}^T\mathbf{r_1}}{\mathbf{r_0}^T\mathbf{r_0}}$，$\mathbf{d_1}=\mathbf{r_1}+\frac{\mathbf{r_1}^T\mathbf{r_1}}{\mathbf{r_0}^T\mathbf{r_0}}\mathbf{d_0}$。此时，注意到，$\mathbf{r_0}^TA\mathbf{d_1}=\mathbf{d_0}^TA\mathbf{d_1}=0$，此式备用。

下面来到第$n$步，根据简单情形，我们利用归纳法构建行进的方向。首先，假设，当$i,j\leqslant n \quad (n\geqslant1)$时，有这些关系成立。$\mathbf{d_i}^T\mathbf{r_i}=\mathbf{r_i}^T\mathbf{r_i}$，$\mathbf{r_i}^T\mathbf{r_j}=0\quad(i\not=j)$，$\mathbf{d_i}^TA\mathbf{d_j}=0\quad(i\not=j)$，$\mathbf{r_i}^TA\mathbf{d_j}=0\quad(i<j)$，$\mathbf{d_i}=\mathbf{r_i}+\beta_{i-1}\mathbf{d_{i-1}}$。

现在来到第$n+1$步，首先可以确定，$\alpha_n=\frac{\mathbf{d_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}$，这样就得到，$\mathbf{r_{n+1}}=\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}A\mathbf{d_n}$。此时，证明$\mathbf{r_{n+1}}$与$\mathbf{r_i}(i\leqslant n)$的正交性。$\mathbf{r_n}^T\mathbf{r_{n+1}}=\mathbf{r_n}^T\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}\mathbf{r_n}^TA\mathbf{d_n}=\mathbf{r_n}^T\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{(\mathbf{r_n}^T+\beta_{n-1}\mathbf{d_{n-1}}^T)A\mathbf{d_n}}\mathbf{r_n}^TA\mathbf{d_n}=\mathbf{r_n}^T\mathbf{r_n}-\mathbf{r_n}^T\mathbf{r_n}=0$，而$\mathbf{r_i}^T\mathbf{r_{n+1}}=\mathbf{r_i}^T\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}\mathbf{r_i}^TA\mathbf{d_n}=0+0,(i<n)$，有了$\mathbf{r_{n+1}}$，就可以构造$\mathbf{d_{n+1}}$，设$\mathbf{d_{n+1}}=\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}$，则$\mathbf{d_n}^TA\mathbf{d_{n+1}}=\mathbf{d_n}^TA(\mathbf{r_{n+1}}+\beta_n\mathbf{d_n})=0$，所以利用$A\mathbf{d_n}=\frac{\mathbf{r_n}-\mathbf{r_{n+1}}}{\alpha_n}$，$\beta_n=\frac{\mathbf{d_n}^TA\mathbf{r_{n+1}}}{\mathbf{d_n}^TA\mathbf{d_n}}=\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}$，$\mathbf{d_{n+1}}=\mathbf{r_{n+1}}+\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}\mathbf{d_n}$，现在，证明$\mathbf{d_{n+1}}$关于$A$的正交性。$\mathbf{d_j}^TA\mathbf{d_{n+1}}=\mathbf{d_j}^TA\mathbf{r_{n+1}}+\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}\mathbf{d_j}^TA\mathbf{d_n}=\frac{\mathbf{r_j}^T-\mathbf{r_{j+1}}^T}{\alpha_j}\mathbf{r_{n+1}}+\beta_n\mathbf{d_j}^TA\mathbf{d_n}=0(j<n)$，对于$j=n$的情形，$\mathbf{d_n}^TA\mathbf{d_{n+1}}=\frac{\mathbf{r_n}^T-\mathbf{r_{n+1}}^T}{\alpha_n}\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}^TA\mathbf{d_n}=-\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\alpha_n}+\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}\mathbf{d_n}^TA\mathbf{d_n}=0$。然后，还有一些辅助的式子的证明，$\mathbf{d_{n+1}}^T\mathbf{r_{n+1}}=\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}^T\mathbf{r_{n+1}}=\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}+\beta_n\mathbf{d_n}^T(\mathbf{r_n}-\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}A\mathbf{d_n})=\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}$。$\mathbf{r_i}^TA\mathbf{d_{n+1}}=(\sum_{j}^{i-1}a_j\mathbf{d_j})A\mathbf{d_{n+1}}=0$。至此，我们证明了所有的结论在第$n+1$步后也成立。这就是CG的形式的由来，关键在于两点，一个是关于$A$正交，另一个是 $\mathbf{d_i}=\mathbf{r_i}+\beta_{i-1}\mathbf{d_{i-1}}$的形式的确立。

$\alpha_n=\frac{\mathbf{r_n}^T\mathbf{r_n}}{\mathbf{d_n}^TA\mathbf{d_n}}$

$\mathbf{x_{n+1}}=\mathbf{x_n}+\alpha_n\mathbf{d_n}$

$\mathbf{r_{n+1}}=\mathbf{r_n}-\alpha_nA\mathbf{d_n}$

$\beta_n=\frac{\mathbf{r_{n+1}}^T\mathbf{r_{n+1}}}{\mathbf{r_n}^T\mathbf{r_n}}$

$\mathbf{d_{n+1}}=\mathbf{r_{n+1}}+\beta _n\mathbf{d_n}$

4. 预条件共轭梯度法(Preconditioner Conjugate Gradient)

按照算法写了代码之后，发现迭代的速度并不是非常快。这与谱半径有关，即最大的特征值绝对值。所以，又有PCG(Preconditioner conjugate gradient)来加速收敛。

PCG的思想是，$A$是对称正定矩阵，所以，可以将$A$写成$A=LL^T$的形式。但是$L$有很多非零项，所以可以设计出某些Preconditioner来近似$A$，这样使得$L^{-1}AL^{-T}\approx I$，而求解$I\mathbf{x}=\mathbf{b}$是非常容易的，所以，这样就使收敛变得更容易了。

这里先推导下引入preconditioner后的CG流程形式是怎么来的。

$A\mathbf{x}=b\Leftrightarrow E^{-1}AE^{-T}(E^T\mathbf{x})=E^{-1}\mathbf{b}\Leftrightarrow \hat{A}\hat{\mathbf{x}}=\hat{b}$

，令$M=EE^T$，CG的流程变为：

$\hat{\mathbf{r}}_n=\hat{\mathbf{b}}-\hat{A}\hat{\mathbf{x}}_n=E^{-1}\mathbf{r}_n$

$\mathbf{\hat{d}_n}=E^{-1}\mathbf{d_n}$

$\hat{\alpha}_n=\frac{\hat{\mathbf{r}}_n^T\hat{\mathbf{r}}_n}{\hat{\mathbf{d}}_n^T\hat{A}\hat{\mathbf{d}}_n}=\frac{\mathbf{r}_n^TE^{-T}E^T\mathbf{r}_n}{\mathbf{d}_n^TE^{-T}E^{-1}AE^{-T}E^{-1}\mathbf{d}_n}$

$\hat{\mathbf{x}}_{n+1}=\hat{\mathbf{x}}_n+\hat{\alpha}_n\hat{\mathbf{d}}_n\Rightarrow \mathbf{x}_{n+1}=\mathbf{x}_n+\hat{\alpha}_nE^{-T}E^{-1}\mathbf{d}_n$

$\hat{\beta}_n=\frac{\hat{\mathbf{r}}_{n+1}^T\hat{\mathbf{r}}_{n+1}}{\hat{\mathbf{r}}_n^T\hat{\mathbf{r}}_n}=\frac{\mathbf{r}_{n+1}^TE^{-T}E^{-1}\mathbf{r}_{n+1}}{\mathbf{r}_n^TE^{-T}E^{-1}\mathbf{r}_n}$

$\hat{\mathbf{d}}_{n+1}=\hat{\mathbf{r}}_{n+1}+\hat{\beta}_n\hat{\mathbf{d}}_n\Rightarrow E^{-T}E^{-1}\mathbf{d}_{n+1}=E^{-T}E^{-1}\mathbf{r}_n+\hat{\beta}_nE^{-T}E^{-1}\mathbf{d}_n$

令$E^{-T}E^{-1}\mathbf{r_n}=M^{-1}\mathbf{r_n}=\mathbf{Z_n}$，$\mathbf{d’_n}=E^{-T}E^{-1}\mathbf{d_n}$，得到，

$\mathbf{Z_n}=M^{-1}\mathbf{r_n}$

$\hat{\alpha}_n=\frac{\mathbf{r}_n^T\mathbf{Z_n}}{\mathbf{d’}_n^TA\mathbf{d’}_n}$

$\mathbf{x}_{n+1}=\mathbf{x}_n+\hat{\alpha}_n\mathbf{Z_{n+1}}$

$\hat{\beta}_n=\frac{\mathbf{r}_{n+1}^T\mathbf{Z_{n+1}}}{\mathbf{r}_n^T\mathbf{Z_n}}$

$\mathbf{d’_{n+1}}=\mathbf{Z_n}+\hat{\beta}_n\mathbf{d’_n}$

常用的Preconditioner，

Jacobi preconditioner：

$A=F+D+F^T$，$M=D$

Incomplete cholesky preconditioner[3]：

$A=F+D+F^T$，$L=FE^{-1}+E$，$LL^T=(FE^{-1}+E)(E^{-1}F^T+E)=FE^{-1}E^{-1}F^T+F^T+F+EE$，为了确定$E$的数值，使$FE^{-1}E^{-1}F^T+F^T+F$的对角项等于$D$。注意$LL^T$并不是完全等于$A$的，有应该是0项的成了非0项。

SSOR preconditioner：

$A=F+D+F^T$，$M=(D+L)D^{-1}(D+L)^T$，

带上参数之后，$M=\frac{1}{2-\omega}(\frac{1}{\omega}D+L)(\frac{1}{\omega}D)^{-1}(\frac{1}{\omega}D+L)^T$。

Preconditioner的关键是需要$\mathbf{Z_n}=M^{-1}\mathbf{r_n}$这一步算起来很容易。上面的几个构建的都是可以直接求解（LDLT，LU，追赶法）的矩阵。

从实践的结果看，加preconditioner的确要快。Incomplete cholesky preconditioner要快于SSOR preconditioner。同一个问题，SSOR要132步收敛，ICP要52步收敛。

5. 总结

接下来进行MultiGrid以及FFT的学习，这两个都可以快速地求解线性方程组而且广泛的应用在很多地方，属于高级技巧，必须掌握。印象中，CUDA里有一个FFT解NS方程的example，可以参考一下。

6. 参考资料

[1] An Introduction to the Conjugate Gradient Method Without the Agonizing Pain

[2] 资料[1]的中文翻译

[3] Fluid simulation for computer graphics

[4] 数值计算方法与算法（第二版）

[5] Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods

Position-based Dynamics

Posted on 2020-06-23

Position-based Dynamics

参考资料

[1] http://graphics.stanford.edu/courses/cs348c/PA2_PBD2016/index.html

记录

Posted on 2020-06-15 Edited on 2020-07-08

记录

1.路由器与手机热点或手机USB共享网络的区别

问题描述：
连接TeamViewer时，当电脑用家里的路由器连接wifi时，无法连上远程电脑。
此时如果手机打开热点，或者手机连上家里的wifi然后将手机连上电脑调成USB共享网络。就能连上远端电脑。
在使用git clone时也出现了类似的状况，用git clone无法下载，出现了错误“port：443“。即使将 .gitconfig 中proxy的IP地址改成代理端口，还是需要手机开热点的情况才能下载。但是，手机用USB共享网络就不行了。
手机USB与手机热点的区别究竟是什么？

手机端口的问题吗

2.类的指针作为参数传递进函数，无法访问到类中的数据

3.Python中类作为函数的参数，作为类元素的测试

问题描述：
Python中的函数能不能用类做参数，类中能不能用其他类做类成员。
做这样一个测试：

class A:
    def __init__(self,a):
        self.x=a   
    def add(self):
        self.x=self.x+1
        
class B:
    def __init__(self,a):
        self.at=a       
    def add(self):
        self.at.x+=2
        
def prints(a):
    print(a.x)     

a=A(2)
prints(a)
b=B(a)

b.add()
prints(a)
prints(b.at)

a.add()
prints(a)
prints(b.at)

输出的结果为：

从测试结果可以知道，类作为参数传递时，相当于是引用。作为类成员时，也是引用。因此面向对象编程会很方便。

4.Gibbs 现象

问题描述：
做傅里叶变换时，非连续的函数最终得到的是理论上是连续的傅里叶函数，但是计算机一离散就会出问题，会有小波纹。
查资料，写code

5.隐式与显式的区别，条件数与收敛速度

问题描述：
一直没有真正理解显式与隐式的区别，需要好好整理理解下
做一个demo学习，一维信号传递，二维图像，旋转加平移。

6.珠链喷泉（Chain fountain）模拟

问题描述：
看能不能用简单的Mass-Spring模拟出来

模拟方案如下：

数据结构：

将弹簧质点存在链表里，使用链表的原因是为了方便对整条链进行操作。

算法

S1. 遍历弹簧质点，计算每个质点受到其他质点作用的力$ f^{spring}_i $，同时给每个质点添加重力$ f^{gravity}_i $。

S2. 对所有质点，计算其距离杯壁的距离$d$，当距离小于厚度$t$时，对质点施加一个指向法向的弹性力。

S3. 更新质点的坐标。进行碰撞处理。对于一直落在地面的粒子，进行重复利用的操作。

具体说明：

S1：没有什么具体要注意的。

S2：这一步涉及到质点与多边形的交互问题。处理的方法是，根据多边形的距离场，定义一个排斥项$F=k(d-t)$。这一排斥力项，一方面可以满足能量守恒、动量守恒。另一方面，可以避免很多复杂的处理。

之前一段时间，我在尝试做PolyBridge桥面要考虑桥面与小球的交互时，想通过多边形引入冲量，然后计算各种速度等，但是从实践的结果来看，这种处理方法很不好，小球根本无法稳定的停在桥面上。

从结果来看，这样简单的处理可以获得很好的结果。所以，以后遇到物体与多边形交互的情形，可以用类似的方法处理。

顺便一提的是，过程中也遇到一个问题，如果链的劲度系数不够的话，会出现链条从多边形边缘穿过的现象。这也符合物理实际，真实的链劲度系数应该非常地大。

这种处理方式应该叫penalty-based system[1]，这种处理方式的问题在于力太强不稳定，力太弱会出现penetrate的现象。

还有一种方式是projection collision reactions[1]，把粒子从碰撞体中拔出来。

S3：这一步有一个注意点，Verlet Integration，在分子动力学模拟中非常常用，Verlet Integration可以通过设置粒子位置来完成碰撞，而带有速度项的更新会引入一些碰撞处理的麻烦，但是，这样一来，物理的真实性就不一定可行了。这里记录它的基本思想。

$x(t+\Delta t)=x(t)+\frac {dx}{dt}\Delta t+\frac {1}{2}\frac {d^2x}{dt^2}{\Delta t}^2+\frac {1}{6}\frac {d^3x}{dt^3}{\Delta t}^3+O({\Delta t}^4)$ $x(t-\Delta t)=x(t)-\frac {dx}{dt}\Delta t+\frac {1}{2}\frac {d^2x}{dt^2}{\Delta t}^2-\frac {1}{6}\frac {d^3x}{dt^3}{\Delta t}^3+O({\Delta t}^4)$ $x(t+\Delta t)+x(t-\Delta t)=2x(t)+\frac {d^2x}{dt^2}{\Delta t}^2+O({\Delta t}^4)$

从式子中可以看到，Verlet积分是local四阶精度，global二阶精度。

除此以外，还有velocity 的Verlet积分，它的形式是：

$v(t+\frac {1}{2}\Delta t)=v(t)+\frac {1}{2}a(t)\Delta t$ $x(t+\Delta t)=x(t)+v(t+\frac {1}{2}\Delta t)\Delta t$ $v(t+\Delta t)=v(t+\frac {1}{2}\Delta t)+\frac {1}{2}a(t+\Delta t)\Delta t$

这种更新形式与Leapfrog方法相似。值得注意的是，时间精度提高一倍，但是加速度的计算开销没有发生变化。

这种方法有两个优点：

一个是time-reversibility，时间回溯性，即向前n步，掉头更新n步可以回头。这一点从

另一个是symplectic nature，辛结构？对于满足辛变换的系统会使用辛算法进行离散求解，辛积分。

这里埋个坑，为什么Euler method是一阶精度，Backward Eular精度是多少，留到问题5中进行解决。

为了让模拟的时间变长，当质点撞到地面的时候，将落在地面上的粒子移到杯子中。连在尾部。这里就体现了链表的优势，操作非常简单。后来想了下，如果用数组做索引也可以。

最终的模拟结果如图：

可以看到，出现了珠链喷泉的现象。从一定程度上反应了，这个现象来源于杯底的作用力。

单位向量的求导

在隐式求解Mass-Spring时，需要对单位向量进行求导。求导的公式如下：

$\frac {\partial {\hat x}}{\partial {x}}= \frac {I-{\hat x}{\hat x}^T}{\left |{x}\right |}$

[1] https://en.wikipedia.org/wiki/Verlet_integration

[2] https://en.wikipedia.org/wiki/Leapfrog_integration

7.Hexo图片、公式插入的问题

问题描述：如何在Hexo中插入图片、公式

不得不说，折腾了好一会才折腾好。

插入图片

图片插入是需要安装一个插件叫hexo-asset-image。

1	npm install hexo-asset-image --save

安装以后，一开始还是无法插入图片。观察以后发现有几个问题。一个是输出信息的时候，“Update linke—>”后面的地址有各种%2E之类的URL编码。这种情况应该是因为文章的名字中有中文（但是Hexo不是支持中文的吗…）。地址中图片名这一块有一个%5C，这个是反斜杠的编码。

之后找了很多教程基本雷同，没能解决问题，但是有一个教程很有意思。他遇到了教程里没有提到的问题，他最终是看hexo-asset-image中index.js的代码解决的。所以模仿了教程的阅读思路，看了下代码，发现应该是

1	var src = $(this).attr('src').replace('\\', '/');

这边有点问题，这个替换没替换成功。所以，看了下md中的图片地址，都是反斜杠连接的，问题应该在这，所以改成了斜杠。问题就没了。

斜杠与反斜杠：反斜杠是win的文件路径。斜杠是unix的文件路径。网络文件的地址都是统一的斜杠。C里反斜杠还有转义的含义。

后来在解决公式问题的时候似乎改了一些东西。之前的问题好像自动消失了。但是可以确定的是之前是斜杠与反斜杠转换这块有问题。

插入公式

直接参考教程即可，很清楚。

https://www.youtube.com/watch?v=hdn6mqubh98&list=PLgPpaTsP_3Dq5KsWd6-8wmjNs3ipnvCU3

格子玻尔兹曼方法(LBM)学习记录

Posted on 2020-06-15

格子玻尔兹曼方法(LBM)学习记录

一、前言

前几天看到一个用LBM模拟卡门涡街（Karman Vortex Street）的demo，被深深地吸引。因为LBM也是目前工业领域模拟比如：多孔介质等比较常用的模拟方法，同时，LBM的理论背景和MPC有一定的渊源，所以决定找些资料学习下LBM背后的原理。

Yiming Xia's Blog

Conjugate Gradient

共轭梯度法（Conjugate Gradient）及其相关算法

1. 引言

2. 梯度下降法

3. 共轭梯度法

4. 预条件共轭梯度法(Preconditioner Conjugate Gradient)

5. 总结

6. 参考资料

Position-based Dynamics

Position-based Dynamics

参考资料

记录

记录

1.路由器与手机热点或手机USB共享网络的区别

2.类的指针作为参数传递进函数，无法访问到类中的数据

3.Python中类作为函数的参数，作为类元素的测试

4.Gibbs 现象

5.隐式与显式的区别，条件数与收敛速度

6.珠链喷泉（Chain fountain）模拟

7.Hexo图片、公式插入的问题

格子玻尔兹曼方法(LBM)学习记录

格子玻尔兹曼方法(LBM)学习记录

一、前言

二、理论模型