差分法と数値積分法
速度や加速度のような「変化」を、飛び飛びの値しか扱えないコンピュータで計算するための、差分法と数値積分法について説明します。
<注意>
このページでは、数式を表示するためにMathJaxを使用しています。ブラウザの設定でスクリプト実行をブロックしている場合、ページ内の数式が正しく表示されません。
微分と積分(連続な関数の場合)
はじめに、高校数学で学ぶ一変数の連続な関数の微分と積分を、簡単に振り返ります。
厳密な定義には立ち入りません。
微分
関数\(\ y=f(x)\ \)で、\(\ x\ \)が\(\ h\ \)だけ増えたときの\(\ y\ \)の増え方の割合\(\ \cfrac{f(x+h)-f(x)}{h}\ \)を、平均変化率と言います。
\(\ h\ \)を限りなく0に近づけたときの値が、その点での変化率です。
定義 (導関数)
\[ f'(x)=\lim_{h\to0}\frac{f(x+h)-f(x)}{h} \]
グラフで見ると、\(\ f'(x)\ \)は点\(\ (x,\ f(x))\ \)での接線の傾きです。
たとえば\(\ (x^n)'=nx^{n-1}\ \)で、\(\ f(x)=x^2\ \)なら\(\ f'(x)=2x\ \)です。
\(\ h\ \)を小さくしていくと、2点を通る直線(割線)は、点\(\ x\ \)での接線に近づいていきます。
積分
関数\(\ f(x)\ \)のグラフと\(\ x\ \)軸の間の、\(\ x=a\ \)から\(\ x=b\ \)までの面積を、定積分と言います。
この範囲を幅\(\ \Delta x\ \)の細い長方形に分けて面積を足し、幅を限りなく0に近づけたものと考えられます。
長方形の数を増やして幅を狭くするほど、長方形の面積の和は、グラフの下の面積に近づきます。
定義 (定積分)
\[ \int_a^b f(x)\,dx=\lim_{n\to\infty}\sum_{k=0}^{n-1}f(x_k)\,\Delta x\qquad\left(\Delta x=\frac{b-a}{n},\ x_k=a+k\Delta x\right) \]
微分と積分は、互いに逆の操作です。
定理 (微分積分学の基本定理)
\(\ F(x)=\displaystyle\int_a^x f(t)\,dt\ \)とすると、\(\ F'(x)=f(x)\ \)である。
たとえば\(\ \displaystyle\int x^n\,dx=\frac{x^{n+1}}{n+1}+C\ \)(\(\ C\ \)は定数)で、微分すると\(\ x^n\ \)に戻ります。
位置・速度・加速度
時刻\(\ t\ \)での位置を\(\ x(t)\ \)とすると、位置を微分したものが速度、速度を微分したものが加速度です。
\[ v(t)=\frac{dx}{dt},\qquad a(t)=\frac{dv}{dt} \]逆に、加速度を積分すると速度が、速度を積分すると位置が求まります。
加速度が一定値\(\ a\ \)のとき(等加速度運動)、初めの位置を\(\ x_0\ \)、初めの速度を\(\ v_0\ \)とすると、次のようになります。
ジャンプで考えてみます。上向きを正、重力加速度の大きさを\(\ g\ \)とすると、高さは\(\ y(t)=v_0t-\cfrac{1}{2}gt^2\ \)です。
速度\(\ v_0-gt\ \)が0になる\(\ t=\cfrac{v_0}{g}\ \)で最高点に達し、その高さは\(\ \cfrac{v_0^2}{2g}\ \)です。
\(\ v_0=4,\ g=0.25\ \)の場合のグラフです。
高さのグラフの接線の傾きが速度、速度のグラフの傾きが加速度です。
速度が0になる時刻に、高さが最高点になっています。
このページでは、数式を分かりやすくするために上向きを正にしています。
SNESの画面のY座標は下向きに増えるので、プログラムでは符号が逆になります。
飛び飛びの値しか扱えない
ゲームのプログラムは、キャラクターの位置や速度を、1フレームに1回だけ計算します。
値があるのは0フレーム目、1フレーム目、2フレーム目、…の時点だけで、その間の値はありません。
プログラムが持っているのは点の値だけで、破線の曲線そのものはどこにもありません。
微分も積分も「幅を限りなく0に近づける」ことで定義されていましたが、飛び飛びの値ではそれができません。
幅は、どんなに小さくしても1フレームです。
音声も同じで、値があるのはサンプリングした時点だけです(サンプリングとデジタルフィルタ)。
そこで、コンピュータでは、微分の代わりに差(差分)を、積分の代わりに和を使います。
これが差分法と数値積分法です。
以下では、\(\ n\ \)フレーム目の値を\(\ x_n\ \)のように書きます。
時間の刻み幅が1フレームなので、「幅で割る」計算は1で割るだけになり、式から消えます。
| 連続な関数 | 飛び飛びの値 |
|---|---|
| 微分\(\ \cfrac{dx}{dt}\ \) | 差分\(\ x_{n+1}-x_n\ \) |
| 積分\(\ \displaystyle\int v\,dt\ \) | 和\(\ \displaystyle\sum v_k\ \) |
| 微分と積分は逆の操作 | 差分と和も逆の操作 |
差分と和が逆の操作であることは、差分を足していくと、途中が打ち消し合って元の値に戻ることから分かります。
\[ x_0+\sum_{k=0}^{n-1}(x_{k+1}-x_k)=x_n \]差分法
3種類の差分
\(\ n\ \)フレーム目の変化率の近似として、次の3つがよく使われます。
| 名前 | 式 | 使う値 |
|---|---|---|
| 前進差分 | \(\ x_{n+1}-x_n\ \) | 今と次 |
| 後退差分 | \(\ x_n-x_{n-1}\ \) | 前と今 |
| 中心差分 | \(\ \cfrac{x_{n+1}-x_{n-1}}{2}\ \) | 前と次 |
なめらかに変わる値では、中心差分がいちばん本当の変化率に近くなります。
たとえば\(\ x_n=n^2\ \)の場合、本当の変化率は\(\ 2n\ \)です。
前進差分は\(\ 2n+1\ \)、後退差分は\(\ 2n-1\ \)とずれますが、中心差分はちょうど\(\ 2n\ \)になります。
曲がっている曲線では、前進差分と後退差分の直線は接線から反対向きにずれ、中心差分の直線は接線とほぼ重なります。
一般に、前進差分と後退差分の誤差は刻み幅に比例し、中心差分の誤差は刻み幅の2乗に比例します。
位置の記録から速度と加速度を求める
差分をもう一度取ると(2階差分)、加速度の近似になります。
\[ (x_{n+1}-x_n)-(x_n-x_{n-1})=x_{n+1}-2x_n+x_{n-1} \]次の表は、ジャンプしたキャラクターの高さを毎フレーム記録し、差分を取ったものです。
| フレーム | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| 高さ | 0 | 3.75 | 7.25 | 10.5 | 13.5 | 16.25 |
| 1階差分(速度) | 3.75 | 3.5 | 3.25 | 3.0 | 2.75 | |
| 2階差分(加速度) | −0.25 | −0.25 | −0.25 | −0.25 |
2階差分が一定なので、等加速度運動で、重力は1フレームに0.25ピクセル/フレームだと分かります。
1階差分を前にたどると、ジャンプの初速度(このあと説明する計算方法では4ピクセル/フレーム)も分かります。
既存のゲームの動きを調べるときにも、この方法が使えます。
座標を毎フレーム記録して差分を取れば、プログラムを読まなくても、初速度や重力、最大の落下速度などの値を推定できます。
- 2階差分が0なら、速度が一定(等速運動)。
- 2階差分が一定なら、加速度が一定(重力など)。
- 2階差分が途中で変わったら、そのフレームで加速度が変わった(ボタンを離した、最大の落下速度に達した、など)。
数値積分法
数値積分法は、加速度から速度を、速度から位置を、毎フレーム足していくことで求める方法です。
これは、速度のグラフの下の面積を、幅1フレームの長方形などで近似していることにあたります。
ここでは、重力\(\ g\ \)を受けて動く物体の、\(\ n\ \)フレーム目の高さを\(\ y_n\ \)、速度を\(\ v_n\ \)とします。
速度の更新\(\ v_{n+1}=v_n-g\ \)はどの方法でも同じで、違うのは高さに足す速度です。
| 方法 | 高さの更新 | 面積の近似 |
|---|---|---|
| オイラー法 | \(\ y_{n+1}=y_n+v_n\ \) | 区間の始めの速度の長方形 |
| 半陰的オイラー法 | \(\ y_{n+1}=y_n+v_{n+1}\ \) | 区間の終わりの速度の長方形 |
| 台形法 | \(\ y_{n+1}=y_n+\cfrac{v_n+v_{n+1}}{2}\ \) | 台形 |
速度が減っていくジャンプの上昇中では、オイラー法の長方形は本当の面積より大きく、半陰的オイラー法の長方形は小さくなります。
台形法の台形は、速度のグラフが直線なら本当の面積と一致します。
プログラムでは、オイラー法と半陰的オイラー法の違いは、2つの処理の順番だけです。
- オイラー法: 速度を位置に足してから、重力を速度に足す。
- 半陰的オイラー法: 重力を速度に足してから、速度を位置に足す。
半陰的オイラー法は、シンプレクティック・オイラー法とも呼ばれます。
ゲームでよく使われるのはこの方法です。
ジャンプの高さで比べる
初速度\(\ v_0=4\ \)ピクセル/フレーム、重力\(\ g=0.25\ \)ピクセル/フレーム²でジャンプさせます。
8.8形式の固定小数点数では、初速度が$0400、重力が$0040です(固定小数点で動かすと同じ重力です)。
連続な式では、16フレーム目に最高点32ピクセルに達し、32フレーム目に着地します。
| 方法 | 最高点 | 着地するフレーム |
|---|---|---|
| 連続な式 | 32 | 32 |
| オイラー法 | 34 | 33 |
| 半陰的オイラー法 | 30 | 31 |
| 台形法 | 32 | 32 |
台形法の点は、連続な式の曲線とちょうど重なっています。
速度が一定の割合で変わる(速度のグラフが直線になる)場合、台形の面積は本当の面積と一致するので、台形法は誤差なく計算できます。
オイラー法と半陰的オイラー法は、どちらも最高点が2ピクセル(\(\ v_0/2\ \))ずれます。
これは、計算方法の近似による誤差で、打ち切り誤差と呼ばれます。
ずれはありますが、ゲームではたいてい、実際に動かしてみて初速度や重力を調整するので、ずれがあっても問題になりません。
重力が一定なら、台形法は \(\ y_{n+1}=y_n+v_n-\cfrac{g}{2}\ \) と同じなので、半分の重力を表す定数を1つ増やすだけで使えます。
飛び飛びの場合の正確な式
数値積分法の結果は、値を足し合わせたものなので、和の公式を使えば正確に求められます。
半陰的オイラー法で、初速度が重力のちょうど\(\ m\ \)倍(\(\ v_0=mg\ \))の場合、速度は\(\ m\ \)フレーム目に0になります。
そのときの高さが最高点です。
同じように、オイラー法では\(\ H=\cfrac{gm(m+1)}{2}=\cfrac{v_0^2}{2g}+\cfrac{v_0}{2}\ \)です。
連続な式との違い\(\ \pm\cfrac{v_0}{2}\ \)は、\(\ \displaystyle\sum_{k=1}^{m}k=\frac{m(m+1)}{2}\ \)と\(\ \displaystyle\int_0^m t\,dt=\frac{m^2}{2}\ \)の違いから来ています。
この式を使うと、ジャンプの高さを先に決めて、初速度を逆算できます。
たとえば、重力0.25で、高さ48ピクセル程度のジャンプを半陰的オイラー法で作りたいとします。
| \(\ m\ \) | 初速度\(\ mg\ \) | 最高点\(\ gm(m-1)/2\ \) |
|---|---|---|
| 19 | 4.75 | 42.75 |
| 20 | 5.0 | 47.5 |
| 21 | 5.25 | 52.5 |
初速度を5.0($0500)にすれば、20フレーム目に47.5ピクセルの高さに達します。
連続な式\(\ v_0=\sqrt{2gH}\ \)で求めた初速度約4.90($04E6)をそのまま使うと、最高点は約45.6ピクセルにしかなりません。
安定性
振動や回転のような動きでは、計算方法の違いが、ずれではなく「動きが壊れるかどうか」の違いになります。
原点のまわりを回る点\(\ (x,\ y)\ \)を考えます。
1フレームに\(\ k\ \)ラジアン回る等速円運動は、\(\ \cfrac{dx}{dt}=-ky,\ \cfrac{dy}{dt}=kx\ \)と書けます(三角関数の「円運動と単振動」)。
片方の座標だけを見ると、ばねにつながれた重りのような単振動です。
これを2つの方法で計算します。
- オイラー法: \(\ x_{n+1}=x_n-ky_n,\quad y_{n+1}=y_n+kx_n\ \)
- 半陰的オイラー法: \(\ x_{n+1}=x_n-ky_n,\quad y_{n+1}=y_n+kx_{n+1}\ \)(更新したあとの\(\ x\ \)を使う)
半径64、\(\ k=1/16\ \)で3周させたものです。破線は半径64の円です。
オイラー法では、1フレームごとに半径が\(\ \sqrt{1+k^2}\ \)倍になります。
1周(約100フレーム)で半径は約78になり、外へ広がり続けます。
半陰的オイラー法では、点は閉じた曲線(わずかに傾いた楕円)の上を回り続け、何周しても広がりません。
周期は\(\ \cfrac{2\pi}{2\arcsin(k/2)}\ \)フレームで、\(\ k=1/16\ \)なら約100.5フレームです。
\(\ k\ \)を2のべき乗分の1にすれば、掛け算は算術右シフトで済みます。
三角関数の表も乗算器も使わずに円運動を作れるので、古くから知られた方法です(ミンスキーの円描画アルゴリズム)。
振動する動きをオイラー法で計算すると、少しずつ振幅が大きくなります。
半陰的オイラー法では振幅がほぼ保たれるので、ゲームではこの順番(速度を先に更新する)がよく使われます。
丸め誤差の積み重なり
固定小数点数で足していくと、1回ごとの値の丸めによる誤差(丸め誤差)も、足した回数に比例してたまっていきます(固定小数点数の「誤差」)。
たとえば、0.1ピクセル/フレームの速度を8.8形式で表すと、$0019(0.09765625)になります。
600フレーム(NTSCで約10秒)動かすと、60ピクセルのはずが約58.6ピクセルしか進みません。
次のような方法で誤差を減らせます。
- 定数を、2のべき乗分の1の組み合わせで表せる値にする(0.25は8.8形式で誤差なく表せる)。
- 小数部のビット数を増やす(16.16形式なら、上の例の600フレーム後の誤差は約0.004ピクセル)。
- 足していく代わりに、\(\ x_0+nv\ \)のように開始の値から毎回計算する。
- 決まった位置で止めたい場合は、最後に目標の値をそのまま書き込む(線形補間)。
刻み幅
数値積分法の結果は、時間の刻み幅によって変わります。
同じ定数のまま、1フレームを2回に分けて半分ずつ計算すると、軌道も少し変わります。
たとえば、60Hz(NTSC)向けの動きを50Hz(PAL)で同じ時間の動きにするために、速度を6/5倍、重力を36/25倍にしたとします。
連続な式の最高点\(\ \cfrac{v_0^2}{2g}\ \)は変わりませんが、半陰的オイラー法の誤差\(\ \cfrac{v_0}{2}\ \)は速度に比例して大きくなります。
そのため、ジャンプの高さはわずかに低くなります。
より精度の高い方法
1回の刻みの中で加速度を何回か計算して平均する、ルンゲ=クッタ法などの方法もあります。
ただし、ゲームでは加速度が一定か単純なことが多く、計算量も増えるので、SNESではほとんど使われません。
指数関数的な近づき方
毎フレーム、目標の値\(\ T\ \)との差の1/8だけ近づける処理を考えます(線形補間の「目標へなめらかに近づける」)。
\[ x_{n+1}=x_n+\frac{T-x_n}{8} \]目標との差を\(\ d_n=T-x_n\ \)とすると、\(\ d_{n+1}=\cfrac{7}{8}d_n\ \)なので、差は\(\ d_n=d_0\left(\cfrac{7}{8}\right)^n\ \)と、等比数列になります。
これは、連続な式\(\ \cfrac{dd}{dt}=-\cfrac{d}{8}\ \)(差に比例した速さで差が減る)をオイラー法で計算したものです。
連続な式の答えは指数関数\(\ d(t)=d_0e^{-t/8}\ \)で、等比数列はその飛び飛び版にあたります。
8フレーム後の値は、等比数列で\(\ (7/8)^8\fallingdotseq0.344\ \)、指数関数で\(\ e^{-1}\fallingdotseq0.368\ \)で、少しずれます。
差が半分になるまでには約5.2フレームかかり、1フレームごとに約1.16dBずつ小さくなります(対数とデシベルの「指数関数的な減衰と直線的な減衰」)。
近づける割合を\(\ c\ \)とすると、差は1フレームごとに\(\ 1-c\ \)倍になります。
\(\ c\ \)が1を超えると目標を通り越して行ったり来たりし、2を超えると差が大きくなっていきます。
これもオイラー法の安定性の問題で、割合を大きくしすぎてはいけない理由です。