彈幕引擎的各種數值積分事情

tryourbreast

21/6: 最後補充了一小段量化誤差的方法




相信大家都幻想過做一個彈幕引擎。
然後用這個引擎做做彈幕遊戲，沒事做了可以模擬重力和天體，還有粒子交互，分子模擬，之類之類。

然而，沒有適切的數學方法，粒子的實際軌跡很快就會跟預定的軌跡偏離甚遠！

這個帖子就是說各種常微分方程的數值積分事情的。



作為出發點，可以舉一個最簡單的例子：模擬平地上的均勻重力。
大家都知道，在這情況下軌跡是拋物線。所以，普通人（初中物理程度）要模擬均勻重力時，就會寫出類似這樣的代碼：

1. x(t) = x(0) + v_x * t // 粒子現在的x值 = 初始x值 + 經過時間 * 初速的x方向的值
   
2. y(t) = y(0) + v_y * t + 1/2 * g * (t^2) // 粒子現在的y值 = 初始y值 + 經過時間 * 初速的y方向的值 + 1/2 * 加速度 * (經過時間)^2
   

复制代码

這種是不懂數學的人的彈幕引擎做法：程序員自己求出軌跡代表的參數方程；即粒子的位置x(t)是一個t的函數。

不過，實際上很多時候都是求不出軌跡的參數方程的。
例如我曾經跟別人討論過一個問題：一個沿正弦曲線（y=sin(x)）以均速行走的粒子，軌跡方程是甚麼？
經過一些運算就發現，這涉及了橢圓積分，而且是要求積分內函數的改變導致積分值為常數，並沒有解析解。

或者，如果你想處理三體重力問題，那麼基本上可以放棄找解析解了。

而這樣的引擎並沒有甚麼卵用，所以我們轉用了另一套思路：
乾脆就不求解析解了，每一步我們用這一步的參數（位置、初速、場強度之類的）算出下一步的值，一步一步走下去。
這思路下的均勻重力的代碼就會變成這樣：

1. // h是步距，兩步之間的時間差
   
2. v_x(t+h) = v_x(t) //  粒子下一步的x方向速度 = 粒子現在的x方向速度
   
3. x(t+h) = x(t) + v_x(t) * h // 粒子下一步的x值 = 粒子現在的x值 + 粒子現在的x方向速度 * 步距
   
4. 
   
5. v_y(t+h) = v_y(t) + g * h // 粒子下一步的y方向速度 = 粒子現在的y方向速度 + 加速度 * 步距
   
6. y(t+h) = y(t) + v_y(t) * h // 粒子下一步的y值 = 粒子現在的y值 + 粒子現在的y方向速度 * 步距
   

复制代码

一般來說，這是彈幕引擎最低必須要有的水平（東方官作的彈幕引擎就是以速度控制彈位置的）。

看到”最低必須水平”了嗎？

因為到了這裡，我們開始步入了一個十分重要的領域：常微分方程的數值解方法。
這當然也是一門藝術。





我們把粒子位置變成矢量w，把剛才的方程改寫並簡化一下：

1. // w是粒子位置矢量，v是粒子速度矢量
   
2. // w' = dw/dt
   
3. w' = v
   
4. v' = (0,g) // 視乎原方程會變成其他函數，即v' = a(t, w(t))

复制代码

現在很明顯能看到這裡有四個積分了，分別對應w和v的兩個方向的值。
進行任何數值運算前，我們都會先把手上的方程（n階常微分方程）拆成n個1階常微分方程，並將原方程進行主項變換，把最高級的導數階級移到左邊（其他式子都是導數鏈），最後由下面一步步進行積分，或者有時用矩陣並行積分。
（其實我們都知道實際存在的只有w，然而要一步步計算w就必須引入w的導數，所以只能把v也當成變量了。
而且在數字積分當中，w和w的導數的關係並非如此簡單。）

那麼，最後的積分便是

1. w(t+h) = w(t) + h * f1(h, a(t), w(t), v(t))
   
2. v(t+h) = v(t) + h * f2(h, a(t), w(t), v(t))
   
3. //f1和f2是由積分方法決定的函數

复制代码

回到算法上，剛才上面我們用的是標準的Euler method，算法是：

