計算機物理学^{∗}

溜渕継博＠静岡大学理学部 tamari@sci.shizuoka.ac.jp

1 第 3 の物理学

前章までにも述べられているように，物理学の体系はいくつかの原理や法則を出発点として構成されていて，物理学は自然界の多様な現象を小数の原理や法則により如何にして統一的に説明するかということをその目的の一つとしている。この際に，物理学ではまずなんらかの理論的な仮説をたてて，その仮説に基づいて予想される現象を実験的に検証するのが普通である。もしも，実験結果が予想と反した場合にはそこで用いられた仮説は誤りであると結論されるし，予想と一致した場合にはその仮説はその現象を説明する理論の “ 一候補” となる。様々な現象にたいして様々な理論が提出され，実験的な検証を経て理論は淘汰されてきた。このようにして仮説のいくつかは原理や法則として採用されて，物理学が体系化されてきたわけである。

自然科学の中でも物理学は最も高度に発達した学問分野の一つといわれ，そこで用いられる方法論や概念の多くは工学や自然科学の他の分野ばかりでなく，社会科学の分野でも利用されて広い応用範囲を持っているが，近年では，物理学の研究で用いられる実験設備や理論体系が高度にしかも先鋭的に進歩したために，かってのように一人の物理学者が実験と理論の両方をこなすことがほとんど困難な状況になっている。現在，物理学は理論物理学と実験物理学に大別されて，研究はそれぞれの専門家によってなされるのが普通である。

ところで近年，「計算機物理学」あるいは「計算物理学」という言葉が理論物理学と実験物理学に並んで市民権を得つつあることをご存知だろうか。後に述べるように，計算機物理学では「紙と鉛筆」だけを用いた狭い意味での理論物理学では手に負えない計算をコンピュータを利用して行い，従来の意味での実験設備を用いた実験物理学では達成できない実験状況をコンピュータ上に出現させることができる。このような意味で，計算機物理学はある面で理論物理学的であるし，実験物理学にかなり似た面をもつ場合もあって，従来の理論と実験という分類の枠をはみ出しているといえる。コンピュータの誕生以来，コンピュータの発達とその物理学への応用は切っても切れない関係にあった。従来でも上に書いた意味でのコンピュータの物理学への応用は盛んであったが，それではなぜ今更「計算機物理学」なのだろうか。その答えはコンピュータの発達と無関係には語れない。

実験でのオンラインデータ処理や測定データの整理などへの利用を除くと，物理学におけるコンピュータの利用は，従来，どちらかというと，理論物理学の補助的な意味合いが強かったと思われる。しかし，最近のコンピュータの発達は著しく，数年前には全く不可能であった計算が最新のコンピュータと最新のアルゴリズム (計算の基礎となる考え方のこと) を利用することで可能になったり，現在不可能と思われている計算も数年後には (あるいは明日にも?) 可能になる場合もある。まだ一部の分野ではあるが，最近では計算機物理学の意味でのコンピュータの利用がかなり実用的なレベルに達しつつあり，この傾向は更に他の分野にも急速に広がっていくことが容易に想像される。このような状況と共に，最近では専らコンピュータを利用して物理学の研究を行う計算機物理学の専門家と呼べる人達が現れ始め，計算機物理学が理論物理学と実験物理学と並ぶ第三の物理学として認知され始めたわけである。

0 ’ ９３年度静岡大学理学部公開講座「自然のしくみ−物理の世界−」第４章より

この講義では最近の計算機物理学の考え方と成果のほんの一端を紹介すると共に，実習を通して計算機物理学の雰囲気を味わっていだければ幸いである。

2 コンピュータの歴史

最近のコンピュータの発達には目を見張るものがあって，持ち運びが可能なノート型のコンピュータでさえ，十年前の超大型コンピュータに匹敵する処理速度をもつものまで現れてきているのは驚きである。ここでは計算機物理学の「道具」であるコンピュータの歴史を簡単に振返ってみよう。

計算を機械に代行させようというアイデアはかなり古く， 17 世紀には「パンセ」や「パスカルの法則」で有名なフランスの哲学者 · 数学者 · 物理学者であったパスカル (1623–1662) が歯車式の計算機械を考案した。また，本格的なものとしてはバベッジ (1791–1871) による「差分機関」や「解析機関」が名高く，特に解析機関は，機械式であることを除くと現代のコンピュータと非常に近い構成をもち，バベッジは計算機の歴史の上で偉大な先駆者として称えられている。しかし，これらはどちらも当時の機械加工技術が十分でなかったなどの理由で実用には至らなかった。

人類初の電子式自動計算機，つまりコンピュータはエッカートらによってアメリカで開発された ENIAC (1942–1946) である。 ENIAC は約 18,800 本の真空管を用いて製作され，弾道力学の計算問題を計算専門家 (7 時間) の 8,400 倍 (3 秒) の速さで解いて当時の人々を驚かせたといわれている。日本最初のコンピュータ FUJIC が誕生したのは ENIAC から 10 年後の 1956 年のことであった。 ENIAC や FUJIC など，真空管を用いた初期のコンピュータは第 1 世代コンピュータと呼ばれ，以降，使用している素子に応じて

  第 1 世代 (  ∼ 1960) 真空管
第 2 世代 (1960 ∼ 1965) トランジスタ
第 3 世代 (1965 ∼ 1970) IC (集積回路)
第 3.5 世代 (1970 ∼ 1980) LSI (大規模集積回路)
第 4 世代 (1980 ∼  ) 超 LSI (超高密度集積回路)

などと分類さている。なお， 1982 年より 10 年計画でスタートした，いわゆる「第 5 世代コンピュータ」開発プロジェクトにおける “第 5 世代” という呼び名は上に挙げた第 1∼4 世代コンピュータとは異なる新しい動作原理にもとづく人口知能の実現を目指して命名されたものである。最近のコンピュータの発達には大きく分けて二つの流れがあって，一つは第 5 世代コンピュータに代表される人間の思考を模倣する人口知能コンピュータの開発であり，もう一つはスーパーコンピュータと呼ばれている超高速科学技術計算用コンピュータの開発である。前者は未だ発展途上で，専用の人口知能コンピュータが広く利用される状況には至っていない。

商用のスーパーコンピュータ (CDC6600) が誕生したのは 1964年のことであった。 CDC6600 の演算速度はおよそ 1MFLOPS (1 秒間に 100 万回の浮動小数点演算が可能) であり， ENIAC に比較して，乗算速度で 3,000 倍程であった。その後のスーパーコンピュータの性能向上は目覚ましく，最近では CDC6600 の 10,000 倍 ( 10 GFLOPS = 10×1024 MFLOPS) を越える演算速度をもつものも実用化されている。したがって，単純計算を行ってみると，最新のスーパーコンピュータは人間の 8, 400 × 3, 000 × 10, 000 ≈ 2.5 × 10^{11} 倍，つまり 2,500 億倍以上の演算能力をもっていることになり，最新のスーパーコンピュータが一秒間に処理する計算を一人の人間が行うとすると約一万年かかってしまうという計算になる。

