計算機シミュレーションの方法と具体例

このモデルによるシミュレーションの身近な例として、待ち行列のシミュレー ションがある。これは、たとえば健康診断のときに身長計、体重計、視力検 査装置などを何台ずつ用意し、どのような順番でならべたら最も順番待ち時間 が短くできるかといったことを検討したりするときに採用されるモデルである。 計算機の中に個々のサービス窓口(体重計)を設定し、サービスを受ける人が 健康診断場に到着する時間間隔がある確率分布にしたがうと仮定してシミュレー ションを行う。窓口の数や、その順番をいろいろ変化させてある人数の人が 健康診断を終えるまでの平均時間を調べ、最適な組み合わせを求めることがで きる。ある人が到着してから次の人が到着するまでの時間間隔は乱数を用いて 決める。次の人はすぐにやってくるかもしれないし、たまたましばらく時間が 空くこともあるだろう。この時間間隔をいわばサイコロを投げて決めるわけで ある。このようなシミュレーション方法のことをモンテカルロ法とよぶ。

計算機の中でサイコロを振ることに相当するのが乱数である。その発生方法と してはある数をつぎつぎと掛けていって、別のある数で割った余りを取る方法など がよく用いられる。このようにして、たとえば 0 から 255 までの整数を同じ 確率で発生させることができる。乱数といってもまったくでたらめな数なので はなく、ある規則にしたがって生成されたものなのである。しかし生成された 数が1の次には2が出やすいといった性質を持っていたりするとこれは乱数とし て採用することはできない。すべての数が同じ確率で出現し、しかも互いに 独立(次に出る数が前に出た数に依存しない)とみなしうるものでなければならな い。さまざまな統計的検定に耐え、かつ計算機で高速に生成できるような乱数 発生方法がいろいろ工夫されている。通常、すべての数が同じ確率で発生する 一様乱数を生成し、これをもとにしてさまざまな確率分布の乱数が導かれる。 たとえば、多くの一様乱数の和をとると、その分布は正規分布に近づくことが 知られている。

モンテカルロ法は、さまざまな物理現象のシミュレーションに利用されている。 たとえば電子計算機の開発者としても知られている von Neumann はモンテカ ルロ法を用いて物質中の中性子の運動を解析した。私自身も、高温電離大気中 の光子の運動をモンテカルロシミュレーションで解析したことがある。中性子 星やブラックホールのまわりのガス円盤のように 10憶度もの温度がある大気 中では光子は電子と衝突する際に電子からエネルギーを受け渡されて高エネル ギーのX線やガンマ線として放射されることになる。光子と電子の速度、運動 方向などを乱数を用いて決め、衝突後の光子が電子との衝突をくり返しながら ガス円盤から飛び出していくまでの過程を計算機で追跡した。1万個程度の光 子を用いると円盤から放射される光子のエネルギー分布(スペクトル)をかなり 正確に求めることができる。シミュレーションから予測されるX線、ガンマ線 のスペクトルをX線天文衛星によって観測されるブラックホール候補天体のス ペクトルと比較することにより、これらの天体の物理状態(温度、密度など)を 知ることができるわけである。

モンテカルロシミュレーションは物質の性質を 探る研究(物性物理学)の分野でも威力を発揮している。たとえば、鉄のような 強磁性体の温度を上げていったときに、最初はある方向に揃っていた原子の スピン磁気モーメントの方向がしだいにランダムになっていき、磁石としての 性質が失われていく。相転移とよばれる、このような過程を原子は2次元格子 状にならんでいるとして、ある原子のスピンの方向(上向き、または下向き)は それと隣接する原子との相互作用のみによって確率的に決まっているという 単純化したモデル(イジングモデル)でシミュレートすることができる。ある原 子のスピン方向が反転するかどうかは、この確率分布にしたがう乱数を発生 することで決める。

