./note006.html

ToC

5. 楕円型方程式

楕円型方程式は、典型的にはラプラス方程式

またはポアソン方程式

のようなものである。上のポアソン方程式はニュートン重力の場合で、

は重力定数、

は物質の質量密度（単位体積あたりの質量）である。ここで、

はとりあえず与えられた空間座標だけの関数とする。空間1次元ではラプラシアンがただの2階微分演算子になってしまい、解くべき方程式は常微分方程式ということになる。

実際にポアソン方程式（とさらに状態方程式を連立させてでてくるレーン・エムデン方程式）を球対称の場合に解くというのは I のほうでやったと思うので、ここでは空間2次元の場合を考える。3次元も考え方は同様である。

工学的な応用では、楕円型方程式の応用として重要なのはラプラス方程式の境界値問題である。これは、構造解析、電場の解析などあらゆるところで現れる。特に複雑な形状の機械部品とか機械全体の中での応力場の解析といったものが重要な応用ということになる。

が、天文・物理のいろんな問題では、複雑な境界値問題というのはあんまり大事ではない。天体は大抵重力でまとまっているわけで、固定された表面（境界）というものが与えられているとは限らないからである。もちろん、数値計算の方法によってはそういうものが出てくることもあるが、とりあえず今日は比較的に単純な境界条件だけをかんがえることにしよう。

ラプラス方程式やポアソン方程式は、基本的には拡散方程式から時間微分の項を落したものと考えることができる。つまり、なんらかの定常状態を記述しているわけである。多くの応用で、実際にそのような定常状態を求めているわけであるが、数値計算という観点からすると、ラプラス方程式やポアソン方程式に適用できる計算法は実は拡散方程式や波動方程式に陰解法を適用する時にも使える。

これは、前回見たように陰解法では空間方向の差分方程式を解く必要がでてくるからである。解くべき方程式は、ラプラス方程式やポアソン方程式の差分化からでてくるものと基本的には同じである。

というわけで、今日は、単純な2次元ポアソン方程式の境界値問題を例に、計算方法について考えていく。

2次元なので、ラプラシアンは

である。例によって変数の方を

と書き、これが

の近似解としよう。また、前回と同じ2次精度の中心差分を使うとすれば、ラプラシアンの差分近似は以下で与えられる。

但し、ここでは

とした。いま、解くべき問題を長方形領域でのポアソン方程式の境界値問題

ということにすれば（以下、面倒なので

の項は 1 とする）、これは2次元配列 u[nx][ny] についての線形連立方程式ということになる。ここで

であり、境界上では

, それ以外では

となる。

5.1. 線形方程式と反復法

さて、というわけで、2次元の長方形領域上でのポアソン方程式の差分近似は、連立線形1次方程式 ( 41 )になるので、これを何らかの方法で解けばいい。

空間1次元の場合には、出てくる行列が3重対角という特別な性質を持ったものであったので、自由度（要素数）に比例する計算量やメモリー使用量で方程式を解くことができた。

2次元の場合はどうなるかというと、そんなにうまくはいかないということがわかる。もちろん、2次元の場合でも、というように添字をふりなおして(以下、添字が1つしかないときにはこのように付け変えたものを考えるということにする)、以下のように行列の形に書くことはできる。

ここで、

は

のベクトル表示である。しかし、

の要素を考えてみると、対角要素、その両側の他に、

だけ離れたところに 0 でない要素が現れる。

このために、計算量としては ( として)になり、1次元のときよりずっと計算量が多い。例えば、つまり正方形領域とすれば計算量は結局になってしまう。

もちろん、これでも一般の場合のガウスの消去法の計算量であるよりは少ないが、しかし非常に多いことには変わりはない。また、3次元の場合にはこの状況はいっそう悪化する。

このために、ガウスの消去法のようにまともに行列を解くのではなく、適当な近似を繰り返していくことで解を求めることはできないかということが昔から考えられてきた。以下、その方法のいくつかを紹介する。これらを反復法という。

5.1.1. ガウス反復

なんらかの方法で適当な近似解（ここで上つき添字は番目の近似解ということで、別にの乗という意味ではない。下につけると他の添字と混乱するので上につけてみた）があるとする。これをもうちょっとましな解にできないかということを考える。

一つの方法は、行列 A を形式的に以下のように分解することである。

ここで、

は対角要素

だけをとりだした行列で、

はそれ以外の残りである。これを使って、以下のような式を考える

は対角行列なので、これは解けていて、実際のプログラムとしては、

 for (i=1; i<nx-1;i++){
     for (j=1; j<ny-1;j++){
         unew[i][j] = 0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j])
                     - rho[i][j]*0.25*deltax*deltax;
     }
 }

ということになる。1で始まって

等で終っているのは、境界は 0 固定だからである。その分の配列を節約することも考えられるが、いじましいしちょっとわかりにくくなるので上の形に書くことにしよう。

この反復で、収束すれば、つまり、となれば、これらが元の方程式を満たしていることは明らかであろう。収束するかどうかは面倒なのでここでは省くが、上のような簡単なポアソン方程式の離散化で妙な非線形項とかがなければ収束するということが証明できる。

実用上は、収束したかどうか判定する必要がある。これには、普通はを計算してそれが適当な設定した値よりも小さくなったらいいことにする。

理論的には、数列は連続する2項間の差が小さくなっても収束に近付いているとは限らないが、ガウス反復の場合には、丸め誤差がなければこれで判定できることが証明できる。

が、この方法には2つ問題がある。一つは収束が非常に遅いことであり、もう一つは u と unew で配列が2ついるのでメモリがたくさん必要になることである。

5.1.2. ガウス・ザイデル法