世界最初の電子計算機となった ENIAC の操作方法は，必要な計算結果を得るために手動で結合線の組み替えを行うという複雑なもので，様々な問題に利用しようとすると不便なものであった。ところで，コンピュータの発達の歴史で忘れることのできない人物は，数学者でありコンピュータ開発の先駆者の一人であったジョン · フォン · ノイマン (〜1957 年) であろう。今日のコンピュータの多くが採用している方式はノイマン方式と呼ばれていて，ノイマン方式のコンピュータでは，計算の手順 (プログラム) を予め記憶装置 (メモリー) に電気的に記憶させておき (プログラム内蔵)，それを逐次実行させる (逐次処理) ので，容易にプログラムを変更することができるという特徴がある。ノイマン方式のコンピュータの出現によりコンピュータの利用範囲が爆発的に増大して，現在のコンピュータ時代が現出したともいえよう。ノイマンが目指していたのは，物理学で，気象現象など自然の法則を数学的に表現した場合によく現れる非線形偏微分方程式 (後述) を数値的に解くことであった。ノイマン方式のコンピュータの原型機は 1949〜50 年に英国と米国で相次いで開発されたが，ノイマンは志半ばにして病死し，結局，彼の存命中には彼が満足できるほどのコンピュータは完成されなかった。しかし，存命中に彼は強力な指導力を発揮し，彼の強固な意志が現在のコンピュータの基礎を確立させたといわれている。後に述べるように偏微分方程式を数値的に解くことは計算機物理学の大きな柱の一つでもあるので，計算機物理学の出現によってノイマンの夢はまさに花開いたといえるだろう。

3 計算機物理学の背景

3.1 モデル (物理模型)

\begin{aligned}-2 & -1.5 & -1 & -0.5 & 0 & 0.5 & 1 & 1.5 & 2 & 2.5 & -1 & 0 & 1 & 2 & 3 & 4 & k = 0.7 & r = 1/(1 − k cos θ) & 焦点 1 & 焦点 2 & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × & × \\ & 図 1: 楕円軌道の例 \end{aligned}自然現象は通常複雑な様相を示すので，

物理学ではそこから普遍的な性質を抽出して，様々な現象を統一的に説明するために，必要に応じたモデル (物理模型) を用いる。モデルといってもプラスチックモデルやファッションモデルなどのように形が目に見えるものではなく，考えている自然現象に本質的と思われる物理の基本法則を組み合わせて，多くの場合，数学を用いて表現したものである。導入したモデルを「解く」ことによって物理的な知見を得ることができ，仮定された物理法則の正当性が試されることになる。

例えば太陽系の惑星の運行を考えるとしよう。太陽系の惑星の運行については，ドイツの天文学者ケプラー(1571–1630)によって彼の師ティコ=ブラーエが生涯をかけて蓄積した膨大な量の火星軌道の観測データに基づいて発見された，いわゆる「ケプラーの法則」 :

1. 惑星は太陽を一方の焦点とする楕円軌道上を運行する。

2. 惑星の面積速度は一定である。

3. 惑星の公転周期の 2 乗は太陽からの平均距離の 3 乗に比例する。

がよく知られている。

惑星の運行に関する最も単純なモデルとしては，太陽と一つの惑星を質点と見なして両者が互いにニュートン (1642–1727) の万有引力の法則^{1} にしたがう力を及ぼしあって運動しているというモデルを採用することができる。この問題は「ケプラー問題」と呼ばれていて，ニュートンによって完成された (古典) 力学の枠内で，微分積分法 (後述) を用いた解析的な方法で，しかも厳密に解くことができ，経験的に発見されたケプラーの法則を理論的に「導く」ことができる。この問題については他の惑星の運行まで考慮すればより精密なモデルが得られ，この精密化したモデルを解くことによってより精密に惑星の運行が予言できる。つまり，ケプラーの法則を “基本法則” として採用しなくとも適当なモデルを用いて力学を適用すれば，観測される惑星の運行は説明できるわけである。

ところで，一般に物理モデルが数学の力を借りても解析的に厳密に解ける場合は非常に少ない。実際，３体以上の多体問題が一般的には解析的に解けないことが数学的に証明されている。上に述べた太陽系の惑星の運行の場合を例にとっても，太陽と一つの惑星の場合のような２体問題は厳密に解けるけれども，モデルを精密化してもう一つの惑星の影響まで考慮するとモデルは３体問題となってしまい，厳密には解くことができなくなってしまうのである。それでは，モデルを厳密に解くことは現象の物理的理解にとって不可欠であろうか? 答えは殆どの場合，「否」である。もちろん，モデルを厳密に解くことができれば，それによって「確実」な結論を導くことができ，それに越したことはないのだが，もともと物理学では理論とその実験的検証は不可分であり，また，実験的検証では測定誤差を完全に排除することが不可能なので，必要な精度でモデルの「近似解」が得られれば十分な場合が多いからである。また，精度は不十分であっても，モデルの定性的な傾向がわかれば仮定した物理法則の妥当性がチェックでき，有用である。

1 二つの質点には両者の質量に比例し，距離の二乗に反比例する引力がはたらく

3.2 摂動論

モデルが与えられたとき，その近似解を得るための方法として代表的なのは「摂動論」と呼ばれるものである。摂動論では，まず，解けるモデルを出発点として，そこに弱い撹乱としての「摂動」が加わった場合を想定する。例えば，惑星の運行の問題では，考えている惑星と太陽の二体問題が解けるモデルであり，他の惑星の影響を摂動として取り入れるのである。摂動の強さを表すパラメタを g としよう。つまり，摂動が弱い場合には g ≪ 1 である。摂動論とは，計算したい物理量 A が

A = a_{0} + ga_{1} + g^{2}a_{2} + ⋯  (1)

のように g についての「べき」の形に表すことができる場合に，係数 a_{0}, a_{1}, a_{2}, ⋯ を逐次，系統的に計算する手法である。このようにある量を小さな変数で展開して近似する手法は，コンピュータで関数の値を数値的に計算するときなどにもしばしば利用される，関数の「多項式近似」と類似している。初等関数の多くは

\begin{matrix}sin x & ≈ & x − \frac{1}{6}x^{3} + \frac{1}{120}x^{5} + ⋯ \\ e^{x} & ≈ & 1 + x + \frac{1}{2}x^{2} + ⋯ \\ log(1 + x) & ≈ & x − \frac{1}{2}x^{2} + \frac{1}{3}x^{3} + ⋯ \end{matrix} \}  (2)

