LINUX UNIX Tips Gnuplot 1

数値計算法・応用数理科学 講義ノート

北里大学 理学部 物理 第１章 はじめに 第２章 逐次計算 第３章 定積分 第４章 常微分方程式 第５章 行列計算 第６章 補追：偏微分方程式 付録１：LINUX（UNIX）お役立ちTips 付録２：Gnuplot簡易マニュアル なお、このテキストおよび例題のプログラムはネット経由で入手可能にしておくので、 適宜に利用してもよい。

1

はじめに

1.1

数値計算とアルゴリズム

「数値計算法」および「応用数理科学」は、おもに物理学の問題において必要となる数 値計算・数値解析を 精度良くかつ迅速に 実行するためのアルゴリズム（算法＝さんぽう、ともいう）と、それを用いて具体的な計 算をおこなうプログラム作成の実際について、講義と演習を行なうものである。 アルゴリズムという用語は９世紀アラビアの数学者アル・フワリズミの名前に由来する。 彼は『Algebr w’al muquabala』というアラビア語で書かれた最初の数学書を著し、２次 方程式の解法を論じたという。そこに何が書かれているかは、容易に相像できると思う。 代数学をalgebraとよぶのも、これに由来している。このように「解を（有限回の手続き で）求める方法＝手順」のことを一般にアルゴリズムという。 数値計算アルゴリズムの世界は、多岐に渡り奥が深いのであるが、ここで取り上げるも のは「初歩的かつ標準的なもの」ばかりなので、是非ともよく理解して自由に使いこなせ るようになってもらいたい。いろいろなアルゴリズムは、例えて言えば 大工さんの道具箱 のようなもので、手になじんだ鉋（かんな）や鋸（のこぎり）を適材適所に使えば、小さ な小物から大きな家屋までも作ることができるのである。 あらかじめ注意をひとつだけ述べておくと、どのアルゴリズムも「万能」ということは あり得ないので、適用条件をよく心得ておくことが大切である。そのためにも、得られた 結果をよく吟味して、「妥当な結果であるか？ もっと改良する余地はないか？」と、い つも考える習慣を身に付けてほしいものである。 「計算機が出した答えに間違いはない」というのは幻想に過ぎない のである。

1.2

プログラミング言語について

この演習では、言語として プログラム言語 C を使用する。過去には、Cに加えてJavaも用いてみたのだが、慣れないうちは一つだけ でも苦痛らしいので、今回は C のみにして、余裕があれば他の言語についても言及する ことにする。 どの言語の場合でも、ある程度まで使いこなせるようになるのは、それなりに大変であ る。人によっては、特定の言語に好き嫌いや得手不得手があったりもする。私のみるとこ

ろでは、これらの「好き嫌い・得手不得手」は、結局は慣れの問題であるように思う。初 めは少し苦労があるかもしれないが、いくつか例題を経験するうちに「こういう時は、こ うする」という技術が身に付くようになり、明るい眺望が開けるようになるはずなので、 しっかり取り組んでほしい。 以下に掲げる[例題]については、模範プログラムを示しているので、よく読んで充分理 解してもらいたい。その上で、[問題] のプログラムを例題に倣って書いてみてほしい。た いていは「例題＋ちょっとした工夫」で解決するはずである。プログラミングは、少しの 工夫次第で良くも悪くもなるところがおもしろいので、諸君もそういう工夫を心がけて大 いに演習を楽しんでもらいたい。 以下では C を使って数値計算し、結果の描画に Gnuplot を用いるのを基本型とする。 C言語を用いて、目的とする計算プログラムが自在に書けるようになれば、この演習の目 的は達せられたといえる。 数値計算では 1. 数式を用いた計算を行う 2. パラメータを入力したり、計算結果を出力する 3. 結果を理解しやすいよう画像表示する（可視化という） などができれば、おおよそ充分である。 これらの目的に限定すれば、それぞれの言語仕様のうち知っておくべきことはそれほど 多くない。物理の学生にとって 数値計算を通じてプログラム言語を学ぶ というのは、物理計算という「動機付け」が明確であるため、その言語を習得する良い近 道ともいえるであろう。 その上で将来、数値計算にとどまらず例えばDB（データベース）やWEBアプリケー ションの作成などに携わるために、これら言語にもっと習熟したいという諸君にも、この 演習で学んだことがさらに深く学ぶ足掛かりとなれば幸いである。

1.3

プログラミング上達法

（１）達人の書いたお手本を読む。 （２）自ら入力し、実行してみる。 （３）何でも試してみる。 （４）コンピュータの身になって考える。 （５）おおよその手順を決めてから書き始める。 （６）変数名はわかりやすく、注釈はこまめに。

1.4

最初のプログラム

最初のプログラムは、数値計算ではないが、定番の「Hello world!」を取り上げよう。こ れは コマンドライン（標準出力）に何かを出力する方法 の演習である。 [例題1-1]（Hello world!） 標準出力に「Hello world of C!」と出力するプログラムを書け。 [解答] これをC で書くと、以下のようになる。 /* hello.c */ #include <stdio.h> main(){ printf("Hello world of C!\n"); } １行目はプログラム名を記述した注釈行で、必ずしも書く必要はない。けれどもプログ ラム名の命名法は大切である。中身を見ないでも内容のよくわかる名前を付けるよう心が けるとよい。 念のため実行までの流れを書いておくと、以下の通りである。Cの場合、コンパイル→ 実行は $ gcc hello.c $ ./a.out とする（./ は「現在のディレクトリにある」という意味）。 [問題1-1] 円周率π = 3.141592· · · を適当な桁数まで出力するプログラムを書け。 ただし、円周率の計算法については、各自で工夫してみよ。円周率の計算アルゴリズム のひとつを次章で学ぶので（[問題2-2]）、ここでは出来るだけ簡単な方法として、各言語 に標準で用意されている数学ライブラリを使って出力する。ただし、数学関数のみならず 各種ライブラリをプログラムで利用するには、一定の「宣言＝手続き」が必要で、これら は書式を知っていないと書けない。書籍・ネットなどで、上記の数学ライブラリの使い方 に関する情報を調査してみよ（C言語の場合は既知のはず）。自分なりの情報検索の具体 的方法を経験・会得するのも本演習の目的のひとつである。

2

逐次計算

この章では 何らかの逐次（繰り返し）計算をして、結果の数値を標準出力に書き出す という問題を考える。

2.1

漸化式

この節の主テーマは「再帰（繰り返し）計算とループの使い方」で、プログラミングの 基本中の基本である。とくに「再帰（recursion）」はアルゴリズムの代表ともいうべく、い ろいろな場面で効果的に使われるものなので、しっかり理解してほしい。 [例題2-1]（フィボナッチ数Fn の計算） フィボナッチ数とは、漸化式 Fn+1 = Fn+ Fn−1, F0 = F1 = 1 で定義される数列であ る。実行時に数値 nを読み込んで、対応する Fn の値を出力するプログラムを書け。 [解答] この問題は「入力を受け取りその値に応じて結果を出力する」という形式のプロ グラムである。「関数サブルーチン」を使った解答例を以下に掲げる。もちろんメイン文 の中で直接計算しても良いのだが、このほうがすっきりしている。このように、作業をい くつかの独立したユニットに分けて記述する 分割して統治せよ というのが「プログラム作成のコツ」のひとつである。これは「自分の頭の中を整理する」 ことでもあるし、ユニットに分けることで「使いまわしが効く」という効用もある。 /* fibo.c */ #include <stdio.h> int F(int n){ int i, F0=1, F1=1, F2; if(n==0) return F0; else if(n==1) return F1; else for(i=0;i<(n-1);i++){ F2=F0+F1; F0=F1; F1=F2; } return F2; } main(){ int n, ans;

