テイラー、フーリエ、球面調和関数


~ Taylor, Fourier, Spherical Harmonic ~






■はじめに

最近の知識の暗黙の前提にある級数展開にとり組んでみましょう。
級数などによる関数の展開は、うまい級数を選んでやると、 少ない基底関数で精度の高い関数の近似を行なってくれます。
級数展開がすらすらとできることは、将来のプログラマの必須条件になるかもしれません。

また、球面調和関数を見るためのプログラムを作って見ました。
qやwでlを、aやsで球面調和関数のmの値を変えることができます。

■テイラー展開(Taylor Series expansion)

関数f(x)=x は2次元グラフではまっすぐなグラフです。関数f(x)=x2は放物線を描いて、f(x)=x3はより複雑なグラフになります。

y=x y=x2 y=x3

これら、べき乗した変数のグラフは拡大縮小しても重なり合うことはありません。 これら変数のべき乗の組は、線形に合成することによって、より複雑なグラフを描く関数を作ります。 例えば、y = x3 - x2 - 2xは下のようなグラフです。

y = x3 - x2 - 2x

変数のべき乗の合成として作られる関数は、一般に

と書けます。これは、正に<1, x, x2, x3, x4, …>を基底ベクトルとするベクトルの姿です。このようにべき乗などの簡単な関数の組を基底ベクトルにつかって一般の関数を展開する方法を級数展開と呼びます。よく見る関数の多くは、べき乗で展開することができます。

そもそもべき級数で関数を展開する理由は、非常に複雑な関数の場合には値を求めるのが難しいので、 べき級数で展開して値を近似的に求めることによって、関数のおおよその形や性質をつかむためです。 それ以外にも、高速に計算するために計算が簡単なべき級数に落とし込むという側面もあります。

さて、a0, a1, a2, …はどのような値でしょうか?f(x)の関数の形が分かっていれば、aiを計算することができます。一度aiを求めておけば、どんなxに対してもすぐにf(x)の値を求めることができます。

aiを求める方法は、次のようになります。f(0)は、xiが必ず0になるので、

となります。従って、a0=f(0)です。

次に、f(x)を微分してみましょう。f(x)の微分f(1)(x)は、

になります。従って、a1 = f(1)(0)です。さらに微分していきましょう。f(x)の2階微分f(2)(x)は、

です。またx=0を代入すれば、a2 = (1/2)f(2)(0)を求めることができます。これを続けていくと、全ての係数は、

と、求めることができます。従って、関数なんてものは、べき級数を使って

と書くことができます。展開係数にx=0での微分を使ったこの展開方法をマクローリン展開(Maclaurin Series expansion)と言います。
さて、今は、f(0)と、x=0を使って係数を求めましたが、別に0でなくてもかまいません。 より一般的に、適当なxの値を使って係数を求める方法をテイラー展開といいます。 テイラー展開の形は次のようになります。

有名なテイラー展開には、次の物があります。

ちなみに、sin(x)を最初の項から順に計算していくと、

y = x y = x - x3/3 y = x - x3/3! + x5/5!

と、項が増えるに従ってだんだんとホントの値に近づいていくようすが見えます。
一般にテイラー展開は展開した位置の精度が高くて、xが展開した場所から離れるに従って、誤差が大きくなります。