のように無限に続く「べき展開^{2} 」の形に表すことができるが，多項式近似というのは， |x| ≪ 1 のときにその「無限級数」を有限項で打ち切って近似的な式を得ることである。正弦関数 sin x について項数を増すにつれて近似が良くなっていく様子を図 2 に示す。

摂動論による計算の歴史は古く， 19 世紀には天王星の軌道が精密に計算され，観測された軌道が計算で求められたものからずれていることがわかり，海王星の発見 (1846) へと導かれたのは有名である。海王星の発見の過程はまさに力学の理論体系と摂動法の組合わせによる勝利といえよう。ところで，海王星の発見の年代からもわかるように精密な摂動計算はコンピュータが発明されるよりはるか以前から盛んに行われていて，そのような数値計算に科学者らが費やした労力は並大抵のものではなかったことは容易に想像できる。実際，計算機科学の先駆者であるバベッジが計算機械の製作を始めたのは，当時の著名な天

2 {正確には，摂動論における (1) 式の展開は「漸近展開」と呼ばれるもので， (2) 式の展開とは区別しなく \atop てはならない。}

-1 -0.5 0 0.5 1 1.5 2 0  0.5  1  1.5  2  2.5  3  3.5  4 sin(x) x x − x^{3}/6 x − x^{3}/6 + x^{5}/120 図 2: sin 関数の多項式近似

文学者であるハーシェルと共に数値計算に用いる数表の校正を行った際にあまりに誤りの多いことに気付いて，人間のミスによる誤りが生じない，機械による自動計算を思い立ったためといわれている。

現在では，種々の軌道計算の方法は既に技術として確立されていて，各種天体現象の予報や，通常の人工衛星から，「ひまわり」や「ゆり」などの静止衛星，さらにボイジャーのように太陽系の外に飛出す人工彗星 (?) にいたるまでの打ち上げや軌道修正に必要な非常に正確な軌道計算が，コンピュータの利用によって可能となっている。

モデルの近似解を微小なパラメタによる展開の形で計算するという広い意味での摂動論は物理学や工学のあらゆる分野で現在でも活発に用いられて，とくに軌道計算の例のように，値が一定値に収束するまで高次の項の計算が可能な場合には，得られた結果は数値的には正確なものと結論でき，あいまいさのない議論が可能となる。

3.3 摂動論の破綻

摂動論はこれまでに様々な問題に適用されて，非常に多くの成果が得られているが，しかし，すべての問題が摂動論で解明され得るかというと，残念ながらそうではない点に注意しなければならない。

摂動論が適用できて，良い結果を与えるのは，摂動の強さを表す結合定数の値が小さい，「弱結合系」と呼ばれている体系のモデルにたいしてである。一方，考えている体系のモデルの結合定数が小さくない場合には，無理に摂動計算を行っても良い結果は得られない。結合定数が大きく摂動論が適用できない体系は「強結合系」と呼ばれている。弱結合系の性質は摂動論により既に殆どが調べられているのにたいして，強結合系の問題はその解明の困難さから未解明の問題として残されているものが多い。つまり，現在物理学で問題となっているものの大部分は強結合系であり，摂動論というほぼ確立された方法は，そのような興味ある問題には残念ながら単純には適用できないのである。

弱結合系であっても単純な摂動論が適用できない場合がある，という点にも注意すべきである。ここでの議論に関連する物理学史の上の最も重要な出来事は，「超伝導」現象の発見とその解明の過程であるので，少し詳しく説明しよう。

超伝導現象というのは，鉛などの金属の温度を低くしていくとある「臨界温度」で電気抵抗が突然に，しかも完全に消失するという現象で， 1911 年にオランダのカメリン · オンネスにより発見された。超伝導現象というのは，通常は原子，分子といったミクロな世界で現れる量子論的な効果がマクロなスケールで現れたものと理解されている。金属の超伝導現象の理論的解明は 1957 年にバーディーン - クーパー - シュリーファー (BCS) の 3 人によってなされたが，現象の発見から解明までに，実に半世紀近くを要しことは注目に値する。通常の金属の性質は，金属中で電気伝導を担う多数の伝導電子が互いに相互作用をしないで自由に動きまわるという自由電子モデルによりかなりよく説明できる。 BCS 理論は，金属中の伝導電子間に，格子振動を媒介とした「弱い」引力相互作用が働くことを前提としていて， g ≪ 1 ，すなわち弱結合系の理論である。この理論によってそれまでに知られていた金属の超伝導に関する殆どすべての実験事実を，半ば定量的に説明することができるようになり^{3} ，彼らは 1972 年のノーベル物理学賞に輝いたのである。 BCS 理論によれば電気抵抗が消失する温度，すなわち臨界温度 T_{c} は

T_{c} = 1.14 θ_{0} e^{−1/g}  (3)

と与えられ，ここに， θ_{0} は物質に固有なパラメータ， g は格子振動を媒介とした伝導電子間の相互作用の大きさを表す結合定数である。ところで， (3) 式の右辺は， g → 0 とすると指数関数の引数が発散してしまうので，たとえ g がどんなに小さくとも T_{c} が (1) 式のような g のべき展開の形で近似することができないことを表してる点に注目していただきたい。量子力学が完成した頃 (1926∼) には，摂動論による方法は方法論としてほぼ確立していて，量子力学はその後，多くの場合摂動論を用いて，様々な問題に適用され大成功をおさめていたのであるが， BCS 理論の結果は，結合定数 g がどんなに小さくても，超伝導状態を摂動論で取り扱うことができないことを意味し，このことが超伝導現象の理論的解明を遅らせたといえる。また，摂動論では決して解明できない現象が自然界に存在するということが認識された意義は非常に大きい。

3.4 計算機実験のアイデア

摂動論が適用できない物理系の性質はどのようにして調べれば良いのであろうか? 一つの解答は，対応するモデルの厳密解を解析的に求めることである。しかし，興味のあるモデルについて厳密解を解析的に求めることは一般的には非常に難しく，残念ながら，今日までに知られている厳密に解けるモデルのなかで，実在する体系のモデルとして相応しいものはほんの一握りのものに過ぎない。

3 最近の話では， 1986 年のベドノルツとミュラーによるいわゆる「酸化物高温超伝導体」の発見は記憶に 新しいところであるが，酸化物高温超伝導体は結晶構造が複雑で代表的な強結合系でもあり，未だその メカニズムは解明されていない。

