介绍常用的时间积分方法，及最终的求解过程。
0 物理系统描述
在物理引擎中，借助牛顿第二运动定律对系统进行描述，即
有时候也会用 \(\boldsymbol{q}\) 来表示模型中节点的位置，那么系统描述即为：
上述方程就是物理引擎中的 ODE 部分。该方程组的解法主要有显式（Forward Euler、Semi-implicit Euler）、隐式（Backward Euler）等。
1、时间积分方法
在仿真计算过程中，已知模型中节点在 \(t\) 时刻的位置、速度等信息，进一步求解其在 \(t+1\) 时刻的速度、位置。
1.1 显式时间积分（Explicit / Forward Euler）
计算方式如下：
在该方法中，由模型中节点的位置 \(\boldsymbol{x}_{t}\) 直接计算得到受力 \(\boldsymbol{f}_{t}\) ，进而可直接计算得到 \(t+1\) 时刻模型中节点的速度、位置。
1.2 半隐式积分（Explicit / Semi-implicit Euler aka. Symplectic Euler）
计算方式如下：
在该方法中，同样由模型中节点的位置 \(\boldsymbol{x}_{t}\) 直接计算得到受力 \(\boldsymbol{f}_{t}\) ，进而可直接计算得到 \(t+1\) 时刻模型中节点的速度、位置。
Tips：Forward Euler 和 Semi-implicit Euler 略微不同，在计算 \(\boldsymbol{x}_{t+1}\) 的时候，一个是用了 \(\boldsymbol{v}_{t}\)，另一个是用了 \(\boldsymbol{v}_{t+1}\)。现在，Forward Euler 用的较少，Semi-implicit Euler 用的较多一些。虽然，Semi-implicit Euler 只差别了一点点，但是准确性上会有本质性的提升。
1.3 仿真流程（显式积分）
在使用显式（Forward Euler or Semi-implicit Euler）进行仿真的时候，仿真流程有如下几个步骤：
- 计算节点受力 \(\boldsymbol{f}_t = \boldsymbol{f}(\boldsymbol{x}_t)\)
- 计算新的速度 \(\boldsymbol{v}_{t+1} = \boldsymbol{v}_t + \Delta t\frac{\boldsymbol{f}_t}{m}\)
- 碰撞检测（此时会更正速度）
- 计算新的位置 \(\boldsymbol{x}_{t+1} = \boldsymbol{x}_t + \Delta t \boldsymbol{v}_{t+1}\) （Semi-implicit Euler）
显式时间积分器的性能缺陷： Easy to explore
关于稳定性 Stability 和爆炸 Explode 问题：
（略）
1.4 隐式积分（implicit Euler）
计算方式如下：
亦或者记作如下的形式
上述系统方程，可以转化成一个非线性的偏微分方程（PDE），通常有两种思路：（1）化简消去 \(\boldsymbol{x}_{t+1}\)；（2）化简消去 \(\boldsymbol{v}_{t+1}\)
（1）隐式积分求解 - 消去 \(\boldsymbol{x}_{t+1}\)
化简消去 \(\boldsymbol{x}_{t+1}\) 后，得到系统方程，即
这是一个关于 \(\boldsymbol{v}_{t+1}\) 的偏微分方程 PDE 。
（2）隐式积分求解 - 消去 \(\boldsymbol{v}_{t+1}\)
化简消去 \(\boldsymbol{v}_{t+1}\) 后，得到系统方程，即
这是一个关于 \(\boldsymbol{x}_{t+1}\) 的偏微分方程 PDE 。
Tips：在这里，消去 \(\boldsymbol{x}_{t+1}\) 而保留 \(\boldsymbol{v}_{t+1}\) 的用意，应该是便于进行碰撞处理。在一些物理引擎中，会先采用 \(\boldsymbol{x}_{t}\) 或 \(\boldsymbol{x}_{t} + \Delta t \boldsymbol{v}_t\) 作为模型中节点的位置，进行碰撞检测。得到碰撞结果后，进一步更新/限制节点的速度 \(\boldsymbol{v}_{t+1}\) ，这样的好处好象是稳定性会好一些，不会出现节点的位置穿过碰撞界限等。
2 物理引擎中 PDE 的求解
如第 1 节中看到，在物理仿真中，通过空间上的离散化，计算得到了模型中节点上的受力 \(\boldsymbol{f}\) ；通过运动方程，描述模型的运动规律，得到了一组 ODE；对模型的运动状态在时间上进行离散，并通过显式/隐式积分器，可以得到一组 PDE。那么，最终，就需要进行 PDE 的求解。（更为复杂的，比如带约束等情况，后面会再整理。这里只描述最基本的模型运动仿真）
（在物理引擎中，会得到各种 ODE、PDE、线性系统、非线性系统，名称上比较混乱。）其中涉及的求解方法也非常多，比如，线性化、牛顿法、Jacobin、CG（Conjugate Gradient）、优化隐式方法等等，非常多样
2.1 线性化（one step of Newton's method）
简单地求解，可以实现一步牛顿法，将 \(\boldsymbol{f}\) 在 \(\boldsymbol{x}_t\) 处一阶泰勒展开，得到如下：
（1）速度上的求解
操作之后，系统变成了线性系统。整理之后，为
（2）位置上的求解
这里不是很确定，还需要再对照资料确认一下。
如上所述，线性化之后会得到如下的线性系统：
对该线性系统的求解，可以采用 Jacobin / Gauss-Seidel iterations，或者 Conjugate gradients 等方法进行求解。
Jacobin iteration 求解线性系统的方法
单步的 Jacobin 迭代过程如下所示：
A = []
x = []
new_x = []
b = []
@ti.kernel
def iterate():
for i in range(n):
r = b[i]
for j in range(n):
if i != j:
r -= A[i,j] * x[j]
new_x[i] = r / A[i,i]
for i in range(n):
x[i] = new_x[i]
Jacobin 迭代的思想就是，每次只让一行（一个点）满足该方程，依次计算完成，就能让所有的点都（曾经）满足方程。多迭代几次，就能接近线性方程组的解了。（大概是这么个意思吧）
2.2 线性化（one step of Newton's method） - 进阶
在上节所述的线性系统中，加入一个系数 \(\beta\) ，那么，就会得到如下线性系统：
其中，当 \(\beta = 0\) 时，为 forward / semi-implicit Euler (explicit integrator)
当 \(\beta = 1/2\) 时，为 middle-point (implicit integrator)
当 \(\beta = 1\) 时，为 backward Euler (implicit integrator)