printf("n="); scanf("%d",&n); ans=F(n); printf("\nF(%d)=%d\n",n,ans); } [問題2-1] n = 0, 1, 2,· · · に対して、階乗n!を計算・出力するプログラムを書け。大き い n のとき、結果が正しいかどうかを確認し、正しい値を出力するにはどうしたら良い か考えてみよ。 [問題2-2] 円周率 π を以下の漸化式によって近似計算せよ（アルキメデスの方法）。直 径=1の円の外接正 n角形の周長を Pn、内接正n角形の周長を pn とすれば、漸化式 P2n= 2Pnpn Pn+ pn , (2.1) p2n= √ P2npn (2.2) が成り立つ。なぜなら、Pn = n· tan(π/n), pn = n· sin(π/n)（確かめよ）を右辺に代入 し、三角関数の倍角公式を用いれば、上記の漸化式が証明される。正三角形（k = 0）の 場合（P = 3√3, p = 3√3/2）から始めて、正 n = 3· 2k 角形の場合までを計算し、各々 の k, n, P, pを出力せよ。 [問題2-3] ２つの正数 a > b があるとき、初期値a0 = a, b0= b として数列an, bn を an+1 = an+ bn 2 , bn+1= √ anbn, (n = 0, 1, 2,· · · ) (2.3) によってつくると、これらは同じ極限値に収束する（a_∞= b_∞）。この極限値を算術幾何 平均といい agM(a, b)であらわす。agM(√2, 1) を計算せよ。

2.2

求解法のいろいろ

この節では 与えられた方程式 f (x) = 0 の解を求める 問題を考える。こうした目的のためのアルゴリズムの代表的なものに、ニュートン法と２ 分法がある。 （１）ニュートン法 方程式 f (x) = 0の近似解を求めるのに、関数 f (x)の微分を用いて、適当な初期値 x0 から逐次的に xn+1= xn− f (xn) f0(xn) (2.4)

を繰り返すと、極限値x_∞が求める解となる。これにより、欲しい精度が得られるまで有 限回の再帰計算をおこなう方法を「ニュートン法」という。 [例題2-2]（平方根の計算１） 問題の関数を f (x) = x2− aとすれば、上記の漸化式は xn+1= 1 2 ( xn+ a xn ) (2.5) となる。与えられた a > 0に対して、その平方根を求めるプログラムを書け。また、繰り 返しの回数と近似精度の関係を調べ、結果（(xn− √ a) vs n）を Gnuplotで描いてみよ。 [解答] a = 2に固定した場合のCプログラムは以下のようになる。aの値を入力させる ように scanf を用いて書き直してみよ。また、ニュートン法は指数関数的に収束する |xn− x∞| ∝ e−c·n, (c =正定数) (2.6) ことが知られている。 /* newton.c */ #include <stdio.h>

double rhs(double x, double a){ return (x+a/x)/2;

}