摂動論が適用できない物理系の性質を調べるもう一つの方法は，何等かの近似方法を見出すことである。良い近似方法が見出せれば，それにより知りたい物理量が近似的に計算でき，体系の性質を定性的に，また，うまくいけば半定量的に説明できる。しかし，摂動論が弱結合系についてほとんど万能であったのにたいして，強結合系の場合には個々のモデルに応じた近似方法を考えなくてはならず，このことが強結合系の性質の解明を遅らせているといえる。また，弱結合の場合でも超伝導の例のように摂動論が適用できない場合があるし，問題によっては下手な近似をするとモデルがもつ物理的な本質を失わせてしまう場合もあり，このようなときは全くお手上げである。

そこで考え付くアイデアが，コンピュータを利用すればモデルを数値的に，明白な近似を用いないで調べることができるのでは，というものである。このアイデアの歴史も古く，コンピュータの父であり， 1945 年にプログラム内蔵方式のコンピュータの概念を確立したノイマンが夢みていたのも，実は，多くの自由度をもつ系の量子力学的問題を解くという物理の問題であった。

主に多自由度系の物理的性質を，コンピュータを用いたモデルの直接計算によって調べることを「計算機実験」または「計算機シミュレーション」と呼んでいる。コンピュータが誕生して間もない 1950 年頃には，既に，イタリア出身の理論 · 実験物理学者フェルミは，若い頃より自身が抱いていた非線形力学の問題を当時のコンピュータを用いて調べることを考えていて， 1955 年には実際にパスタとウラムと共同して計算機実験を行ったが，得られた結果は彼の予想に反したものであった。そして，このことが，現在では「ソリトン」と呼ばれている，非線形波動現象の発見の糸口となったのである。それ以降も，計算機実験の結果が常識的な予想を覆し，新しい現象の発見へと導かれたことが相次ぎ， 1980 年代頃より計算機物理学が理論物理学と実験物理学に並ぶ第３の物理学として認知されるようになったのである。

計算機実験の利点は以下のように要約されよう：

1. 適切なモデルを設定することにより，着目している現象を覆い隠してしまう様々な要因を排除できるので，物理的な本質を捉えやすい。

2. 温度や圧力など，パラメタを自由に調節できる。したがって，実現しようとすると膨大な経費を要するような極限条件下や，現実の系では実現不可能な条件での物理的知見を得ることができる。

3. 自由度が比較的小さな場合には正確な結果が得られるので，解析的に厳密解が得られない問題にたいしての，近似理論の検討や改良のための規準となり得る。

4 微分方程式

4.1 力学と常微分方程式

既に前の章でも述べられているように，物理学は自然現象の裏に隠れた基本原理や基本法則を探り，自然 (あるいは自然界) のより深い理解を得ることを目的としていて，そのための「共通語」としては数学が使われる。とくに，物理学の基本法則は「微分方程式」で表現される場合が多い。そもそも，これから説明する微分法というのは，あの有名なニュートンが自らが体系化した力学の理論を表現する際に考案した^{4} ものといわれていて，このことからも物理学と数学の関係の深さは理解されるだろう。「微分法」について理解するために，まず，我々が日常的に使っている「速さ」の概念を数学的^{5} に表現することにしよう。

速さとは一定の時間にどれだけ進むかを表す概念であるので，同じ時間内ならば長い距離を進んだ方が，また，同じ距離ならばかかった時間が短いほど「速さ」は大きい，つまり速いということになる。そこで，「平均の速さ」を

(進んだ距離)/(かかった時間)  (4)

と定義すれば，我々の日常感覚に合ったものとなることがわかる。実際我々は，新幹線の「こだま」が静岡から東京まで約 180 Km を 1 時間 30 分かけて走行したとすると，そのときの平均時速は 180/1.5 = 120 (Km/h) であると瞬時に答えることができるだろう。しかし，ここで定義した平均の速さは運転席のスピードメータに表示される速さとは明らかに異なっている点に注意してほしい。スピードメーターに表示されるのは時々刻々変化する「瞬間の」速さであり，ある有限の時間間隔の「平均の」速さではない。それでは「瞬間の」速さはどのように数学的に表せばよいのだろうか。

簡単のため，物体が直線上だけを移動している場合を考えて，任意の時刻 t における物体の位置を，直線上にとった目盛 (位置座標と呼ぶ) で表し，その値を x(t) とすると，時間間隔 Δt の間に物体が直線上を進んだ距離は位置座標の変化

Δx = x(t ＋Δt) − x(t)  (5)

の大きさで表される。そこで，時間間隔 Δt の間の物体の「平均の速度」を

\frac{Δx}{Δt}  (6)

で定義する。ここで，物体が直線上を正の向き (目盛が増加する方向) に進んだときには Δx は正で，負の方向に進んだときには Δx は負となることに着目すれば， Δx/Δt の正負が物体の進む向きを表すことに注意しよう。ところで，時間間隔 Δt を限りなく小さくすると，その間の位置の変化 Δx の大きさも限りなく小さくなるが，多くの場合， Δx/Δt が一定の極限値に近付く^{6} ことは想像できるだろう。この極限値のことを t の関数 x(t) の t における「微分係数」と呼んで， dx/dt ^{7} と書く。 dx/dt は t に依存する関数なので， x(t)

4 {微分法の発明については同時期に数学者のライプニッツとの先陣争いがあったのだが，ここでは触れな \atop い}

5 数学的と言っても，高校生以上の方には十分理解できる内容なので我慢していただきたい。

6 極限値が存在しない場合もあるが，ここでは扱わない。

7 微分係数にたいする記号 dx/dt はライプニッツからの伝統である。ニュートンは ẋ を用いた。

の「導関数」とも呼ばれる。関数 x(t) からその導関数 dx/dt を求めることを， “x を t で微分する” という。したがって，速度 v_{x}(t) は位置座標 x を時刻 t で微分することによって得られ，

v_{x}(t) = \frac{dx}{dt} ≡ lim_{Δt→0} \frac{x(t + Δt) − x(t)}{Δt}  (7)

となる。こうして定義された速度 v_{x}(t) の「大きさ」が時刻 t における瞬間の「速さ」である^{8} 。

同様に，時間間隔 Δt の間に物体の速度が

Δv_{x} = v_{x}(t ＋Δt) − v_{x}(t)  (8)

だけ変化したとすると，物体のその間の平均の「加速度」は

\frac{Δv_{x}}{Δt}  (9)

となり，加速度 a_{x}(t) は速度 v_{x}(t) を時刻 t で微分して

a_{x}(t) = \frac{dv_{x}}{dt} ≡ lim_{Δt→0} \frac{v_{x}(t + Δt) − v_{x}(t)}{Δt}  (10)

と書ける。なお， (7) 式のように，速度は位置の時刻に関する導関数なので，加速度は位置の時間に関する「二階の導関数」であり，