これはジェット機のまわりの空気の運動などをシミュレートする際に用いられ るモデルで、空間を格子に分けて各格子点における密度、速度、温度などの時 間発展を計算するという手法である。たとえば全地球的な天気予報をする場合、 地球を経度方向、偉度方向、高さ方向それぞれについて数百個の格子に分割し、 それぞれの点での現在の密度等の情報を与えて、少しだけ時間がたったときの 密度などを大気の運動を記述する偏微分方程式(基本的には数値風洞実験で用 いられるのと同じナビエストークス方程式であるが地球の回転の効果が含めら れている)を数値的に解くことによって次々と求めていく。方程式が偏微分方 程式になるのは、空間的に固定されたある領域を考えたとき、大気の運動によっ てその領域に入ってくる空気と出ていく空気があり、その領域の内と外との間 の物理量の差が重要になってくるからである。たとえば、南から暖かい空気が 流れこんでくれば、いま考えている領域の温度は上昇する。時間的な変化を計 算するときに空間的な変化も考慮する必要があることから、方程式が偏微分で 表現されるのである。偏微分方程式をコンピュータで解くにはいろいろ難しい 点があって以前は数値的不安定性によって計算が止まってしまうということが よくあった。この10年の間にさまざまな改良が施され、現在では、それぞれの 問題に応じてどのような方法を選択すればよいかがよくわかってきている。 基本的には空間微分を差分におきかえる、すなわち隣接する格子点の物理量の 差をとることによって空間微分を計算し、時間方向に積分していく。

3次元空間を格子に分けて計算する場合、1方向について100個の点を取ると 全体では 100 X 100 X 100 = 100万個の点が必要になる。1点について微少 時間後の物理量を求めるのに 100回の演算が必要とし、時間方向について たとえば 10秒きざみで計算して 1日 (86400秒)追跡するとすると全部で 8640000000000 回の演算が必要となる。これが、数値天気予報などに、 きわめて高速な計算機が必要とされる理由である。1秒間に1000万回の 計算しかできない通常のワークステーションでは、明日の天気を計算するのに 10日かかることになって実用にならない。最近では 1000 X 1000 X 1000 格子 の計算を可能にすることを目標として次世代のシミュレーション用計算機の 開発が進められている (注：2002年に完成した地球シミュレータを使えば このような計算が可能になった)。

連続モデルによるシミュレーションの例として我々のグループが行っている、 天体プラズマの数値シミュレーションを紹介しよう。天体プラズマでは磁場が 運動を支配する重要な要素になることがある。たとえば太陽活動が 11年の周 期を持っていることはよく知られているが、これは磁場によって起きている。 活動の激しい時期には表面に多くの黒点があらわれ、また、フレアという 爆発現象が起きて、オーロラや磁気あらしといった形で、地球にもその影響が 及んでくる。このような太陽表面の活動現象は、太陽内部の対流層で、発電作 用（ダイナモ）によって生成された磁場が、表面に浮上してくることによって 生じていると考えられている。磁場が浮上するといってもイメージがわきにく いかもしれないが、水素やヘリウムのように軽いガスをつめた風船が大気中を 浮き上がっていくのと同様に、磁束管中のガス密度がまわりよりも小さいと、 浮力の作用を受けて、太陽大気中を浮上していくのである。このような磁束が 太陽表面から顔を出した部分が活動領域のもとになる浮上磁場領域である。 我々は太陽対流層からコロナに至る領域を 64 X 64 X 150 の格子に分割し、 太陽内部の対流層の中で捻られた磁束菅が表面に浮上していく様子の3次元 シミュレーションを行った。その結果、日本の太陽X線観測衛星「ようこう」 で観測された太陽コロナの構造によく似た磁気ループが形成されること、 磁気ループの根元に垂直方向の強い磁場を持つ太陽黒点に類似した構造が できること、また、それらの黒点の運動が観測される運動と一致していること などを示すことができた。また、捻れた磁束菅と磁束菅の相互作用のシミュレー ションも行ってみた。ふたつの磁束菅がぶつかったところで磁力線のつなぎか え(リコネクション)がおこり、磁束菅に沿った方向に高速ジェットが飛び出す ことがわかった。これらのシミュレーション結果は「ようこう」が観測した 太陽フレアや X線ジェットなどの観測結果とよく一致している。

磁束菅の浮上という現象は太陽大気だけではなく、銀河系の中でも起きている。 シミュレーションによって、最初銀河の回転方向を向いていた磁力線が銀河 面をつらぬいて湾曲し、ガスは磁力線の谷底に落下して集められることがわかっ た。たとえばオリオン座における星形成の母体となっている巨大分子雲は、 このメカニズムによって形成された可能性がある。さらに、シミュレーショ ン結果は、この高密度領域の上に銀河面に垂直に長く延びた密度の高い領域が できることを示している。銀河面におけるガス分布を観測すると実際にこのよ うな構造が見られ、シミュレーション結果は現実の銀河の姿をよく反映してい る。その他、ブラックホールに落下するガス円盤からジェットが形成される 様子など、さまざまな天体現象を計算機の中に再現することができるように なっている。