■ホーナーの方法(Horner's method)

あんまり数学的なことを言っていると疲れるので、少し実際に使う時のテクニックも取り扱っておきましょう。最近の GPU 等では、積和計算

を高速に処理することができます。テイラー展開も積和計算を有効に使えれば高速化につながります。テイラー展開を変形すると、

とかけます。つまり、テイラー展開は、積和計算を入れ子にした形に変形できます。より具体的には、

として、何回か積和計算すれば、テイラー展開を好きな精度で実現することができます。この方法を入れ子乗算(nested multiplication)もしくはホーナーの方法と言います。

■フーリエ級数(Fourier series)

関数がべき乗で展開できたのなら、他の関数でも展開したくなるのが人情というものです。展開する基底ベクトルを整数倍の周期を持った三角関数(フーリエ級数)にした関数の展開がフーリエ展開です。フーリエ展開は、区間で区切られた関数、特に区間の端でつながった周期的な関数を展開するのに良く用いられます。その根底には、「同じ形を繰り返している周期をもった波はどんなに複雑なものも、単純な波たくさん足し合わさってできている」という思想が隠されています。

例えば、-XからXの範囲の値が分かっている関数f(x)を考えて見ましょう。この関数は三角関数を使って、

のように展開することができます。まぁ、どうして展開できるのかといえば、フーリエが証明したからで、ここでは展開できるということを事実として受け取って先に進みましょう。
あと、c0に掛けられた、わざとらしい1/2も、後々の係数を求めるときに綺麗な値になるように人為的に調整されています。
ちなみに、展開する三角関数は、項が増すに従って、間隔が短くなる波になります。 例えば、X=1の時のフーリエ級数のそれぞれの関数は次のようになります。

Y = cos(πx) Y = cos(2πx)
Y = sin(πx) Y = sin(2πx)

それでは、展開係数cn, snを求めましょう。ここで、三角関数に関して、次の式が成り立つことを利用します。

この式は、三角関数の和の公式を使って、証明することができます。これらを使うと、それぞれの係数は、

と求められます。

ためしに符号関数(xが負のときに-1正の時に+1)の関数をフーリエ展開してみましょう。計算をすると、cnの値は0になります。snは、

になります。従って、グラフを描くと、Xを1として、下のようなグラフがかけます。

Sinの1次の項まで Sinの3次の項まで Sinの5次の項まで Sinの7次の項まで

フーリエ級数は「波」を重ね合わせるので、展開の形も波打ったものになります。特に今回のようにx=0で不連続になるような関数は苦手で、不連続点の周辺では収束が悪くなります。

■球面調和関数(Spherical Harmonic)

いままでに、テイラー展開、フーリエ展開をしてきましたが、 級数 Fn(x)をいくつも持ってきて、それらを使えば、

と、関数を展開でき、さらに、この級数 Fn(x)が、

という性質を持っていれば、展開係数を

という形で計算することができます。
フーリエ展開;< 1, cos(πx/X), cos(2πx/X), ..., sin(πx/X), sin(2πx/X), ...> は、その様な級数の1つです。
そのような関数の1つで有名なものにルジャンドルの(陪)関数があります。ルジャンドル陪関数Pnm(x)は、 -1から1の範囲で有効な関数で、ルジャンドル関数Pn(x)から

として定義されます。
ちなみに、ルジャンドル陪関数 Pnm(x)の添え字、nとmは、0<=m<=nでなければなりません。n<mの時は、定義から0になります。

ルジャンドル陪関数は、次の性質を持っています(これらは、前の定義の式から計算することができます)。

この関係式から、しちめんどくさい微分を計算しなくても、ルジャンドル陪関数を求めることができます。
具体的には、1つめの式から、P00(x), P11(x), P22(x)…を求めて、次に2つめの式からP10(x), P21(x), P32(x)…を求めて、その後に3番目の式からそれ以外を求めます。
この計算方法は、特にコンピュータを使って求める場合に標準的な方法になるでしょう。

具体的にいくつかのルジャンドル倍関数を求めると、

になります。

さて、ルジャンドル陪関数の性質はいくつか扱いましたが、ルジャンドル陪関数を使った展開はどうなるのでしょうか?ルジャンドル陪関数の積分の性質は次のようになります。

nとlについては等しくなければ0という展開に使いやすい性質を持っています。
但し、これは、mの値は同じ時にだけ成り立つ式です。 mが違うときについては、一般には複雑な式になります。
といっても、mが違うときのことを考える必要はありません。なぜなら我々はるジャンドル倍関数をそのままでは使わないからです。
ルジャンドル関数を使うときには、フーリエ級数を併用します。
ルジャンドル陪関数にフーリエ級数を掛けて、新しい級数Ylm(θ,φ)を作ります。

このYlm(θ,φ)が球面調和関数です。いままでにxとして扱っていた変数は、cosθになりました。
cosθを引数とするので、ルジャンドル陪関数の引数の範囲は-1から1になります。これは、ルジャンドル陪関数の定義域と等しくなっています。
球面調和関数は、3次元球面の角度方向に関する展開をするために考えられました。半径1の3次元球面を考えてみましょう。球面上の点は

と、θとφの2つの変数で書くことができます。地球儀に対応させてみれば、緯度がθ、軽度がφになります。

このようなθとφでかける変数は、例えば地表の高さh(θ,φ)があります。地表の高さh(θ,φ)は、エベレストの山頂なら8850mとなるでしょう。また、マリアナ海溝なら-11000mになります。このようなθとφの関数を展開するのに適しているのが球面調和関数です。

球面調和関数は、主に物理学の分野で有名な級数で、電磁気学におけるアンテナが受け取る電磁波の志向性(どの向きからどのくらいの電磁波を受けるか)や、量子力学での原子中の電子の軌道を求めるのに使います。電子の軌道を求める問題は、量子力学の山場のひとつで、多くの学生が球面調和関数にたどり着く前に数式の複雑さに睡魔に降伏してしまいます(だめじゃん)。原子核の周りを電子が飛んでいる絵を見かけたことがあると思います。実はこの電子が飛んで行る軌道は球面調和関数(の2乗)なのです。こう思うと球面調和関数に親しみがわくのではないでしょうか?(ホントか?)球面調和関数は現在の物理学では量子力学で主に勉強するので、より詳しく知りたければ、量子力学岩波基礎物理シリーズ (5)の量子力学などの量子力学の本を読んでみるといいと思います(どうして、わざわざこの本を出したかって?それは私の名前が載っているから…)

球面調和関数同士の積分は、

になります。この関係式を用いると、3次元球面上で定義された関数f(θ,φ)は球面調和関数によって、

と、展開されます。球面調和関数のそれぞれの具体的な形は次のようになります。見て分かるとおりlが増えるに従ってどんどん複雑な関数になっていきます。

(Robin Green, Spherical Harmonic Lighting: The Gritty Details から引用)

球面調和関数を使って、丸まったオブジェクトを展開すると、最初は丸いオブジェクトですが、lが増えるに従って本来の形に近づいていきます。

(Robin Green, Spherical Harmonic Lighting: The Gritty Details から引用)

■最後に

さて、こんな難しい球面調和関数ですが、何の役に立つのでしょうか?
最近の使い道は、光源が動く場合の静的なオブジェクトへのライティングです。
光源の方向を(θ,φ)としたときに、グローバルイルミネーションなどを使って、ライティングした結果f(θ,φ)をピクセルごとに計算するアルゴリズム(BRDF)をあらかじめ決めておきます。
この複雑な計算はリアルタイムには実行できないので、明るさの球面調和関数に関する展開係数を各ピクセルごとにテクスチャとして、格納しておきます。
レンダリング時には、球面調和関数を頂点単位(ピクセル単位でもできるかも)に求めて、テクスチャとして格納された係数をかけてライティングを計算します。
早い話、遮蔽マップの1つなのですが、光源が動いたときの変化が比較的なだらかであり、拡散光を少し変更した程度のライティングなら、テクスチャの枚数が少なくても実現できるところがもてはやされる理由でしょうか。





もどる

imagire@gmail.com