a_{x}(t) = \frac{d^{2}x}{dt^{2}}  (11)

とも書く。こうして，位置の時間変化率である速度と，速度の時間変化率である加速度が微分を用いて表現されたことになる。

ところで，ニュートンの運動の第 2 法則 (運動方程式) によれば，物体の加速度 (ベクトル) a⃗ はそのとき物体に加えられている力 (ベクトル) f⃗ に比例し，物体の質量 m に反比例する。つまり，式で表せば

ma⃗ = f⃗  (12)

となる。直線 (x–軸) 上の運動に限定すると， (12) 式は

ma_{x} = f_{x}  (13)

に置き換えてよい。ここに， a_{x} と f_{x} は物体の加速度の x–成分と物体にはたらく力の x– 成分である。時刻 t における物体の位置を x(t) で表したので，時刻 t における物体の加速度 a_{x}(t) は (11) で与えられる。また，一般に，物体にはたらく力 f_{x} は物体の位置 x と時刻 t に依存することに注意すると，物体にはたらく力は f_{x}(x(t)， t) と書ける。したがって，これらをまとめると，結局，直線上を運動する物体の運動方程式は

m\frac{d^{2}x(t)}{dt^{2}} = f_{x}(x(t), t)  (14)

8 {運動が直線の上に制限されていないときは，一般に「速度」をベクトルとして v⃗(t) で表し，「速さ」を \atop その大きさとして v(t) と表すのが普通である。}

となることがわかる。このように導関数が関係した方程式を一般に「微分方程式」と呼び，力学の目的の一つは，与えられた初期条件 (例えば， t = 0 における x と v_{x} の値) にたいして (14) 式のような微分方程式を解いて，任意の時刻における位置 x(t) を求めることである。なお，ここの例での x(t) のように，求めたい物理量が一変数関数である場合の微分方程式はとくに「常微分方程式」と呼ばれている。

4.2 微分方程式の数値解法 (オイラー法)

ある種の微分方程式については，解析解を得る具体的な方法が確立している。しかし，そのような場合でも，手順が膨大で実質的には人が手 (と紙と鉛筆) で解くことが不可能であったり，一般的には，解析解が得られないのが普通である。最近ではコンピュータの能力が向上したために，比較的簡単な微分方程式については解法の手順を予めコンピュータに与えておき，「数式処理システム」^{9} を用いることでコンピュータによっても解析解を得ることができるが，一般的には，コンピュータを利用した数値解に頼らざるを得ない。

微分方程式を数値的に解くとはどのようなことだろうか。また，どのようにして数値解を得れば良いのだろうか。

具体的な例として，地表からの高さ h の地点から，水平となす角 θ ，速さ v_{0} で投げられた質量 m の物体の運動について考えてみよう (図 3)。ここで，物体には一様な重力と，空気の粘性による抵抗力がはたらくとしよう。空気抵抗の力の大きさを見積もることは，物体の速さや大きさなどに依存する複雑な要素が絡んで，一般的には非常に難しい問題であるけれども，小石程度の物体が空気中を自由落下する場合については，抵抗力の大きさは概ね，物体の速さ v の二乗に比例することが知られている。抵抗力の向きは，常に速度 v⃗ と逆向きである。そこで，空気中を運動する物体のモデルとして，物体にはたらく力 f⃗ が次式で与えられるとしよう。

f⃗ = mg⃗ − γvv⃗  (15)

ここで，右辺第一項は一様な重力，第二項は大きさが速さの二乗に比例する空気抵抗の力で， g⃗ は重力加速度 (ベクトル)， γ は比例定数である。

図 3 のように x–座標と y–座標を設定すると^{10} ， (12) 式の運動方程式は両辺を m で割り， x–， y– の各成分に分けることで

\begin{matrix}a_{x}(t) & = & f_{x}(t)/m \\ a_{y}(t) & = & f_{y}(t)/m \end{matrix} \}  (16)

となる。 a_{x}, a_{y} は加速度ベクトル a⃗ の成分である。また， (15) 式の力 f⃗ の成分は，

\begin{matrix}f_{x}(t) & = & −γv(t)v_{x}(t) \\ f_{y}(t) & = & −mg − γv(t)v_{y}(t) \end{matrix} \}  (17)

9 一種の人工知能。プログラミング言語としては Mathematica や REDUCE が有名である。

10 平面内の運動なのでもう一つの座標 (z) は考える必要がない。

[figure: 0; 5; 10; 15; 20; 0; 10; 20; 30; 40; 50; θ; v⃗_{0}; r⃗; v⃗; g⃗; x; y; h; v⃗_{0} = 20(m/s); θ = π/5(rad); 軌跡; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×] 図 3: 放物体の軌跡

となる。ここに， g ≃ 9.8(m/s^{2}) は重力加速度の大きさ， v_{x}, v_{y} は速度ベクトル v⃗ の成分，また，速さ v は

v(t) = \sqrt{v_{x}(t)^{2} + v_{y}(t)^{2}}  (18)

である。 f_{x} と f_{y} のそれぞれは， v を通して v_{x} と v_{y} の両方に依存することに注意しよう。題意より，「初期条件」は，時刻 t = 0 のとき，

\begin{matrix}x(0) & = & 0, & v_{x}(0) & = & v_{0} cos θ \\ y(0) & = & h, & v_{y}(0) & = & v_{0} sin θ \end{matrix} \}  (19)

である。問題は，以上の条件より，任意の時刻 t における x(t) と y(t) を得ることである。空気抵抗が無い場合，すなわち， γ = 0 のときには，この問題はコンピュータを使わずとも解析的に簡単に解くことができて，その解は

\begin{matrix}x(t) & = & v_{0}t cos θ \\ y(t) & = & −\frac{1}{2}gt^{2} + v_{0}t sin θ + h \end{matrix} \}  (20)

であることが知られている。しかし， γ ≠ 0 の場合にはこの問題は簡単には解析的に解くことができない。そこで，コンピュータの出番となるわけである^{11} 。

この問題を数値的に解くには， 4.4.1 で説明した微分法の導入と逆の手順を辿ればよい。微分を定義している (7) 式に着目していただきたい。 (7) 式のように，微分係数 dx/dt は比 Δx/Δt の Δt → 0 の極限値として定義されているので，十分に小さな Δt をとれば，たとえ Δt の値が有限であっても

\frac{dx}{dt} ≃ \frac{Δx}{Δt}  (21)

11 {この問題は， θ = π/2 の場合には 1 次元運動になり，結果は少し複雑な形になるものの解析的に解くこ \atop とができる。挑戦してみよ。}

と近似でき，原理的には， Δt の値を小さくとる程近似は良くなることがわかる。そこで，適当な小さな値として Δt を定めて， 0, Δt, 2Δt, ⋯ といった離散的な時刻