main(){ int n;

double a, x, ans, error;

a=2; ans=sqrt(a); x=a;//初期値 for(n=0; n<10; n++){ x=rhs(x,a); error=x-ans; printf("x(%d)= %14.12f, error=%14.12f\n",n+1,x,error); } } （２）２分法 方程式f (x) = 0を解く別法を考える。中間値の定理によれば、連続関数f (x)において a < b のときf (a) < 0, f (b) > 0 となっていれば、少なくとも一つの解xがa < x < bに 存在する。このことを使って解を以下のように求める方法を「２分法」という。初期値を

a0 = a, b0 = bとして、数列an, bn を次のように構成する。すなわち中点cn= (an+ bn)/2 における関数の正負を調べて f (cn) < 0のときは an+1= cn, bn+1= bn とし、 f (cn) > 0のときは an+1= an, bn+1= cn とする。 以下これを繰り返せば、回を追うごとに区間縮小されて f (x)の零点すなわち解に近づく のである。もちろん途中で f (cn) = 0 となれば、それが解となる。 [例題2-3]（平方根の計算２） 同じく f (x) = x2− aとして、平方根を２分法で計算するプログラムを書け。繰り返し の回数と近似精度の関係を調べ、ニュートン法の場合と比較してみよ。 [解答] aの値によって、初期値の取り方に注意が必要である。a > 1の場合のプログラ ム例を挙げる。 /* bisection.c */ #include <stdio.h>

double f(double x, double a){ return (x*x-a);

}

main(){ int n;

double a, left, right, center;

a=2; left=0;

right=a;// a>1 のとき \sqrt{a}<a である for(n=0;n<10;n++){ center=(left+right)/2; if(f(center,a)<0) left=center; else right=center; } printf("x=%lf\n",center); } [問題2-4]（デカルトの正葉線） 次の方程式によって定義される曲線を描け。 x3− 3xy + y3= 0 (2.7) 注意：この曲線は直線 y = x に関して対称である。また、区間0 < x < 41/3 において多 価（３価）となっている。なお、この曲線は「デカルトの正葉線」と呼ばれている。

[問題2-5]（プランク分布） 熱輻射に関するプランクの分布則は P (x) = Ax 3 ex_{− 1}, (x = hν kBT , A =定数) (2.8) で与えられる。P (x)のx 微分を実行すると P0(x) = Ax 2 (ex_{− 1)}2 (3(e x_{− 1) − xe}x_{) (2.9)} であるから、分布の最大値を与えるxm > 0は方程式f (x)≡ (3 − x)ex− 3 = 0を満たす。 これを数値的に解いて xm の値を求めよ。（ヒント）xm= 2.8· · ·。 この結果から、例えば「太陽（T = 6000 K）から最もたくさん放射される光の振動数 （波長）」や「宇宙背景輻射（T = 3 K）が出す主要な電磁波の振動数（波長）」を求めて みよ。 [問題2-6] 強磁性体の磁化 m は、イジング模型の平均場近似によれば、m に対する自

己無撞着方程式（self consistent equation）

m = tanh(Km), (K = zJ/kBT ) (2.10) を解いて求められる。ここで z は最近接スピン数、J は交換相互作用エネルギーである。 これから臨界温度はK = 1のとき、すなわちTc = zJ/kB となることがわかる。パラメー タを Tc = zJ/kB = 1 に選び、区間0 < T < 1を N 分割し、温度T とそのときの磁化 m の値を出力するプログラムを書け。磁化m vs温度 T のグラフ（m(0) = 1, m(1) = 0） を、Gnuplotを用いて描いてみよ。

3

定積分

与えられた関数 f (x)を、区間a≤ x ≤ b について定積分することを考える。定積分値 を精度良く数値計算するには、関数f (x) の振る舞いや、端点a, b の取り方に応じて工夫 を必要とする場合がある。

3.1

台形公式、シンプソン公式

まず定積分の数値計算法の代表的な２つ、台形公式とシンプソン公式から始めよう。 （１）台形公式 区間 [a, b] をN 等分して、h = (b− a)/N, xn= a + nh（n = 0, 1, 2,· · · , N）とすれば （x0 = a, xN = b）、区分求積法によって定積分は S = ∫ b a f (x)dx = lim N→∞SN, SN = N_∑−1 n=0 h 2 · (f(xn) + f (xn+1)) (3.1) と書かれる。この極限をとる代わりに、有限の N での値 SN で近似したものを台形公式 という。台形公式は 各微小区間で、関数を直線（１次関数）で近似する ことに相当しており、計算誤差は h2 _{の程度であることが知られている。} （２）シンプソン公式 台形公式は「微小区間を直線で近似する」ものであったが、対してシンプソン公式は 隣り合う２つの微小区間で、関数を放物線（２次関数）で近似する ことによって、近似の精度を上げる公式である。３点 (xn−1, f (xn−1)), (xn, f (xn)), (xn+1, f (xn+1))を通る関数を２次関数で近似すると f (x) = f (xn−1) (x− xn)(x− xn+1) (xn−1− xn)(xn−1− xn+1) + f (xn) (x− xn−1)(x− xn+1) (xn− xn−1)(xn− xn+1) + f (xn+1) (x− xn−1)(x− xn) (xn+1− xn−1)(xn+1− xn) (3.2) と表されるので、これを区間xn−1≤ x ≤ xx+1で積分する。簡単のため等間隔xn−xn−1= xn+1− xn= hの場合で考えれば、積分結果は ∫ xn+1 xn−1 f (x)dx = h 3(f (xn−1) + 4f (xn) + f (xn+1)) (3.3)

となる。よって、全区間a≤ x ≤ bを偶数のN 等分に分割して、h = (b−a)/N, xn= a+nh とすれば S = ∫ b a f (x)dx = lim N→∞SN, (3.4) SN = h 3(f (x0) + 4 (f (x1) +· · · + f(xN−1)) + 2 (f (x2) +· · · + f(xN−2)) + f (xN)) が得られる。シンプソン公式の計算誤差は h4 の程度であることが知られている。 [例題3-1]（アークタンジェント） 被積分関数を f (x) = 1/(1 + x2) として、区間0≤ x ≤ 1で数値積分するプログラムを 書け。台形公式とシンプソン公式の両方で計算して、結果を較べてみよ。ちなみに答えは tan−1(1) = π/4となるはずである。 [解答] 以下の通り。ちなみに trapezoid（トラペゾイド）とは台形のこと。 /* sekibun.c */ #include <stdio.h> #include <math.h> double f(double x){ return (1/(1+x*x)); } main(){ int n, N=100; double h, x, S1, S2; h=1/(double)N; /* trapezoid */ S1=0; for(n=0; n<N; n++){ x=h*(double)n; S1+=(h/2)*(f(x)+f(x+h)); } printf("S1=%f\n", S1); /* simpson */ S2=0; for(n=1; n<N; n+=2){ x=h*(double)n; S2+=4*f(x); }

for(n=2; n<N; n+=2){ x=h*(double)n; S2+=2*f(x); } S2+=f(0)+f(1); S2=(h/3)*S2; printf("S2=%f\n", S2); /* exact */ printf("S =%f\n", atan(1)); } [問題3-1]（計算精度） [例題]のプログラムを分割数N が可変になるように書き換え、計算精度(SN− π/4) vs N を調べて、結果をGnuplotで描いてみよ。 [問題3-2]（楕円積分） 第１種完全楕円積分（0≤ k < 1）は K(k) = ∫ π/2 0 dθ √ 1− k2_sin2_θ (3.5) で定義される。これを数値積分して、そのグラフ（横軸 k、縦軸K(k)）を描け。 [問題3-3]（ベッセル関数のゼロ点） ０次のベッセル関数J0(x) は定積分 J0(x) = 1 2π ∫ 2π 0 cos(x sin θ)dθ (x≥ 0) (3.6) で定義される。この関数のゼロ点 0 < x1 < x2 <· · · のうち、最初の２つを数値的に求め よ。（ヒント）2 < x1 < 3, 5 < x2< 6 であることを用い、２分法によって計算する。

3.2

特異積分、広義積分

積分区間の端で被積分関数が無限大になるとか、積分区間が無限区間となる場合を特異 積分・広義積分という。こういう場合の数値積分に有効な方法として、「２重指数関数の 方法」（Double Exponential Method略してDE）がある。例を取り上げて、具体的に説 明しよう。 定積分 ∫ 1 0 dx √ 1− x2 = π 2 (3.7) は、左辺の不定積分が sin−1x であることを知っていれば計算するまでもなく答えがわか るのであるが、これを数値計算する場合には積分の右端（x = 1）で被積分関数が発散す

るので、何らかの工夫が必要である。試みに「端点を除いて」台形公式・シンプソン公式 を使い、その計算精度を調べてみるとよい。 このような積分の端点での特異性を「置換積分（変数変換）によって無限遠に押しやろ う」というのが、以下で紹介する計算法の要点である。変数変換 x = tanh (_π 2sinht ) (3.8) によって、区間 0 ≤ x ≤ 1は区間 0 ≤ t < ∞ に移される。これにより端点 x = 1 近傍 は「無限に長い領域」に引き伸ばされ、またtに関して等間隔な分割は（端点に近づくに 従って）より細かな区間分割に対応していることもわかる。こうして、特異積分を「非等 間隔区分求積法」によって精度よく近似する方法を「２重指数関数の方法」といい、高橋 秀俊・森正武によって考案された。 [例題3-2]（アークサイン） 例として挙げた定積分 S = ∫ 1 0 dx √ 1− x2 (3.9) を２重指数関数法で計算するプログラムを書け。変数変換(3.8) のヤコビアンは dx dt = π 2 cosht cosh2(π₂sinht) (3.10) で与えられることに注意し、適当な T に対して区間 0≤ t ≤ T で台形公式を適用せよ。 [解答] あらかじめ被積分関数を、変数変換した結果の形 S = ∫ _∞ 0 π 2 cosht cosh(π₂sinht)dt (3.11) に書いておくと良い。 /* de.c */ #include <stdio.h> #include <math.h> #define pi (4.0*atan(1.0)) double f(double t){ return cosh(t)/cosh((pi/2)*sinh(t)); } main(){ int n, N; double dt, t, S; N=500;

dt=0.01; /* arcsin(1.0) */ S=0; for(n=0; n<N; n++){ t=dt*(double)n; S+=(pi/2)*(dt/2)*(f(t)+f(t+dt)); } printf("S=%lf\n", S); /* exact */ printf("exact=%lf\n", pi/2); } [問題3-4] 次の定積分を数値計算せよ。これをガウスのレムニスケート積分という。 ω = ∫ 1 0 dx √ 1− x4 (3.12) 結果は ω = 1.3110287· · · となるであろう。 (参考) ガウスはこれが算術幾何平均 agM(√2, 1) （[問題2-3]参照）と agM(√2, 1) = π 2ω (3.13) の関係にあることを、数値計算を通して発見したという（高木貞治『近世数学史談』）。 [問題3-5] 次の定積分（ガウス積分と呼ばれる）を数値計算で確かめよ。 ∫ _∞ −∞e −x2 dx =√π (3.14) この場合の変数変換としては x = sinh (_π 2sinht ) (3.15) を用いて、区間 −T ≤ t ≤ T に対して台形公式を適用するとよい。 （注）微分を実行して dx dt = π 2cosht· cosh (_π 2sinht ) (3.16) を用いた置換積分を数値的に計算する。 [問題3-6]（エアリー関数） 微分方程式 d2Φ dx2 = x· Φ, Φ(+∞) = 0, Φ 0_{(+∞) = 0 (3.17)}

を満たす関数 Φ(x)をエアリー関数といい、光学や量子力学において重要な特殊関数のひ とつである。この関数のグラフは x > 0 で減衰、x < 0 で振動する： Φ(x)∼ 1 2x −1/4_exp(₋2 3x 3/2 ) (x > 0) (3.18) ∼ |x|−1/4sin ( 2 3|x| 3/2₊π 4 ) (x < 0) (3.19) このような振る舞いの概略を、数値計算によって確かめてみよう。それには、エアリー関 数の積分表示 Φ(x) = √1 π ∫ _∞ 0 cos ( xu + u 3 3 ) du (3.20) を用いる（確かめよ）。２重指数関数の変換 u = exp(π· sinh t), (−T1≤ t ≤ T2) (3.21) によって積分を実行し、関数 Φ(x)のグラフを描いてみよ。

