差分方程式
前の値から次の値を決める式(漸化式・差分方程式)の解き方と、その値が落ち着くか、振動するか、大きくなり続けるか(安定性)の見分け方について説明します。
<注意>
このページでは、数式を表示するためにMathJaxを使用しています。ブラウザの設定でスクリプト実行をブロックしている場合、ページ内の数式が正しく表示されません。
差分方程式とは
毎フレーム、前のフレームの値から次の値を計算する処理は、ゲームのいたるところにあります。
このような「前の値から次の値を決める式」を差分方程式(漸化式)と言います。
\(\ x_{n+1}\ \)を決めるのに\(\ x_n\ \)だけを使うものを1階、\(\ x_n\ \)と\(\ x_{n-1}\ \)を使うものを2階の差分方程式と言います。
このページでは、\(\ f\ \)が掛け算と足し算だけでできている(線形な)場合を扱います。
差分法と数値積分法で扱った、毎フレーム速度を足す計算や、差の1/8ずつ近づける計算も差分方程式です。
差分方程式の解き方が分かると、何フレーム後にどうなるかを、1フレームずつ計算しなくても求められます。
1階の差分方程式
掛け算だけの場合
\[ x_{n+1} = a\,x_n \]1回ごとに\(\ a\ \)倍されるので、\(\ x_n = a^n x_0\ \)です(等比数列)。
\(\ a\ \)の値によって、振る舞いが大きく変わります。
| \(\ a\ \)の値 | 振る舞い |
|---|---|
| \(\ 0 < a < 1\ \) | 同じ符号のまま0に近づく。\(\ a\ \)が0に近いほど速い |
| \(\ -1 < a < 0\ \) | 符号が毎回入れ替わりながら0に近づく(振動しながら落ち着く) |
| \(\ \lvert a\rvert > 1\ \) | 大きくなり続ける(発散する) |
| \(\ a = 1\ \) | 変わらない |
0に近づく速さは、\(\ |a|^n\ \)が半分になるまでの回数\(\ \log(1/2)/\log|a|\ \)で表せます。
対数とデシベルで言えば、1回ごとに\(\ 20\log_{10}|a|\ \)dBずつ小さくなります。
定数を足す場合
\[ x_{n+1} = a\,x_n + b \]\(\ x = ax + b\ \)を満たす値\(\ x^* = b/(1-a)\ \)(不動点)からのずれを\(\ d_n = x_n - x^*\ \)とすると、\(\ d_{n+1} = a\,d_n\ \)になります。
定理
\[ x_n = x^* + (x_0 - x^*)\,a^n, \qquad x^* = \frac{b}{1-a} \]
\(\ a \neq 1\ \)のとき、\(\ x_{n+1} = ax_n + b\ \)の解はである。\(\ |a| < 1\ \)なら、\(\ x_n\ \)は\(\ x^*\ \)に近づく。
例 (目標へなめらかに近づける)
\(\ x_{n+1} = x_n + (T - x_n)/8\ \)は、\(\ a = 7/8,\ b = T/8\ \)の場合で、不動点は\(\ x^* = T\ \)である。
目標との差は、毎フレーム7/8倍になる(線形補間の「目標へなめらかに近づける」)。
例 (空気抵抗のある落下)
毎フレーム、速度に重力\(\ g\ \)を足してから\(\ k\ \)倍する(\(\ 0 \lt k \lt 1\ \))と、\(\ v_{n+1} = kv_n + kg\ \)になる。
速度は\(\ v^* = kg/(1-k)\ \)に近づくので、最大の落下速度を決めなくても、自然に一定の速さに落ち着く。
2階の差分方程式
\[ x_{n+1} = p\,x_n + q\,x_{n-1} \]\(\ x_n = r^n\ \)の形の解を探すと、\(\ r\ \)は次の2次方程式(特性方程式)の解になります。
\[ r^2 = p\,r + q \]定理
特性方程式の解を\(\ r_1,\ r_2\ \)(\(\ r_1 \neq r_2\ \))とすると、一般の解は\(\ x_n = c_1 r_1^n + c_2 r_2^n\ \)である。
\(\ c_1,\ c_2\ \)は、\(\ x_0\ \)と\(\ x_1\ \)から決まる。
特性方程式の解が実数なら、1階の場合と同じように近づいたり発散したりします。
解が複素数(\(\ r = |r|e^{\pm i\theta}\ \)、複素数)のときは、\(\ x_n\ \)は\(\ |r|^n\ \)倍されながら、1回に角度\(\ \theta\ \)ずつ進む正弦波になります。
\(\ |r| < 1\ \)なら揺れながら小さくなり、\(\ |r| = 1\ \)ならいつまでも同じ大きさで振動し、\(\ |r| > 1\ \)なら揺れがどんどん大きくなります。
図は、どれも1周期16回(\(\ \theta = 2\pi/16\ \))の場合です。
正弦波を作る
\(\ |r| = 1\ \)になるように\(\ p = 2\cos\theta,\ q = -1\ \)とすると、次の式で正弦波が作れます。
\[ x_{n+1} = 2\cos\theta \cdot x_n - x_{n-1}, \qquad x_0 = 0,\ x_1 = \sin\theta \]このとき\(\ x_n = \sin n\theta\ \)になります。
三角関数の表を引かずに、掛け算1回と引き算1回で次の値が求まります。
ただし、固定小数点数で計算すると、\(\ 2\cos\theta\ \)の丸めや計算の丸めで\(\ |r|\ \)がわずかに1からずれ、長く続けると振幅が少しずつ変わっていきます。
安定性
ここまでの結果は、次のようにまとめられます。
定理 (安定性)
線形な差分方程式の解が、どんな初めの値からでも0(または不動点)に近づくのは、特性方程式のすべての解の絶対値が1より小さいときである。
複素平面で言えば、特性方程式の解が、すべて半径1の円(単位円)の内側にあれば安定です。
1つでも外側にあれば、値は大きくなり続けます。
図の×は、このページのグラフの例の解の位置です。円の内側の青は落ち着き、円の上の緑は振動が続き、外側のオレンジは発散します。
差分法と数値積分法の「安定性」で、オイラー法で円運動を計算すると半径が大きくなり続けたのも、計算を差分方程式として見ると、解の絶対値が\(\ \sqrt{1+k^2} > 1\ \)になっているためです。
入力のある差分方程式
毎回、外から値\(\ u_n\ \)が加わる場合を考えます。
\[ y_n = a\,y_{n-1} + u_n \]音の処理では、\(\ u_n\ \)が入力のサンプル、\(\ y_n\ \)が出力のサンプルで、これはIIRフィルタそのものです(サンプリングとデジタルフィルタ)。
入力に、1回だけ1でほかは0の値(インパルス)を入れると、出力は\(\ 1,\ a,\ a^2,\ \dots\ \)と続きます。
出力が過去の出力に依存しているので、インパルスを1回入れただけでも、出力はいつまでも続きます。
エコー
DSPのエコーは、\(\ D\ \)サンプル前の出力に係数\(\ g\ \)を掛けて、入力に足し戻しています(実際にはその間にFIRフィルタも入ります)。
\[ y_n = u_n + g\,y_{n-D} \]1つの音は、\(\ D\ \)サンプルごとに\(\ g\ \)倍されながら繰り返されます。
\(\ |g| < 1\ \)なら小さくなって消えていきますが、\(\ |g| \ge 1\ \)にすると、音が大きくなり続けて割れてしまいます(エコーを使う)。
特性方程式は\(\ r^D = g\ \)なので、解の絶対値はどれも\(\ |g|^{1/D}\ \)で、\(\ |g| < 1\ \)が安定の条件になります。
入力のある差分方程式の性質は、Z変換を使うと、まとめて調べられます。
整数で計算するときの注意
差分方程式を整数(や固定小数点数)で計算すると、小さな値が丸めで消えて、計算どおりにならないことがあります。
目標100へ差の1/8ずつ近づける計算を、差を右に3ビットシフトして整数で行うと、差が7以下になったところで1/8が0になり、93で止まってしまいます。
小数で計算すれば、目標に限りなく近づいていきます。
このような「動かなくなる範囲」を作らないためには、小数部を持たせて計算するか(固定小数点数)、差が小さくなったら目標の値をそのまま書き込みます。
反対に、丸めの向きによっては、0に落ち着くはずの値が、0の周りで小さく振動し続けることもあります。