t_{n} = nΔt,  (n = 0, 1, 2, ⋯)  (22)

についてだ̇け̇の位置や速度などを考えることにしよう。速度が位置の時間微分であり，加速度が速度の時間微分であることに注意して， (21) 式の近似を用いれば，速度と加速度の x–成分 v_{x}(t_{n}) と a_{x}(t_{n}) は

\begin{matrix}v_{x}(t_{n}) & ≃ & \{x(t_{n}) − x(t_{n−1})\}/Δt \\ a_{x}(t_{n}) & ≃ & \{v_{x}(t_{n}) − v_{x}(t_{n−1})\}/Δt \end{matrix} \}  (23)

と近似できる。 y–成分 a_{y}(t_{n}) と v_{y}(t_{n}) についても同様である。 (23) 式を用いて， (16) 式の運動方程式を整理すると，結局，解くべき近似的な方程式は

v_{x}(t_{n+1}) = v_{x}(t_{n}) + \frac{1}{m}f_{x}(t_{n})Δt  (24)
v_{y}(t_{n+1}) = v_{y}(t_{n}) + \frac{1}{m}f_{y}(t_{n})Δt  (25)
x(t_{n+1}) = x(t_{n}) + v_{x}(t_{n+1})Δt  (26)
y(t_{n+1}) = y(t_{n}) + v_{y}(t_{n+1})Δt  (27)

(28)

と書ける。このような方程式は「差分方程式」と呼ばれているが，微分方程式を差分方程式で近似する方法は一意的でない点に注意しよう。 (24) ∼ (27) で採用したのはその中でも最も簡単なもので，オイラー法と呼ばれているものである。

オイラー法によって常微分方程式の数値解を具体的に計算する手順 (アルゴリズム) についてまとめてみよう。

1. 与えられた初期条件により， n = 0 の場合について (24) 式と (25) 式の右辺を計算する。

2. 初期条件で与えられた x(0), y(0) と 1 の結果と (24)， (25) 式から定まる v_{x}(t_{1}), v_{y}(t_{1}) を (26)， (27) 式に代入して， x(t_{1}), y(t_{1}) を得る。

3. 初期条件の代りに，直前の段階で得られた x(t_{n}), y(t_{n}), v_{x}(t_{n}), v_{y}(t_{n}) を用いて 1.， 2. と同様の手順で， x(t_{n+1}), y(t_{n+1}), v_{x}(t_{n+1}), v_{y}(t_{n+1}) を得る。

4. 必要なだけ，手順 3. を繰り返す。

このアルゴリズムに従ったプログラムの例を章末に付録として付けておくので，興味のある人は参照して欲しい。

実際にプログラムを実行して，グラフにプロットした結果が図 4 に示されている。

[figure: 0; 5; 10; 15; 20; 0; 10; 20; 30; 40; 50; x; y; h; v_{0} = 20(m/s); θ = π/5(rad); γ = 0.00; ◇; ◇; ◇; ◇; ◇; ◇; ◇; ◇; ◇; ◇; ◇; ◇; ◇; ◇; 0.05; +; +; +; +; +; +; +; +; +; +; +; +; +; 0.10; □; □; □; □; □; □; □; □; □; □; □; □; □; 0.20; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; ×; 0.40; △; △; △; △; △; △; △; △; △; △; △; △; △; △; △] 図 4: 空気抵抗がある場合の放物体の軌跡

ここで述べた方法は直ちにケプラー問題にも適用することができるし，さらに自由度が多い場合についても，基本的には同様である。多くの粒子からなる体系にたいして，古典力学の運動法則に従って粒子を動かして，その体系の様々な性質を調べる計算機シミューレーションの手法を「分子動力学法」という。 1980 年代以後，スーパーコンピュータの発達と連動して，分子動力学法についての新しい方法論が次々と提案され，この方法を用いて調べることのできる現象や分野が大きくひろがってきた。分子動力学法の適用範囲は，「物性物理学」でミクロな原子や分子の集合体としての物質の性質を調べる場合から，「宇宙物理学」で銀河の集合としての超マクロな銀河団の形成の過程を調べる場合に至るまで非常に広い分野に及んでいる。

4.3 ‘場’ と偏微分方程式

第 1 章でも述べられているように，現代物理学で最も基本的な概念は「近接作用」と「場」 12 の考え方である。場というのは，任意の時刻および位置で，定義される物理量で，電磁気学における「電場」や「磁場」^{13} が代表例である。 4.4.1 では話を簡単にするために直線上の運動を考えたので，物体の位置を定めるのに一つ (一次元) の座標が必要であったが，一般には，位置を定めるのに縦 · 横 · 高さといった三つ (三次元) の「座標の成分」が必要となる。したがって，「電場」や「磁場」のような場の量は，一般に，時刻と位置座標の三つの成分を併せて，四つの変数に依存することがわかる。場というと難しく聞こえるかもしれないが，身近な例では，天気予報で示される気圧配置図を思い出してほしい。気圧配置図というのは，ある時刻に気圧の等しい地点を結んだ「等圧線」を示したものである。つまり，気圧配置図はある時刻の各地点の気圧を示しているのでこれは「気圧の場」を二次元的に表現したものといえる。

12 「界」と呼ぶこともある

13 「電界」や「磁界」とも呼ばれる

多変数関数 (変数の個数が二つ以上の関数) において，ただ一つの変数のみを変化させて，その変数について微分する，つまり，「偏導関数」を求めることを「偏微分」するという。偏導関数は (7) で定義した導関数の自然な拡張となっていて，例えば，「x と y の二変数関数 f(x, y) の x についての偏導関数」は

\frac{∂f(x, y)}{∂x} = lim_{Δx→0} \frac{f(x + Δx, y) − f(x, y)}{Δx}  (29)

と定義される。偏導関数に関係する方程式は「偏微分方程式」と呼ばれている。力学をはじめ，電磁気学，量子力学，流体力学など，物理学の基本法則の多くは偏微分方程式で表されているので，物理学上の問題の多くは，結局，偏微分方程式を解く問題に帰着する。

場の偏微分方程式を数値的に解くには，常微分方程式を数値的に解くのと基本的には同じアイデアを用いればよく， (21) 式を導入したときと同様に「微分」を「差分」で近似すればよい。力学の問題では位置や速度などの物理量が t のみの関数であったので，時間だけを離散的に考えたが，場の量は位置にも依存するので，場の偏微分方程式を差分方程式に近似するには，時間を離散的にとるだけでなく，空間についても格子に分ける必要がある。