3.3

応用 ― 超伝導の BCS 理論

超伝導のBCS（Bardeen-Cooper-Schrieffer）理論では、ギャップ方程式とよばれる積分 方程式 1 = V0NF ∫ ~ωD 0 dξ Etanh ( E 2kBT ) , (E =√ξ2_{+ ∆}2_{) (3.22)} を解いて、臨界温度やエネルギーギャップが計算される。 いま、エネルギーをすべてデバイ・エネルギー~ωD を単位に用いて無次元化し、 ξ ~ωD → x, E ~ωD → E, ∆ ~ωD → ∆, kBT ~ωD → T (3.23) と改めて表わすことにすれば、（無次元の）結合定数を V0NF = g と書いて、ギャップ方 程式は 1 = g ∫ 1 0 dx E tanh ( E 2T ) , (E =√x2_{+ ∆}2_{) (3.24)} となる。以下、これを解くことを考えよう。 [例題3-3]（臨界温度） 臨界温度 Tc は、そこでギャップ ∆ = 0となるので 1 = g ∫ 1 0 dx x tanh ( x 2Tc ) (3.25) を解いて求められる。Tc をg の関数（0 < g < 1）として計算するプログラムを書け。ま た Gnuplotを用いてTc をg の関数としてグラフに描け。 [解答] 以下の通り。

/* bcs.c */ #include <stdio.h> #include <math.h> double sekibun(double T){ double x, dx; int j, N=100; double f[N+1], bcs; dx=1/((double)N); f[0]=1/(2*T); for(j=1; j<=N; j++){ x=dx*(double)j; f[j]=tanh(x/(2*T))/x; } bcs=0; for(j=0; j<N; j++){ bcs+=(f[j]+f[j+1])*dx/2; } return bcs; } main(){ int j, k, max=50; double g;

double left, half, right; double residue; /* coupling constant */ g=0.1; for(j=0; j<10;j++){ left=0; right=1; for(k=0; k<max; k++){ half=(left+right)/2; residue=1-g*sekibun(half); if(residue>0){

right=half; }else{ left=half; } } printf("g=%3.2f Tc=%8.5f\n",g, half); g+=0.1; } } [問題3-7] 絶対零度 T = 0におけるエネルギーギャップ∆(0)を、結合定数gの関数と して、区間 0 < g < 1において計算し、そのグラフを描け。 [問題3-8] 結合定数 g を適当な値（例えば g = 0.2）に固定して、エネルギーギャップ ∆(T ) を温度（0 < T < Tc）の関数として計算し、そのグラフを描け。 [問題3-9] 絶対零度におけるエネルギー・ギャップ ∆(0)と臨界温度Tc との比は、結合 定数 g が小さいとき（BCS の弱結合極限） ∆(0) Tc = 1.764· · · (3.26) という一定値になるという。このことを数値的に確かめよ。

4

常微分方程式

4.1

微分方程式の差分近似

１階の微分方程式 dx dt = f (x, t) (4.1) を、初期条件 x(0) から出発して数値的に解く問題を考える。 左辺の微分は dx dt = limh→0 x(t + h)− x(t) h (4.2) で定義されるから、十分小さい h を用いれば 微分を差分で置き換えてもよい であろう。このことから、t = nh, xn= x(nh)と書いて（x のn乗ではない！） xn+1= xn+ h· f(xn, nh) (4.3) という微分方程式の差分近似方程式が得られる。これは初期値 x(0) = x0 _{から逐次的} （iterative）かつ陽（explicit）に計算するスキームで「オイラー法」とよばれている。こ の計算法は簡単ではあるが、精度が悪く（オーダーh2 の誤差がある：これを１次精度と いう）、また問題によっては誤差の累積によって発散する結果を与えるなど、問題点が多 いことでも知られている。 そこで普通は、次のような２次精度の計算法がよく使われている。 1. 修正オイラー法 k = f (xn, nh), xn+1= xn+ h· f(xn+h 2k, nh + h 2) (4.4) 2. ホイン（Heun）法 k1 = f (xn, nh), k2 = f (xn+ hk1, nh + h), xn+1= xn+ h 2(k1+ k2) (4.5) [問題4-1]（線形常微分方程式） 次の微分方程式をいろいろな計算法で数値的に解き、計算スキームによる精度の違いを 調べてみよ。 (1) ˙x =−x, x(0) = 1, (2) ¨x =−x, x(0) = 1, x(0) = 0˙ (4.6)

4.2

２階の常微分方程式

力学の問題によく登場する２階の常微分方程式 d2x dt2 = f (x, t) (4.7) の場合には、dx/dt = pと置き、１階の連立微分方程式 dx dt = p, dp dt = f (x, t) (4.8) に変形して考えればよい。この場合によく使われる計算法を挙げておこう。 （１）リープ・フロッグ法 最初は、式(4.7)の右辺がtに依らない場合（f (x)）に使われるもので、リープ・フロッ グ（Leap Frog かえる跳び）法とよばれる。これは２次精度の計算法である。 pn+1/2= pn+ h 2f (x n_{), x}n+1_{= x}n_{+ hp}n+1/2_{, p}n+1 _{= p}n₊h 2f (x n+1_{) (4.9)} なお、これを変形して４次精度にした計算法として 「４次の吉田法」があり、以下の [例題4-1] の解答で紹介する。 （２）ベルレ法 これも f がt に依らない場合の解法だが、p を経由しないで２階のまま xn+1= 2xn− xn−1+ h2f (xn) (4.10) によって計算する。直前の値 xn だけでなく、もうひとつ古い情報 xn−1 も記憶しておか ねばならない。 リープ・フロッグ法とベルレ法は、ともに時間反転（h→ −h）に関して対称な計算法 であることは興味深い。（初期でなく）末期条件から始めて、時間を遡る場合にも、同じ 計算法が使えるのである。 （３）ルンゲ・クッタ法 一方で、式(4.7)の右辺が時間に依存する場合（f (x, t)）は、これらの計算法が使えな い。その場合にも有効な計算法として、古くから使われている方法が「ルンゲ・クッタ （Runge-Kutta）法」で、４次精度の計算法である。 この方法は(xn, pn)から(xn+1, pn+1)を求めるのに、つぎの４段階の手続きを経る。す なわち x1= hpn, p1= hf (xn, nh), (4.11) x2= h(pn+ p1 2), p2 = hf (x n₊x1 2 , nh + h 2), (4.12) x3= h(pn+ p2 2), p3 = hf (x n₊x2 2 , nh + h 2), (4.13) x4= h(pn+ p3), p4= hf (xn+ x3, nh + h); (4.14)

と計算しておいて、最後に xn+1= xn+1 6(x1+ 2x2+ 2x3+ x4); (4.15) pn+1= pn+ 1 6(p1+ 2p2+ 2p3+ p4); (4.16) と、平均された増分を加える。これを繰り返して、時間発展を計算するのである。 [例題4-1]（非線形振動）非線形バネの運動方程式 dx dt = p (4.17) dp dt =−ω 2_{x− βx}3 _(4.18) を、パラメータ ω2, β および初期条件を適当に選んで数値的に解け。 [解答] これは正準力学系であるので、４次吉田法という解法が特に優れている。これは リープ・フロッグ法を（時間間隔を変えて）３回繰り返して、４次精度で計算する方法で ある。その時間間隔は h1 = h 2− 21/3, h2=−h 21/3 2− 21/3, h3= h1 (4.19) で与えられ、行きつ戻りつしながら（h1 = h3 > 0, h2 < 0, h1+ h2+ h3 = h）解いて いく。 以下に掲げたプログラムは３つの計算法（lf2, ys4, rk4）からひとつを選んで実行し、結 果をファイル（ *.dat ）に書き出す。それを Gnuplotに渡して位相空間における軌道を 描いてみよ。また、エネルギー保存の精度を比較してみよ。 /* osc.c */ #include <stdio.h> #include <math.h> #define dt (0.01) #define omega2 (1.0) #define beta (0.4) #define f(x) (-omega2*x-beta*pow(x,3))