さて、実は、上の2つの問題は一度に解決することができる。上のプログラムを、以下のように変えてしまうのである。

 for (i=1; i<nx-1;i++){
     for (j=1; j<ny-1;j++){
         u[i][j] = 0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j])
                    - rho[i][j]*0.25*deltax*deltax;
     }
 }

なにをしたかというと、 unew というのをやめにして u にしただけである。つまり、一回の反復の中では u の要素を順に変えていくわけだが、その時に、既に新しい値が求まっていればそっちを使おうというだけのことである。新しい値のほうが、多分より正確であろうと考えられるからである。

形式的には、これは以下のように書くことができる。行列を、対角成分、非対角成分の上半三角分、残りの3個に分ける ()と、

ということになる。

5.1.3. 収束の速さ

収束の速さをガウス反復について検討する。面倒なので1次元として、

に対する差分近似

を考える。ここでは、収束の速さが問題なので、

としてよくて、

である。さて、これは線型の式なので固有値と固有ベクトルがあって、ある固有値

と固有ベクトル

が

という関係を満たすと考えられる。固有関数がどういうものかは自明ではないが、等間隔格子なのでフーリエ級数解が固有関数になっていることを期待すると、境界条件もいれて

である。ここからの収束の議論では

は虚数単位で

が格子のインデックスとする。

はフーリエモードのインデックスである。これを 49 にいれて整理すると

となることがわかる。これは、

を与えると

の時にもっとも1に近く、収束が遅いことがわかる。また、反復回数は

程度必要であることもわかる。

なお、式 48 は、前節で扱った拡散方程式の陽解法差分近似とみなすこともできることに注意して欲しい。なので、上のように天下りに固有関数を仮定しなくても、前節と同様に線型差分方程式であることから解を求めることもできる。線形差分であることから指数関数型に限られるため、周期解はフーリエ級数解なので、2つの方法は数学的には等価である。

また、陽解法の安定性条件を満たすように、となっている。収束までは時間が1より十分大きくなる必要があるので、反復回数が程度必要となる。

ガウス・ザイデル法についても同様に必要な反復回数を求めることができる。次に述べる SOR についても(多分)できる。

5.1.4. SOR

というわけで、ガウス・ザイデルはガウス反復よりはましだが、まだもうちょっとなんとかならないものかという気もする。実際に近似解の挙動を見るとわかるが、どちらの方法でも、誤差の減り方にはパターンがあって、1方向から真の解に近付いている。

というわけで、以下のようなことが考えられる：ガウス・ザイデル法で、修正量に適当な係数を掛けて水増ししてやればもうちょっとうまくいくのではないか？

こんな安直な方法で大丈夫なのかと思うであろうが、うまくいってしまうのがすばらしいところである。プログラムとしては、ガウス・ザイデルの時の式をちょっと変えて、

        u[i][j] += omega*(0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j])
                    - rho[i][j]*0.25*deltax*deltax -u[i][j])

というだけで済む。実際的な問題は、

の値をどうとるかということである。理論的には

でないといけないことがわかっていて、さらに本当はもちろん線形のポアソン方程式とかなら最適な

の値が計算できるが、ここではそこまでは立ち入らないことにする。

この方法を SOR (Successive Over Relaxation, 本によっては Simultaneous ... となっているものもあるかもしれない)という。

5.2. もっと高度な解法

天文等で出てくるポアソン方程式の場合には、特にうまくを決められれば SOR で非常にうまくいくことが多い。が、まあ、世の中はそういう問題ばかりでもないので、それ以外の方法を紹介だけしておく。

5.2.1. ADI

ADI は alternating direction implicit の略である。ここでの基本的な考察は、「3重対角なら簡単に解ける」ということである。これを使って、方向と方向を交代で解いているうちになんとかなって欲しいというのが ADI の考え方ということになる。

もうちょっとちゃんと書いてみると、例えば方向に解くときには、

となる。ここで、 new とつけたのが更新後の値で、これを

から順に解いていく。

方向が終ったら次は

方向をやる。

ただし、実際に使う上では、 SOR と同じように、修正量の補正をしたほうが収束が速い。 SOR とは逆に修正量をすこし減らす必要があり、どれくらい減らすべきかを決めるのにはまたややこしい理屈がいる。

ただし、適当に修正量を決めても最適な SOR よりも遅くはない程度で収束するという利点がある。

5.2.2. 共役勾配法

共役勾配法は、実は、現在大規模な行列を解くのにはもっとも広く使われている方法である。これは、特にいろいろな前処理と組み合わせることで、例えば拡散方程式の定常問題で拡散係数が空間に強く依存する場合にも収束させられるとか、刻みを場所によって変えるような高度な方法でも同様になんとかできるとかいう利点があるからである。

しかし、これは安直にプログラムを書いただけではこれまでの紹介してきた方法に比べて速くなくて、いろいろ高度な変形をして初めて速くなる。また、ほとんど共役勾配法だけについて議論した良い教科書がいくつもあるので、ここでは省略する。ここでは、とりあえずそういう名前のものがあるということくらいをちょっと憶えておいてほしい。

とはいえ、極めて重要な方法ではあるので、これについては、もうちょっと数学的な準備をした後でもう一度やることにする。

ポアソン方程式を解く方法の一つで、応用上は重要なのはフーリエ変換を使った方法である。いま、固定境界をやめて周期境界条件でいいことにすると、ポアソン方程式の場合、をフーリエ級数展開して、出てきた係数に ( は波数)を掛けてから逆フーリエ変換すればポアソン方程式の解になっているわけである。固定境界なら関数だけで展開すればよい。フーリエ変換の計算量は多次元でも格子点の総数をとしての程度なので、これは反復法でいうと回数がせいぜいで済むということに相当する。