場の偏微分方程式を数値的に解く問題として比較的身近なのは，天気の「数値予報」であろう。数値予報というのは，大気循環や気温変化などのモデルを「流体力学」や「熱力学」に基づいて導入し，観測値を初期条件にして未来の気象を予測するものである。数値予報に用いられるモデルは，基本的には流体力学や熱力学の基礎方程式である偏微分方程式をもとにしていて，日本付近の局所的な比較的詳しいデータをもとにした短期予報を目的にしたものから，全地球的な長期予報を目的としたものまでが開発されていて，更に改良が試みられている。日本国内の気象観測データは地域気象観測システム「アメダス」により収集されている。アメダスというのは，全国約 1,300 の地点に設置された無人自動観測所から毎正時に電話回線を通じて雨量，風向，風速，気圧，温度などの観測データを収集するシステムのことで，集められたデータは天気予報の基礎データなどとして利用されている。毎日の天気予報は，現在でも最終的には予報官の経験と勘に頼らざるを得ないようであるが，予報官は数値予報が出す客観的なデータを重要な参考としており，台風の進路予報など，近年の天気予報の精度の向上にたいする数値予報の役割は非常に大きいものといえる。

5 確率過程のシミュレーション

われわれの身の回りにある様々な物質が，非常に多く (10^{24} ∼ 10^{25} 個程度) の原子や分子から構成されていることを考えると，そのように多くの自由度をもった「巨視的な系」の様相はたいへんに複雑であろうと想像されるかもしれない。しかし，どこで水道の蛇口をひねっても出てくる水に大差はないように見えるし，百円硬貨を見ればどれも同じに見える。このことはどのように理解すればよいのだろうか。実は，自由度が大きい方が話が簡単になる場合があるのである。確率や統計に関係したものがそうである。例えば男女比について考えてみよう。われわれが日常関係している，家族，同級生，同僚など，比較的小さな集団ではその構成員の男女比はまちまちであろう。しかし，県内の全住民とか，日本人全体といったより大きな集団では，総人数に占める男性 (または女性) の割合は 1/2 に近づく。もともと性別は X–Y 染色体に関係して確率的に決まることが知られていて， 1/2 というのはちょうど男子 (または女子) が生まれる確率に等しい^{14} 。このように，大きな集団では平均値 (ここでは男性または女性の割合) が「確定」してしまうことを，確率論では「大数の法則」と呼んでいる。巨視的な系の性質が「温度」や「圧力」といった小数の巨視的な系に固有のパラメタで確定してしまうのは，まさに大数の法則が成り立つからである。物理学の分野で，巨視的な量の間の関係を調べるのは「熱力学」である。また，巨視的 (マクロ) な系の性質を，原子や分子といった微視的 (ミクロ) な性質から導こうとする分野を「統計力学」という。統計力学では巨視的体系の「微視的状態」が非常に多いことに着目して，その複雑さの中に確率的要素を取り込むのである。

巨視的な系といっても，それは原子や分子といった力学や量子力学にしたがう多数の微視的な系から構成されているのであるから，力学なり量子力学なりを正しく適用すれば正しい結果が得られるはずである。実際，分子動力学法によっても色々な巨視的な系の物理的な性質が精力的に調べられている。一方，統計力学の確率的要素を直接取り込んだ計算機シミュレーションの方法として「モンテカルロ法」^{15} も開発されている。モンテカルロ法は分子動力学法と並ぶ代表的な計算機シミュレーションの方法で，分子動力学法が運動方程式を直接解く決定論的な方法であるのにたいして，モンテカルロ法は (擬似) 乱数を用いた確率論的な方法である。

6 計算機物理学の守備範囲

計算機物理学の守備範囲は物理学の殆ど全ての分野にわたっているが，ここでは，最近の成果の中からいくつかを筆者の独断と偏見でとりあげ，スライドやビデオを交えて解説する。予定している内容は以下の通り (順不同):

1. 流れのシミュレーションとその可視化

2. 数値予報

3. 量子系の統計力学

4. 格子ゲージ理論と専用計算機

5. 非線形力学におけるソリトン

6. 宇宙物理学

14 実際には，さまざまな要因で男女比は正確には 1 : 1 にはならない。

15 {統計力学に限らず，乱数を用いた計算機シミュレーションの方法を総称して，「モンテカルロ法」と呼ぶ \atop 場合もある。}

7 最後に

これまで述べてきたように， 10 GFLOPS 級のスーパーコンピュータがかなり自由に使える時代を迎えて，計算機物理学は基礎から応用の時代へと入った。流体力学的計算や建造物の構造解析など，既に実用の域に達している分野もあるが，量子系の統計力学の場合のように，現在のコンピュータの性能ではまだ不十分な分野も存在する。今後も，コンピュータの能力は急速に進歩することが予想され，現在行われている計算はより短時間で処理されるようになり，応用の範囲が広がっていくことは疑いない。また，場合によっては現在全く予期されていない分野の出現も有り得るだろう。この講義を機に，現在急速に発展している計算機物理学の世界に興味をもつ方が一人でも多く現れれば望外の幸せである。

参考文献

1. 唐木幸比古, “スーパーコンピュータの現状と展望”, bit 臨時増刊「スーパーコンピュータと大型数値計算」 , 共立出版 (1987).

2. “計算機物理学 · 講習会テキスト”, 日本物理学会 (1989).

付録 a プログラミング例

図 4 を描くためのデータを出力する全く同じアルゴリズムを FORTRAN と BASIC の二種類の「言語」で実現した。

付録 a.1 FORTRAN の場合^{16}

PROGRAM shot
DIMENSION r(2),v(2),a(2)
CALL init(t,r,v,a,dt,ncalc,g,c)
CALL prtout(t,r,v,a,0)
CALL euler(t,r,v,a,0.0,1,g,c)
CALL prtout(t,r,v,a,1)
DO 10 i = 1,100
  CALL euler(t,r,v,a,dt,ncalc,g,c)
  CALL prtout(t,r,v,a,1)
  IF(r(2) .LE. 0.0) GOTO 999
10 CONTINUE
999 END
C
SUBROUTINE init(t,r,v,a,dt,ncalc,g,c)
DIMENSION r(2),v(2),a(2)
g = 9.81
t = 0
prper = 0.2
dt = 0.001
*      WRITE(*,*) ' プリントの間隔 (s), 時間の刻み巾 (s) = '
*      READ(*,*) prper,dt
ncalc = prper/dt + 0.5
height = 10
v0 = 20
th = 40
*      WRITE(*,*) ' 初期条件：高さ (m), 速さ (m/s), 角度 (degree) = '
*      READ(*,*) height,v0,th
th = th/180*3.141592
WRITE(*,*) ' 空気抵抗の係数 (1/m) = '
READ(*,*) c
v(1) = v0*cos(th)
v(2) = v0*sin(th)
r(1) = 0
r(2) = height
END
C
SUBROUTINE prtout(t,r,v,a,n)
DIMENSION r(2),v(2),a(2)
IF(n.EQ.0) THEN
  WRITE(*,'(/2(1X,2A/))')