void lf2(double x,double p,double *xn,double *pn,double h){ double ph;

ph=p+(h/2)*f(x); *xn=x+h*ph;

*pn=ph+(h/2)*f(*xn); }

void rk4(double x,double p,double *xn,double *pn, double h){ double x1,x2,x3,x4,p1,p2,p3,p4; x1=h*p; p1=h*f(x); x2=h*(p+p1/2); p2=h*f(x+x1/2); x3=h*(p+p2/2); p3=h*f(x+x2/2); x4=h*(p+p3); p4=h*f(x+x3); *xn=x+(x1+2*x2+2*x3+x4)/6; *pn=p+(p1+2*p2+2*p3+p4)/6; } main(){ double x,p,xn,pn; double d,d1,d2; double energy; int solver; int n,max=800; FILE *fp;

/* 4th order yoshida parameters */ d=pow(2,1.0/3.0);

d1=1/(2-d); d2=-d/(2-d);

printf("Which algorithm? (1:Runge-Kutta, 2:Leap Frog, 3:Yoshida)"); scanf("%d",&solver);

/* initial condition */ x=2.0;

p=0.0;

energy=p*p/2+omega2*x*x/2+beta*pow(x,4)/4; printf("First energy=%10.6f \n",energy);

/* file open */ switch(solver){ case 1:

fp=fopen("rk4.dat","w"); break; case 2: fp=fopen("lf2.dat","w"); break; case 3: fp=fopen("ys4.dat","w"); break; } fprintf(fp,"%10.7f %10.7f\n",x,p); /* time development */ for(n=1; n<max; n++){ switch(solver){ case 1:/* Runge-Kutta */ rk4(x,p,&xn,&pn,dt); x=xn;p=pn; break;

case 2:/* Leap Flog */ lf2(x,p,&xn,&pn,dt); x=xn;p=pn; break; case 3:/* Yoshida */ lf2(x,p,&xn,&pn,dt*d1); x=xn;p=pn; lf2(x,p,&xn,&pn,dt*d2); x=xn;p=pn; lf2(x,p,&xn,&pn,dt*d1); x=xn;p=pn; break; } /* write data */ fprintf(fp,"%10.7f %10.7f\n",x,p);

}/* end of time development */

/* check of energy conservation */

energy=p*p/2+omega2*x*x/2+beta*pow(x,4)/4; printf("Final energy=%10.6f \n",energy); }

