介紹常用的時間積分方法，及最終的求解過程。
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)