+   '   t       x       y      v       a   ',
+   '    vx      vy      ax      ay  ',
+   '  (s)     (m)     (m)   (m/s)  (m/s^2)',
+   '  (m/s)   (m/s)  (m/s^2) (m/s^2)'
ELSE
  va = sqrt(v(1)**2+v(2)**2)
  aa = sqrt(a(1)**2+a(2)**2)
  WRITE(*,'(1X,F5.2,8F8.2)') t,r,va,aa,v,a
ENDIF
END
C
SUBROUTINE euler(t,r,v,a,dt,ncalc,g,c)
DIMENSION r(2),v(2),a(2)
DO 10 i = 1, ncalc
  r(1) = r(1) + v(1)*dt

16 一桁目の * を消去すると，各種パラメタをキーボードより入力できる。

        r(2) = r(2) + v(2)*dt
        vv = sqrt(v(1)**2 + v(2)**2)
        a(1) = - c*vv*v(1)
        a(2) = - g -c*vv*v(2)
        v(1) = v(1) + a(1)*dt
        v(2) = v(2) + a(2)*dt
        t = t + dt
        IF(r(2) .LE. 0.0) RETURN
   10 CONTINUE
      END
æ

付録 a.2 BASIC の場合^{17}

'PROGRAM shot
DIM r(2), v(2), a(2)
CALL Init(t, r(), v(), a(), dt, ncalc, g, c)
CALL Prtout(t, r(), v(), a(), 0)
CALL Euler(t, r(), v(), a(), 0, 1, g, c)
CALL Prtout(t, r(), v(), a(), 1)
LET i = 0
DO
  LET i = i + 1
  CALL Euler(t, r(), v(), a(), dt, ncalc, g, c)
  CALL Prtout(t, r(), v(), a(), 1)
LOOP UNTIL i > 100 OR r(2) < 0
END
SUB Euler (t, r(), v(), a(), dt, ncalc, g, c)
FOR i = 1 TO ncalc
  LET r(1) = r(1) + v(1) * dt
  LET r(2) = r(2) + v(2) * dt
  LET vv = SQR(v(1) ^2 + v(2) ^2)
  LET a(1) = -c * vv * v(1)
  LET a(2) = -g - c * vv * v(2)
  LET v(1) = v(1) + a(1) * dt
  LET v(2) = v(2) + a(2) * dt
  LET t = t + dt
IF r(2) <= 0 THEN EXIT SUB
NEXT i
END SUB
SUB Init (t, r(), v(), a(), dt, ncalc, g, c)
LET g = 9.81
LET t = 0
LET prper = .2
LET dt = .001
'INPUT "プリントの間隔 (s), 時間の刻み巾 (s) = "; prper,dt
LET ncalc = prper / dt
LET height = 10
LET v0 = 20
LET th = 40
'INPUT "初期条件：高さ (m), 速さ (m/s), 角度 (degree) = "; height,v0,th
LET th = th / 180 * 3.141592
INPUT "空気抵抗の係数 (1/m) = "; c
LET v(1) = v0 * COS(th)
LET v(2) = v0 * SIN(th)
LET r(1) = 0
LET r(2) = height
END SUB
SUB Prtout (t, r(), v(), a(), n)
IF n = 0 THEN
  PRINT

17 一行目を除いて，一桁目の ' を消去すると，各種パラメタをキーボードより入力できる。

         PRINT "   t       x       y      v       a   ";
         PRINT "    vx      vy      ax      ay  "
         PRINT "  (s)     (m)     (m)   (m/s)  (m/s^2)";
         PRINT "  (m/s)   (m/s)  (m/s^2) (m/s^2)"
         PRINT
       ELSE
         LET va = SQR(v(1) ^2 + v(2) ^2)
         LET aa = SQR(a(1) ^2 + a(2) ^2)
         PRINT USING "##.##"; t;
         PRINT USING " ####.##"; r(1); r(2); va; aa; v(1); v(2); a(1); a(2)
       END IF
       END SUB
       æ

付録 b 計算機実習

実習では実際に全員にパソコンに触れてもらって，計算機物理学の雰囲気を味わっていただきたい。内容的には高度なものを避けているので，自宅や学校，会社などでパソコンを利用できる方は，ここでの実習で用いたプログラムをもとにして自らプログラミングを行い，色々と試みると良いだろう。

課題: 以下の問題についてシミュレーションを行う前に考えてみよ。また，シミュレーションの結果と比較し，予想が外れた場合はその原因を考えよ。

付録 b.1 モンキーハンティング

高さが h で，水平方向に ℓ だけ離れた二つの地点のそれぞれから，二つの同等な物体を同時に，一方は静止状態からそっと離し，他方は水平方向へもう一方の物体に向かって速さ v_{0} で放った。

1. 空気抵抗が無視できるとき，どうなるか。

2. 空気抵抗の大きさが物体の速さの二乗に比例するとき，どうなるか。

3. 空気抵抗の大きさが物体の速さに比例するとき，どうなるか。

付録 b.2 ケプラー問題

1. 人工衛星が円軌道上を周回しているときに，軌道に沿って加速した。軌道はどのように変化するか。

2. 人工衛星が円軌道上を周回しているときに，軌道と垂直に外側に加速した。軌道はどのように変化するか。

3. 地球と木星の軌道データを用い，木星が地球の軌道に及ぼす影響について調べよ。

付録 b.3 電荷分布と静電界

二次元 (x, y) 平面上に分布した電荷がつくる電界について調べよ。なお， q(x, y) は位置 (x, y) に置かれた電荷 q を表す。

1. 1(1, 0), −1(−1, 0)

2. 1(1, 0), −4(−1, 0), 3(0, 1)

3. 1(1, 0), −1(−1, 0), 1(−1, −1), −1(1, −1)

付録 b.4 モンテカルロ法

1. [0, 1] の一様乱数を二つずつ生成して i 番目の組みを x_{i}, y_{i} とする。生成した組の総数を N とし， x_{i}^{2} + y_{i}^{2} < 1 (i = 1, 2, ⋯, N) を満たす組の数を M とする。 N を大きくすると， 4M/N はどのような値に近づくか。

2. [0, 1] の一様乱数を十個ずつ生成してそれらの和を計算する。 i 番目の組の和を t_{i} とすると， t_{i} は乱数としてどのような度数分布をするか。

参考文献

1. H. Gould and J. Tobochnik, “An Introduction to Computer Simulation Method Applications to Physical Systems”, Addison-Wesley (1988).