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