[問題4-2]（ファン・デル・ポール方程式） 真空管の発振現象を記述する方程式として提案された d2x dt2 =−x + ε(1 − x 2₎dx dt (4.20) をファン・デル・ポール方程式という。右辺第２項の x2 の項がなければ線形な減衰振動 方程式となる。ここでは ε > 0なので、x2 < 1 のとき増幅、x2> 1 のとき減衰の効果を もたらす。その結果、x(t) は有限の振幅で振動するようになる。これが発振現象である。 ２階の方程式なので、p = dx/dtとして、連立方程式 dx dt = p, dp dt =−x + ε(1 − x 2_{)p (4.21)} に書き換え、ルンゲ・クッタ法を適用して、これを解け。ただし ε = 0.3とし、初期条件 x(0), p(0) は適当に設定せよ。 [問題4-3]（２重振り子） ２重振り子はたいへん不規則な運動をすることで知られている。そのラグランジアンは、 鉛直下方からの角度を θ1, θ2 として L =m1+ m2 2 ` 2 1θ˙12+ m2 2 ` 2 2θ˙22+ m2`1`2θ˙1θ˙2cos(θ1− θ2) + (m1+ m2)g`1cos θ1+ m2g`2cos θ2 (4.22) で与えられるから、運動方程式として ˙ θ1= m2`22p1− m2`1`2p2cos(θ1− θ2) m2`2₁`2<sub>2</sub>[m1+ m2sin2(θ1− θ2)] , (4.23) ˙ θ2= (m1+ m2)`21p2− m2`1`2p1cos(θ1− θ2) m2`21`22[m1+ m2sin2(θ1− θ2)] , (4.24) ˙ p1=−(m1+ m2)g`1sin θ1− m2`1`2θ˙1θ˙2sin(θ1− θ2), (4.25) ˙ p2=−m2g`2sin θ2+ m2`1`2θ˙1θ˙2sin(θ1− θ2) (4.26) を得る。ここで、第３、４式の右辺にある θ˙1θ˙2 は第１、２式を代入するものとする。パ ラメータ m1, m2, `1, `2, g の値および初期条件を適当に（例えば m1 = m2 = 1, `1 = 1, `2 = 0.8, g = 9.8など）定めて、θ1, θ2 の時間発展を調べよ。

5

行列計算

物理の問題には、連立方程式を解いたり固有値を求める問題に帰着する場合が多くある。 この章では、これら「行列を扱う数値計算」の基礎的なアルゴリズムを紹介する。

5.1

ガウスの消去法

ポアソン方程式とは ∇2_{φ =−ρ (5.1)} のような偏微分方程式をいう。例えば、右辺のρは電荷分布、左辺の φは静電ポテンシャ ルで、与えられた電荷分布から生じる静電ポテンシャルを求めるにはこの偏微分方程式を 適当な境界条件の下で解く必要がある。 簡単のため、１次元の場合を考えよう。２階微分を２階差分で近似すると φj+1+ φj−1− 2φj h2 =−ρj (5.2) となる。境界を j = 0とj = N + 1とすればφ0, φN +1 は境界条件から既に決まっている ので、ρj（j = 1,· · · , N） が与えられたとき、連立方程式を満たす内点での値φ1,· · · , φN を求める問題となる。 式(5.2)を行列形式に書けば      2 −1 0 · · · 0 −1 2 −1 · · · 0 · · · · 0 · · · 0 −1 2           φ1 φ2 · · · φN     =      h2ρ1+ φ0 h2ρ2 · · · h2ρN + φN +1      (5.3) となる。右辺は既知のベクトルであるから、形式的には左辺の行列の逆行列がわかれば、 問題は解けたことになる。ただし現実問題としては、逆行列の計算は計算量が多いので行 なわず、諸君もよく知っている消去法（行列の基本変形）で解く。 [例題5-1] 一様な電荷分布 ρ0 の場合の１次元ポアソン方程式（0 ≤ x ≤ L）を数値的に解け。こ の場合は φ(x) = ρ0 2 x(L− x) (5.4) が厳密解である。 [解答] ガウスの消去法を用いた数値解法のプログラムを以下に示す。ここで「gauss」と 名づけられたサブルーチンが消去法を実行している部分である。行列は N× N で、対角

成分の配列 diagと非対角左下成分の配列subおよび非対角右上成分の配列supとからな

る。このような行列を３重対角行列（tridiagonal matrix）という。右辺のベクトルは配列

b に収められ、計算結果が上書きされる。

Cでは、配列は j = 0,· · · , N − 1と番号付けられるので、b[j-1]=bj となっていること に注意する必要がある。ちょっと不便なのだが、致しかたない。

/* poisson_1d.c */ #include <stdio.h> #define N 100

/* ３重対角行列のガウス消去法 */

gauss(double diag[], double sub[], double sup[], double b[]){ int i;

double t;

for(i=0; i<(N-1); i++){ t=sub[i]/diag[i]; diag[i+1]-=t*sup[i]; b[i+1]-=t*b[i]; }

b[N-1]/=diag[N-1];

for(i=(N-2); i>=0; i--){

b[i]=(b[i]-sup[i]*b[i+1])/diag[i]; } } main(){ int i, j; double h=0.01;

double diag[N], sup[N], sub[N], b[N]; double phi[N+2], rho[N+2];

/* 境界条件と右辺の定義 */ phi[0]=0; phi[N+1]=0; for(j=1; j<=N; j++){ rho[j]=1.0; b[j-1]=h*h*rho[j]; } b[0]+=phi[0]; b[N-1]+=phi[N+1]; /* 行列の定義 */

for(i=0; i<N; i++){

diag[i]=2; sub[i]=-1; sup[i]=-1; }

/* 行列方程式を解く：配列 b に解が上書きされる */

gauss(diag, sub, sup, b); /* 出力 */

for(j=1; j<=N; j++){ phi[j]=b[j-1]; } for(j=0; j<=(N+1); j++){ printf("%d %lf\n", j, phi[j]); } } [問題5-1] N × N の行列 A =         2 −1 0 · · · 0 −1 2 −1 · · · 0 · · · · · · · · · 0 · · · −1 2 −1 0 · · · 0 −1 2         (5.5) の逆行列をN = 2, 3, 4の場合に計算してみよ。どんな規則性があるか、わかるだろうか？ （ヒント） 逆行列の各行列要素を行列式|A|（= N + 1 となる）倍すれば、すべて整数 になる。

5.2

LU

分解

前節のプログラムは３重対角行列にしか使えなかったが、一般には扱う行列が３重対角 であるとは限らない。そのような場合によく用いられるのが、ここで紹介する「LU分解」

の方法である。LUとは Lower triangular matrix（下三角行列）およびUpper triangular

matrix（上三角行列）を意味し、任意の行列A をふたつの三角行列の積 A = LU (5.6) の形に書き表そうというものである。もしこれが出来れば、 Ax = b → x = U−1L−1b (5.7) と解くことができる。三角行列の逆は（後述するように）代入操作だけで実行できるので、 LU 分解できれば、解けたも同然 なのである。一般に A の逆行列が存在するとき、すなわち行列式 |A| 6= 0 ならば、行列 A はLU 分解可能である。 例として、前述のような３重対角行列の場合を考えよう。このとき、      d1 f1 0 · · · 0 e2 d2 f2 · · · 0 · · · · · · · · · 0 · · · 0 eN dN     =      1 0 0 · · · 0 `2 1 0 · · · 0 · · · · · · · · · 0 · · · 0 `N 1     ·      a1 u1 0 · · · 0 0 a2 u2 · · · 0 · · · · · · · · · 0 · · · 0 0 aN      (5.8)

と置いて、両辺の各成分を等値すれば d1 = a1, f1= u1, ej = `jaj−1, dj = `juj−1+ aj, fj = uj, (j = 2,· · · , N − 1) eN = `NaN−1, dN = `NuN−1+ aN を得る。よって、順番に a1 = d1, u1 = f1, (5.9) `j = ej/aj−1, aj = dj − `juj−1, uj = fj, (j = 2,· · · , N − 1) (5.10) `N = eN/aN−1, aN = dN − `NuN−1 (5.11) と、L, U の各行列要素が求められる。したがって Ax≡ LUx = b → x = U−1L−1b (5.12) と解が求まる。ここで、上（下）三角行列の逆行列はまた上（下）三角行列となるのだが、 わざわざ逆行列を求めずとも 代入操作だけで解くことができる という点（これが LU 分解する理由）が重要である。 このことを、一般の下三角行列の場合で説明しよう。L = (`jk) に対してLy = b の解 は`jk = 0 (j < k) であるから `11y1 = b1, `21y1+ `22y2 = b2, · · · · `N 1y1+· · · + `N NyN = bN を上から順に y1= b1/`11, y2= (b2− `21y1)/`22, · · · · yN = (bN − `N 1y1− · · · − `N N−1)/`N N と、解いていけば良いのである。上三角行列の場合も同様である。このような解法を、と くに「前進消去・後退代入法」と呼ぶことがある。 一般の（３重対角でない）行列の場合には、(5.8)のように簡単ではないが、それでも LU 分解のアルゴリズムが存在する。その説明は長くなるので省略して、直接にLU 分解 のプログラムを示そう。

/* lu.c */

#include <stdio.h> #include <math.h> #define N 3

double ludcmp(double a[N][N], int ip[N]){ int i, j, k, ii, ik;

double t, u, det, weight[N];

det=0.0;/* 行列式 */ for(k=0; k<N; k++){

ip[k]=k; u=0.0;

for(j=0; j<N; j++){

t=fabs(a[k][j]); if(t>u) u=t; }

if(u==0) return 0; /* u=0 ならLU分解不能 */ weight[k]=1.0/u;

}

det=1.0;/* 行列式の初期値 */ for(k=0; k<N; k++){

u=-1.0;

for(i=k; i<N; i++){ ii=ip[i]; t=fabs(a[ii][k])*weight[ii]; if(t>u){ u=t; j=i; } } ik=ip[j]; if(j!=k){ ip[j]=ip[k]; ip[k]=ik; det=-det; } u=a[ik][k]; det*=u;

if(u==0) return 0;/* u=0 ならLU分解不能 */ for(i=k+1; i<N; i++){

ii=ip[i];

a[ii][k]/=u; t=a[ii][k]; for(j=k+1; j<N; j++){

} } }

return det; }

void solve(double a[N][N], double x[N], double b[N], int ip[N]){ int i, j, ii;

double t;

for (i=0; i<N; i++){ ii=ip[i]; t=b[ii]; for(j=0; j<i; j++){ t-=a[ii][j]*x[j]; } x[i]=t; }

for(i=N-1; i>=0; i--){ t=x[i]; ii=ip[i]; for(j=i+1; j<N; j++){ t-=a[ii][j]*x[j]; } x[i]=t/a[ii][i]; } } main(){/* ３×３ 行列の例 */ int i, j, ip[N];

double det, a[N][N], x[N], b[N]; /* 行列の定義 */

a[0][0]=2; a[0][1]=-1; a[0][2]=1; a[1][0]=-1; a[1][1]=2; a[1][2]=-1; a[2][0]=1; a[2][1]=-1; a[2][2]=2; /* 右辺の定義 */ b[0]=1; b[1]=0; b[2]=0; /* 解く */ det=ludcmp(a, ip); if(det==0){ printf("Unable to LU decompose!\n"); }else{ solve(a, x, b, ip);

for(j=0; j<N; j++){ printf("x[%d]=%lf\n", j, x[j]); } } } [問題5-2] 次の６×６行列の行列式および逆行列を計算せよ。 A =           2 −1 0 0 0 0 −1 2 −1 0 0 0 0 −1 2 −1 −1 0 0 0 −1 2 0 0 0 0 −1 0 2 −1 0 0 0 0 −1 2           （ヒント）行列式=3となる。逆行列を３倍すると、その行列要素はすべて整数になる。 [線形最小２乗法] LU 分解を用いる応用問題をひとつ示そう。物理実験のいろいろな場面で 実験データ (xj, yj)Nj=1 に理論式 y = f (x) をフィットさせる という問題が発生する。 とくに関数f (x) がパラメータ ak と関数 φk(x)（k = 1,· · · , n）を用いて f (x) = n ∑ k=1 akφk(x) (5.13) と表される、という「線形問題」を考えよう。 一般にデータ数はパラメータ数より多い（N > n）ので、実験誤差等によって実験値は 理論曲線に完全に一致することはない。この場合、理論値と実験値との差の２乗和 Q(a1,· · · , an) = N ∑ j=1 (yj− f(xj))2= N ∑ j=1 (yj− a1φ1(xj)− · · · − anφn(xj))2 (5.14) を最小にするようにパラメータa1,· · · , anを決めよ、という問題を線形最小２乗法という。 極値問題であるから、ak で偏微分して（k = 1,· · · , n） ∂Q ∂ak =−2 N ∑ j=1 (yj− a1φ1(xj)− · · · − anφn(xj))φk(xj) = 0 (5.15) より、a1,· · · , an は連立方程式 n ∑ `=1  ∑N j=1 φk(xj)φ`(xj)   a_`= N ∑ j=1 yjφk(xj), (k = 1,· · · , n) (5.16)

を解いて得られる。すなわち、行列形式で書けば M a = b, mk`= N ∑ j=1 φk(xj)φ`(xj), bk= N ∑ j=1 yjφk(xj) (5.17) を解く問題に帰着する。そして、この方程式は LU分解法により解ける。 [問題5-3] 自然数x が大きいとき、漸近的に log(x!) = log(√2πx ) + x log x− x (5.18) が成り立つ（対数は自然対数）。これをスターリングの公式という。 5個のデータ log(1!) = 0, log(10!) = 15.1044, log(100!) = 363.739, log(1000!) = 5912.13, log(10000!) = 82108.9, をもとに、スターリングの公式

log(x!) = a0+ a1log x + a2x + a3x log x (5.19)

の係数a0, a1, a2, a3 を最小2乗法により推定し、公式の値をどの程度再現するかを調べて みよ。 ［注意］ 本章で述べた「消去法・LU分解」などのアルゴリズムを使うには、上掲ソースを自分 のプログラムに取り込んで書くのも良いが、 定評のあるライブラリをインクルードして使う というのも良い考えである。お勧めは （１）奥村晴彦『Ｃによるアルゴリズム事典』：    http://oku.edu.mie-u.ac.jp/∼okumura/algo/ （２）LAPACK：    http://www.netlib.org/lapack/ （３）Numerical Recipes：    http://www.nr.com/ などである。ただし、ライブラリ化が面倒であったり、使用法などを学習する必要がある。 また、（３）は小額ではあるが有料である。

6

補追：偏微分方程式

6.1

２次元ポアソン方程式

１次元の場合は前章で取り上げたので、ここでは２次元のポアソン方程式 ∇2_{φ =−ρ, ∇}2₌ ∂2 ∂x2 + ∂2 ∂y2 (6.1) を解くためのアルゴリズムのひとつ「緩和法」を紹介しよう。 緩和法（逐次過緩和法：SOR=Successive Over Relaxation）

２次元ポアソン方程式を差分化すると φj+1,k+ φj−1,k− 2φj,k h2 + φj,k+1+ φj,k−1− 2φj,k h2 =−ρj,k (6.2) となる。緩和法では、これを次の漸化式 φn+1_j,k = φn_j,k+α 4 ( φn_j+1,k+ φn+1_j−1,k+ φn_j,k+1+ φn+1_j,k−1− 4φn_j,k+ h2ρj,k ) (6.3) による繰り返し計算によって（近似的に）解く。３次元の場合への拡張も容易に想像でき るであろう。 注意してほしいのは、 φn+1_j−1,k, φn+1_j,k−1 は既に計算した n + 1 での値を使う という点である。この逐次計算を充分な回数繰り返せば、φn+1+ φnとなり、このとき右 辺の括弧内 = 0 すなわちポアソン方程式が満たされるだろう、というのである。ここで α は「加速パラメータ」とよばれる係数で、最適な値に選べば収束が速くなるが、普通は 1 < α < 2の適当な値が使われる。この解法は、仮想的な時間t のもとでの拡散型方程式 ∂φ ∂t =∇ 2_{φ + ρ (6.4)} を離散化したものと考えられ、充分な時間が経てば平衡状態に達するというのは物理的に も妥当なことである。 [例題6-1] 正方形の領域内に一様な電荷分布 ρ0 がある場合のポアソン方程式を境界条 件 φ(境界) = 0 のもとで解き、結果を Gnuplotで描画（等高線、立体図）せよ。 [解答] 以下の通り。繰り返し計算の回数 max を変えて、緩和の様子を確認してみると おもしろい（例えば、最大値 phi[N/2][N/2]の推移を調べよ）。 /* poisson_2d.c */ #include <stdio.h> #define N 100 main(){

int j, k, n, max=1000; double phi[N+2][N+2]; double h=0.01, rho0=1;

double sum, residue, alpha=1.5; /* 初期化 */ for(j=0; j<=(N+1); j++){ for(k=0; k<=(N+1); k++){ phi[j][k]=0; } } /* SOR */ for(n=0; n<max; n++){ for(j=1; j<=N; j++){ for(k=1; k<=N; k++){ sum=phi[j+1][k]+phi[j-1][k]+phi[j][k+1]+phi[j][k-1]-4*phi[j][k]; residue=sum+h*h*rho0; phi[j][k]+=(alpha/4)*residue; } } } /* 出力 */ for(j=0; j<=(N+1); j++){ for(k=0; k<=(N+1); k++){ printf("%d %d %lf\n", j, k, phi[j][k]); } } } [問題6-1] （２次元ポアソン方程式２） 例題と同じ設定（一様電荷、境界値ゼロ）ではあるが、領域が正方形ではなく右上４分 の１が欠けた L 字形をしている場合のポアソン方程式を緩和法によって解け。

6.2

拡散方程式

こんどは実時間 tを含む偏微分方程式の例として、拡散方程式 ∂P ∂t = D∇ 2_{P (6.5)} を取り上げよう。この方程式は線形なので、形式的に P (t + ∆t) = eD∆t∇2P (t) (6.6) と解けることに注意する。 以下は簡単のため空間１次元の場合で説明しよう。∆tが小さいことを使えば eD∆t∇2x _{= 1 + D∆t} ∂ 2 ∂x2 (6.7) と近似できる。したがって、時間・空間を差分化して P_jn+1= P_jn+D∆t h2 ( P_j+1n + P_jn₋₁− 2P_jn) (6.8) を得る。これは１次元拡散方程式の陽解法を与える。しかしながら残念なことに、この陽 解法による計算法はあまり満足できる結果を与えないので、通常は以下のような陰解法が 用いられる。 式(6.6)を書き直し e−D∆t∇2P (t + ∆t) = P (t) (6.9) として、式(6.7)と同様の∆tが小さい近似を行えば P_jn+1−D∆t h2 ( P_j+1n+1+ P_jn+1₋₁ − 2P_jn+1 ) = P_jn (6.10) を得る。これは前章「ガウスの消去法」で扱ったのと類似の連立方程式系であるから、行 列計算によって解ける。ただし、いつも同じ行列が出てくるので、この場合には 最初に一度だけ逆行列を計算しておく のが賢い。 [例題6-2] １次元拡散方程式を周期的境界条件 P (x + L, t) = P (x, t) (6.11) のもとで解け。差分化は、LをN 等分して、周期的境界条件を P_Nn = P₀n, P_{N +1}n = P₁n (6.12) とするのが見やすい。計算すべき独立な変数は P_jn（j = 0, 1,· · · , N − 1）である。

初期条件は、例えば P (x, 0) = P0 ( 1−|2x − L| L ) (6.13) とせよ（３角形）。 [解答] 解くべき差分方程式は、 −D∆t h2 P n+1 j−1 + ( 1 +2D∆t h2 ) P_jn+1−D∆t h2 P n+1 j+1+ = P n j (6.14) であるから、d = 1 + 2D∆t/h2, s =−D∆t/h2 と置いて、行列形に表せば      d s 0 · · · s s d s · · · 0 · · · · · · s 0 · · · s d           P₀n+1 P₁n+1 · · · P_Nn+1₋₁     =      P₀n P₁n · · · P_Nn₋₁      (6.15) となる。よって、この行列方程式を LU 分解によって解くプログラムを書けばよい。 "dat" 0 ₅₀₀ 1000₁₅₀₀ 2000₂₅₀₀ 3000 3500₄₀₀₀ 45000 1020 3040 5060 7080 90100 0 50 100 150 200 時々刻々の(n, j, P_jn)の計算結果を出力しておき、ファイルをGnuplotに渡して３次 元的に表示する（上図および付録２を参照）ためのプログラムは、以下のようになる。 /* diffusion.c */ #include <stdio.h> #define N 100

double ludcmp(double a[N][N], int ip[N]){ /* P.28 lu.c の ludcmp と同じ */

void solve(double a[N][N], double x[N], double b[N], int ip[N]){ /* P.29 lu.c の solve と同じ */

}

main(){

int i, j, k, max=5000, ip[N]; double D, dt, h, d, s, P0, dummy; double det, a[N][N], P[N+1], PN[N+1]; /* パラメータ */ D=1.0; dt=0.001; h=0.1; d=1+2*D*dt/(h*h); s=-D*dt/(h*h); /* 初期条件 */ P0=200; for(j=0; j<=N; j++){ dummy=(double)(2*j-N); dummy=dummy/(double)N; P[j]=P0*(1-fabs(dummy));/* ３角形 */ } /* 行列の定義 */ for(j=0; j<N; j++){ for(k=0; k<N; k++){ a[j][k]=0; } } for(j=0; j<N; j++){ a[j][j]=d; } for(j=0; j<(N-1); j++){ a[j][j+1]=s; a[j+1][j]=s; } a[0][N-1]=s; a[N-1][0]=s; /* LU 分解 */

det=ludcmp(a, ip);/* det != 0 を確認しておく */ /* 時間発展 */

for(i=0; i<max; i++){ if((i%500)==0){ for(j=0; j<N; j++){ printf("%d %d %lf\n",i, j, P[j]); } printf("\n");/* 空白行を入れる */ }

solve(a, PN, P, ip); PN[N]=PN[0];/* 周期境界条件 */ for(j=0; j<=N; j++){ P[j]=PN[j]; } } } [問題6-2] （Gnuplotでアニメーション） 上記の例題の場合、出来れば 時々刻々の P_jn をアニメ描画したい と思うのは当然であろう。例題で得たのと同様のデータ・ファイルを用いれば、Gnuplot によるアニメ表示ができる。どうすれば可能かを調べて実行せよ。 Gnuplotを使う以外にも、画面操作を扱うライブラリ（XLib, GLなど）を組み込んで、

実行する方法がある。また直接に画面操作できる言語（Java, Basic, Perl, Python, Ruby など）もあるので、試してみるとよい。 高次元の拡散方程式 ２次元・３次元の場合には e−D∆t∇2P (t + ∆t) = P (t) (6.16) において e−D∆t(∇2x+∇2y)_{= e}−D∆t∇2x· e−D∆t∇2y_{, （２次元の場合） (6.17)} e−D∆t(∇2x+∇2y+∇2z)_{= e}−D∆t∇2x· e−D∆t∇2y· e−D∆t∇2z_{, （３次元の場合） (6.18)} であるから、前述の１次元解法を向きを変えて交互に実行すればよいことがわかる。これ を「ADI法」（ADI=Alternating DIrection）という。

付録１．

LINUX

（

UNIX

）お役立ち

Tips

1. ファイル名入力補完 長いディレクトリ名や長いファイル名を、キーボードから打ち込む場合がよくある。 綴りを間違えてやり直すという経験をした人も多いだろう。こんなとき「↑」キー を使えば履歴が出るので、労力軽減になり便利である。 ところで、上記の機能は使っていても「ファイル名入力補完」機能は知らない人が 多いようだ。名前がうろ覚えで不確かなときは「途中まで書いてタブキーを押す」 とよい。シェルが自動的に候補となる名前を補ってくれるので、とても便利である。 候補が複数ある場合は途中で止まるので、そこから続きを入力して改めてタブキー を押せばよい。 2. バックグラウンド・ジョブ ターミナル（xterm やkterm）をひとつ開いて、何かのコマンドを実行する場合を 想定してみよう。例えば、emacsを起動した場合、新しく作業用ウィンドウが開く。 こんなとき、編集作業中に別のツール（例えば gccや gnuplot）を使いたくなるこ とがよくある。こういう場合は、起動の際に $ emacs newton.c &

のように、最後に「&」を付けておけば、ターミナルが引き続き使えるので、emacs を終了したり、新しいターミナルを開かずにすむ。 3. リダイレクション 標準出力に出力された数値データを、ファイルに格納したいということがよくある。 こういう場合には、例えば $ ./a.out > magnet.dat としてやればよい（ファイル名は自由）。ファイルの読み書きを C 言語で実現する 方法は[例題4-1]の解答に示されているが、急ぐ場合はこの方法が簡単である。 4. emacs のショートカット 皆さんはエディタとして emacs を使っているが、「メニューをマウスで操作する」 ひとが多い。これでも良いが、いろいろなショートカット（近道）を知っていると、 もっと楽に使うことができる。メニュー欄の括弧内にも書かれているが、いくつか を紹介しよう。 ファイル保存 （save）Ctrl+x Ctrl+s 終了  Ctrl+x Ctrl+c

カット＆ペースト（cut＆paste） Ctrl+w Ctrl+y コピー＆ペースト（copy＆paste） Esc+w Ctrl+y シェルモード Esc+! または Esc+x shell

付録２：

Gnuplot

簡易マニュアル

（１）起動法、終了法 Gnuplotを起動するには、コマンドラインから $ gnuplot と入力すると、端末プロンプトがgnuplot>と入力待ちになるので、以後ここにコマンド 入力する。 終了は gnuplot> exit または gnuplot> Ctrl+D （コントロールキーを押したまま D を押す） とする。 （２）描画法 （例１）既知関数を描画する gnuplot> plot x**2 gnuplot> plot sin(x)

gnuplot> splot exp(-x**2-y**2) （例２）データから描画する

１行ごとに

  xの値 [スペースまたはタブ] yの値

の並んだファイル graph2d.dat に対して gnuplot> plot "graph2d.dat"

とすれば、２次元グラフが作成される。 同様に、

  xの値 [スペースまたはタブ] yの値 [スペースまたはタブ] z の値

の並んだ（ただしブロックごとに空行が必要）ファイル graph3d.datに対しては

gnuplot> splot "graph3d.dat" とすればよい。

目盛りの付け方・軸の名前・点の種類などの「描画スタイル」は、オプションで指定で きるので、詳しくはオンラインマニュアル

や参考書で調べてほしい。 （３）印刷・出力

[紙に出力]

プリンタで印刷するには

gnuplot> set terminal postscript

gnuplot> set output "| lpr -P lp1 [lp2, lp3, lp4]"

と設定しておいて、描画する。情報演習室には４台のプリンタがあるので、その内のひと つを選ぶ。

[ファイルに出力]

画像をファイルとして出力するには、例えば gnuplot> set terminal postscript eps gnuplot> set output "filename.eps"

とすれば、EPS（Encapsulated PostScript）ファイルとして出力される。例えば、_TEXファ イル中に \begin{center} \includegraphics[width=10cm,clip]{filename.eps} \end{center} と記述すれば _TEX文書の中に画像を挿入できる。 この他にもいろいろなオプションがあるので、調べてみると良い。