1. w(t+h) = w(t) + h * w'(t) // Euler method
   

复制代码

這是最原始的算法，準確度也是最低的，每步的誤差正比於h^2（用泰勒展開證明的），總誤差就正比於h^1了，所以又叫1階算法。
因此，所有稍為複雜的演算都不會用這種算法。一般來說，積分誤差在幾秒間就會有點明顯，除非把步距縮得很短。
彈幕是直線均速（或均加速）的話，倒是沒有問題。

然而，彈幕引擎限制我們只能選h = 1/60（每秒60幀），步距不由得我們調整，只能改進算法了。
（假設效率瓶頸位並不在位置運算上）
到了這裡，你有以下三個選擇：

a. 取中值，即 w'(t+h/2) 或 (w'(t) + w'(t+h) / 2
b. 保留多次結果，即同時考慮 w'(t), w'(t-h), w'(t-2h), w'(t-3h), ...
c. 意識到彈幕引擎是經典力學simulator

a選項會導致Runge Kutta系列，b選項會導致multistep系列，c選項會導致symplectic integrator（Euler-Cromer，Verlet和leapfrog）。



a. 取中值

常識告訴我們，如果選的h足夠小，w'(t)在t和t+h的平均值應該是個很好的估算，這樣可以把總誤差降到h^2：

1. w(t+h) = w(t) + h * (w'(t) + w'(t+h)) / 2     // Implicit Euler method
   
2. w(t+h) = w(t) + h * w'(t+h/2)                     // Explicit Euler method，又稱RK2，RK是Runge-Kutta的縮寫
   

复制代码

這樣看起來是沒問題的，除了一點：別忘了，w' = v，而v的計算中又有a，a是個t和w的函數；如果我們想算出w(t+h)就得知道w'(t+h)，然而如果我們想算出w'(t+h)，又得先知道w(t+h)！
這是一個死循環，而且如果我們知道了w(t+h)的值，我們根本就不需要算了（喂

不過這是當然的事，因為我們本來就只有w(t)的值，由一個點出發的積分法目前只有最簡單的Euler method。
也就是說，我們還是要從Euler method起步，有次序的一次次進行所有變數的積分：

1. // Implicit Euler method
   
2. v_0(t+h) = v(t) + h * a(h, w(t))
   
3. w_0(t+h) = w(t) + h * v(t)    // 以Euler method得到一對暫存矢量
   
4. 
   
5. v(t+h) = v(t) + h * ( a(h, w_0(t+h)) + a(h, w(t)) / 2      // w(t) = w_0(t)
   
6. w(t+h) = w(t) + h * (v_0'(t+h) + v'(t)) / 2    // v(t) = v_0(t)
   
7. 
   
8. // RK2
   
9. v_0(t+h/2) = v(t) + h/2 * a(h, w(t))
   
10. w_0(t+h/2) = w(t) + h/2 * w’(t)
    
11. 
    
12. v(t+h) = v(t) + h * a(h, w_0(t+h/2))
    
13. w(t+h) = w(t) + h * v_0(t+h/2)
    

复制代码

一般來說，如果我們想達到h^n的最終精度，就需要把加速度的函數a計算n次，和把v及w計算n次。
（這十分合理，付出多少收獲多少嘛）

RK2是RK4的簡化版，RK4可以達到h^4的精度（h^4在幾乎所有情況都很足夠了。RK4可以再進一步改成h^5精度的算法，然而通常並沒有必要），是個萬用算法，也是各種RK4變種的根基：

1. // RK4
   
2. v_0(t+h) = v(t) + h * a(h, v(t), w(t))   //假設a也是v的函數，方便說明
   
3. w_0(t+h) = w(t) + h * v(t)
   
4. 
   
5. v_1(t+h) = v(t) + h/2 * a(h/2, v(t) + v_0(t+h)/2, w(t) + w_0(t+h)/2)
   
6. w_1(t+h) = w(t) + h/2 * ( v(t) + v_0(t+h)/2 )
   
7. 
   
8. v_2(t+h) = v(t) + h * a(h/2, v(t) + v_1(t+h)/2, w(t) + w_1(t+h)/2)
   
9. w_2(t+h) = w(t) + h/2 * ( v(t) + v_1(t+h)/2 )
   
10. 
    
11. v_3(t+h) = v(t) + h * a(h/2, v(t) + v_2(t+h), w(t) + w_2(t+h))
    
12. w_3(t+h) = w(t) + h * ( v(t) + v_2(t+h) )
    
13. 
    
14. v(t+h) = v(t) + 1/6 * (v_0(t+h) + 2*v_1(t+h) + 2*v_2(t+h) + v_3(t+h) )
    
15. w(t+h) = w(t) + 1/6 * (w_0(t+h) + 2*w_1(t+h) + 2*w_2(t+h) + w_3(t+h) )
    

复制代码

看起來很複雜對吧？這的確能算是一個缺點，然而其實也不算甚麼缺點，因為可以先寫好公式實現時看著寫，捉急點的也可以用測試驗證嘛。
（我大學其中一個物理教授只教我們Implicit Euler method，因為他說RK4容易實現錯誤…（
另一個缺點是單次運算多，雖然可以讓步距放寬點，但是彈幕引擎就是得每秒60幀，不能多不能少，RK4的單次大量運算容易造成效能瓶頸，視乎情況需要改用RK2。

有一些RK4的變種是同步進行一個h^3精度和h^4精度的算法，比較兩個結果的差用以調整步距，做成adaptive step size（可變步距）的RK算法。
這個就更複雜了，我還沒看懂（喂

不過，如果你要做的是彈幕引擎，而不是拿各種常微分方程來玩的話，那麼RK2和RK4也許幫不了多少忙。
因為RK4模擬經典力學的一個最大問題是，能量是永遠不會守恆的，一定會逐漸減少。
物理模擬裡沒有人會用RK4，都是選的選項c。



b. 保留多次結果

RK4雖然是個泛用算法，然而問題很明顯：複雜，而且涉及大量運算。
這些問題都是由於我們只有w(t)一個值，所以舉步艱難，要試水好幾遍才敢跳下去。

一個很明顯的另類做法，就是拿著w(t), w(t+h), w(t+2h), ... 去算出下一個值：

1. // Adams–Bashforth methods
   
2. // 函數內的t是可以隨意變換的，把t往前或後移h的倍數就可以校齊已有數據
   
3. w(t+h) = w(t) + h * w'(t) // 1階，同Euler method
   
4. w(t+2h) = w(t+h) + h * ( 3/2 * w'(t+h) - 1/2 * w'(t) ) // 2階
   
5. w(t+3h) = w(t+2h) + h * ( 23/12 * w'(t+2h) - 4/3 * w'(t+h) + 5/12 * w'(t) ) // 3階
   
6. w(t+4h) = w(t+3h) + h * ( 55/24 * w'(t+3h) - 59/24 * w'(t+2h) + 37/24 * w'(t+h) - 3/8 * w'(t) ) // 4階
   
7. // 由此類推...
   

复制代码

這樣，我們只需要記住前面幾個結果，便能同樣輕易達到h^4的精度，以空間換取了時間。

不過我們可以做得更好。

1. // Adams–Moulton methods
   
2. // 下一步的值不再只存在於左邊，以此換取比階數高一級的精度
   
3. w(t) = w(t-h) + h * w'(t) // 0階
   
4. w(t+h) = w(t) + h * ( 1/2 * w'(t+h) + 1/2 * w'(t) ) // 1階，等價於梯形法則
   
5. w(t+2h) = w(t+h) + h * ( 5/12 * w'(t+2h) + 2/3 * w'(t+h) - 1/12 * w'(t) ) // 2階
   
6. w(t+3h) = w(t+2h) + h * ( 3/8 * w'(t+3h) + 19/24 * w'(t+2h) - 5/24 * w'(t+h) + 1/24 * w'(t) ) // 3階
   
7. w(t+4h) = w(t+3h) + h * ( 251/720 * w'(t+4h) + 646/720 * w'(t+3h) - 264/720 * w'(t+2h) + 106/720 * w'(t+h) - 19/720 * w'(t) ) // 4階
   
8. // 由此類推...

复制代码

我們可以把兩個方法並行採用，用兩次運算得出高一級的精度：

1. // 4階，精度h^5
   
2. // 先用Adams–Bashforth method算出近似值
   
3. v_0(t+4h) = v(t+3h) + h * ( 55/24 * v'(t+3h) - 59/24 * v'(t+2h) + 37/24 * v'(t+h) - 3/8 * w'(t) )
   
4. w_0(t+4h) = w(t+3h) + h * ( 55/24 * w'(t+3h) - 59/24 * w'(t+2h) + 37/24 * w'(t+h) - 3/8 * w'(t) )
   
5. // 再用Adams–Moulton method算出更準確的值
   
6. v(t+4h) = v(t+3h) + h * ( 251/720 * v_0'(t+4h) + 646/720 * v'(t+3h) - 264/720 * v'(t+2h) + 106/720 * v'(t+h) - 19/720 * v'(t) )
   
7. w(t+4h) = w(t+3h) + h * ( 251/720 * w_0'(t+4h) + 646/720 * w'(t+3h) - 264/720 * w'(t+2h) + 106/720 * w'(t+h) - 19/720 * w'(t) )
   

复制代码

最後用多少精度的算法視乎於實際需求。

不過，推導那些系數的方法有點複雜，嗯，有點複雜。
還是有需要時找別人算好的數字代入就好了。



c. 更聰明的做法

如果我們認同粒子是物理物件，遵循物理法則行走，那麼它們必須遵照牛頓的定律。
這件事我們剛才的算法並沒有考慮在內，而它會導致一些重要的特性，例如能量守恆和時間可逆性（正反方向可隨意行走），還有我們能求出更高階級的導數（儘管不在式子中）。

在牛頓定律的限制下開發積分的算法，會導致一系列的sympletic integrator。
（仔細來說，Hamiltonian formalism會導致跟p = mv和q（位置）有關的兩條一階偏微分方程，如果q與p無關，會變成常微分方程。）

1階的sympletic integrator會產生Euler-Cromer method（精度h^1）：

1. // 留意運行順序，由上至下
   
2. v(t+h) = v(t) + a(t, w(t))
   
3. w(t+h) = w(t) + a(t, v(t+h))
   
4. //或者
   
5. w(t+h) = w(t) + a(t, v(t))
   
6. v(t+h) = v(t) + a(t, w(t+h))
   

复制代码

2階的sympletic integrator會產生Verlet或leapfrog（leapfrog又稱Velocity verlet），精度h^2。
Verlet的特點是運算只需要位置w和加速度a，沒有了速度v，對於彈幕引擎來說也許有用：

1. // 留意運行順序，由上至下
   
2. // Verlet
   
3. w(t+h) = 2 * w(t) - w(t-h) + h^2 * a(t, w(t))
   
4. // Velocity Verlet (leapfrog)
   
5. w(t+h) = w(t) + h * v(t) + 1/2 * h^2 * a(t, w(t))    // s = u*t + (1/2)*a*t^2
   
6.     v(t+h) = v(t) + h * (a(t, w(t)) + a(t+h, w(t+h))) / 2
   
7.    

复制代码

更高階的推導一般都非常複雜，平常的演算2階就很足夠了。

這裡有兩個走向，一個是引入更高級數的導數：

1. // 留意運行順序，由上至下
   
2. // 精度h^4
   
3. // a'(t, w(t))有時可能算不出來，但普遍情況下轉移主項後做個求導就行。
   
4. w_0(t+h) = w(t) + h * v(t) + 1/2 * h^2 * a(t, w(t)) + 1/6 * h^3 * a'(t, w(t))
   
5. v_0(t+h) =  v(t) + h * a(t) + 1/2 * h^2 * a'(t, w(t))
   
6. v(t+h) = v(t) + h * (a(t, w(t)) + a_0(t+h, w(t+h))) / 2 + (a'(t, w(t)) - a_0'(t+h, w(t+h))) / 12
   
7. w(t+h) = w(t) + h * (v(t) + v(t+h)) / 2 + h^2 * (a(t, w(t)) - a_0(t+h, w(t+h))) / 12
   

复制代码

另一個是拿數個精心安排的leapfrog打出一套組合拳：

1. // Forest-Ruth algorithm
   
2. // 留意運行順序，由上至下
   
3. // 精度h^4
   
4. // define c0 = 1/(2-2^(1/3)), c1 = 1 - 2*c0
   
5. // 表達方式簡約化，每運算一次leapfrog就呼叫一次leapfrog函數
   
6. v0, w0 = leapfrog(t+h*c0, v(t), w(t), a(t, w(t)))
   
7. v1, w1 = leapfrog(t+h*c1, v0(t), w0(t), a(t, w0(t)))
   
8. v, w = leapfrog(t+h*c0, v1(t), w1(t), a(t, w1(t)))
   
9. // 這個方法可以輕易拓展到h^8或以上，精度h^2n的系數會是個長度為n陣列：（c0, c1, ... ,c_n），運行時由c0進行leapfrog到達c_n，再由c_(n-1)倒回去c0
   
10. // 加速度的運算次數比同階級RK少一次
    

复制代码

Sympletic integrator的優點有很多，除了剛才所述的能量守恆和時間可逆性外，對於震盪運動可以保持絕對的穩定性。
當然還有上述的，運算加速度次數比RK少一次的事。

不過在大步距下，sympletic integrator的表現會比RK4遜色，所以在遊戲和3DCG上一般用的RK4，一般標準推介的也是RK4。

然而．現在你知道各算法的優缺點，就可以自己決定怎麼處理你的子彈了。

另外，multistep method和RK系列都能改動以獲得sympletic integrator的特性，然而複雜度就更高了，除非你是做學術研究的，否則絕對沒有需要使用。



總結：
如果只是直線勻速彈的話，Euler method就足夠了。
更複雜的軌跡可以用2階算法，極端情況下（高加速度+高滯留時間）4階也極其足夠。

常微分方程的數值積分算法主要只有上述3種，比起積分的數值積分算法要簡單得多。

實現上，我的建議是先寫好Euler method，RK2，RK4，Verlet或leapfrog，和Forest-Ruth這麼幾個算法，然後跟據彈型和軌跡選擇一個合乎該彈需求的。
數值積分本來就有各種不同的算法，所以積分用的函數是應該想換就換的，覺得不妥的話換個函數名就行了。

如果你編程語言可以做並行運算，你可以寫出這樣的代碼：

1. // State-variable form
   
2. u(t+1) = u(t) + h * ( A * u )
   
3. //  u = (w,v)^T，是(w,v)的4維矢量的轉置
   
4. //  A =(
   
5.   //  {v_x(t),0,0,0},
   
6.   //  {0,v_y(t),0,0},
   
7.   //  {0,0,a_x(t, w(t)),0},
   
8. //  {0,0,0,a_y(t, w(t))}
   
9. //  )

复制代码

剩下來的好像也沒有甚麼好說的了，彈幕引擎遇到極端情況的機會微乎其微，除了算法實現出錯就沒甚麼了。

…對了！請務必記得你的彈幕引擎是補幀而不是掉幀的，因為不少算法遇到持續變化的步距時會出問題！





補充：
數值積分的錯誤是可以量化的：我們只能捕捉最高次導數的值，途中就會產生誤差。
經典力學一般只處理二階常微分方程，所以誤差大約等比於三次導數的值，即加速度的導數。
（證明的方法也是泰勒展開。）

不過，還有些系統本來就是混沌的，例如重力模擬，相近的初始條件最終一定會偏離甚遠。
我們會用Lyapunov time來表示混沌的程度：Lyapunov time被定義為，兩個無限相近的起始點，距離被放大e倍所需的時間。

太陽系模擬的Lyapunov time約是5千萬年（5千萬次公轉）。化學和水動力學的模擬時間大概是分鐘或秒鐘級。

步距當然也跟誤差成正比，步距足夠小的時候，兩者關係便由準確度的h^n表示，n為常數。

song_5007

感谢科普，慢慢吃

nialv7

如果是做物理引擎的话可能需要这样。弹幕真的需要这样的精确度吗。

不过真是学到了好多。

Capt.Murasa

看來加入不均衡引力之後，就應該不再是常微分方程了吧
要動用那些Laplace Transform﹑Bessel Functions...
再來可以有物質，暗物質，正負引力場，速度....也是位置的泛數...又可以玩那個泛函分析
三次導數.....又有啥可用....誒~~?

為個彈幕軌道...又激起我對數學的無限興趣

收藏~~~!!!

.void。

嗯。。。。。先收藏吧