計算物理 WS

行列の反転とラプラス方程式の解法

行列の反転

1

前回までは常微分方程式を解いて粒子の運動を調べてきたが、次に取り組む問題は、Laplace 方程式の境界値問題を解くことである。ラプラス方程式は、定常状態の熱伝導方程式の解や電極にかこまれた真空中の静電場の問題、重力場を求める問題など、物理にしばしば現れる方程式である。まずその準備として、巨大な連立一次方程式を数値的に解く方法（行列反転と等価）について学ぶ。連立一次方程式を解く作業は多くの計算物理に現れる作業で、きわめて重要である。

1.1

行列の反転で連立一次方程式を解く

今回取り扱うのは、vector の表記では以下のような連立一次方程式である。

b = Ax(1)

ここで b は定数を成分とする列ベクトル、A はやはり定数を成分とする行列である。この方程式は成分をあからさまに書くと、


b1a11 a12 b2 
 a21 a22 . = .
..
 ..  ..
.
bn
an1 an2

x1· · · a1n
· · · a2n   x2 .
..
 .. . ..  . xn· · · ann(2)とあらわされる。数学的にこの方程式に解があるためにはまず det A がゼロでない値を持たなければならない。ゼロではない場合には逆行列 A−1 が存在し、解は A−1 を用いて

x = A−1 b(3)と書くことができる。
n が小さい問題であれば、直接この逆行列の成分を計算し、解を求める方が効率的であるが、これから取り扱うような n の大きな問題では、逆行列を求めるのではなく、反復によって式 (1) を満たす x を探す方法が計算を早く進めることができる。今回は反復法を学ぶ。

1.2

ヤコビの反復法

ヤコビの反復法は最も単純な反復法である。手順は• まず試しの値を解として x にいれ、それを x(0) とする。全部ゼロとかでよい。

(0)• 各行の方程式を xk に関して解く。このとき、他の xi はすべて最初の試行の値 xi である。すなわち(1)x1 =..
.

(1)

xk =..
.

{({(a11 −1 b1 −(0)(0)(0))}(0))}+a12 x2 · · · + a1k xk · · · + a1n−1 xn−1 + a1n x(0)n..
.

akk −1 bk − ak1 x1 + ak2 x2 · · ·(0)(0)· · · + akn−1 xn−1 + akn x(0)
n..
.
{
((0)
(0)
(0)
(0)
(1)
−1
xn = annbn − an1 x1 + an2 x2 · · · + ank xk · · · + ann−1 xn−1
(1)(4))}(2)• 前のステップで求まった xk を用いて同じことをすると xk が求まる。これをどんど(i)ん繰り返して xk を計算していく。
(i)

(i+1)

• xk と xkがあらかじめ設定した誤差の範囲内に収まったら、繰り返しのループを抜け、その値を答えとする。
次に実際のプログラムを見てみる。→トップページから。プログラムは以下の 3 変数の連立 1 次方程式を解くコードになっている。








63 −1 2x1
 8  =  1 5 −4   x  2 21 −4 7x3

1.3(5)

ヤコビ法の収束条件

反復法ではいつも必ず解が求まるわけではない。問題によっては答えが収束せず、発散していく場合がある。今回の例のプログラムではうまく解けるが、たとえば




63 −1 2x1 8 = 15−4x2

x3210 −4 7(6)とすると、すでに解けなくなる。ヤコビ法が（実はこのあとのガウス−ザイデル法も）収束するための十分条件が知られており、|aii | >n∑|aij |(7)j=1,j̸=iであれば必ず収束する（十分条件なのでこうでなければ収束しないというわけではなく、これに近いものは収束する場合が多い）。この条件は対角要素が非対角要素よりも大きいという意味で、対角優位な行列、という言い方をする。物理の問題では多くはこの条件を満たすものが多い。

1.4

ガウス-ザイデル法

ガウス−ザイデル法はヤコビの反復法を少し賢くし、解の収束を早めた方法である。以下に手順を示す。

• まず試しの値 x(0) を設定する。
(0)• 各行の方程式をまず x1 に関して解く。このとき、他の xi はすべて最初の試行の値 xi(0)である。ヤコビ法との違いは x1 の次の x2 に関して解く時に、x1 ではなく、すぐ前(1)(1)に計算した x1 を用いることである。つまり xk を計算するときには j < k である j(1)に対しては xj を使っていくことによって収束を早めるのである。すなわち{({(−1(1)b1 −x1 = a11..
.
(1)
xk =(0)akk −1 bk − ak1 x1 + ak2 x2 · · ·..
.

−1x(1)
n = ann{(0)(0))}(0))}+a12 x2 · · · + a1k xk · · · + a1n−1 xn−1 + a1n x(0)
n..
.
(1)

(1)

· · · + akn−1 xn−1 + akn x(0)
n..
.
(1)
(1)
(1)
(1)
bn − an1 x1 + an2 x2 · · · + ank xk · · · + ann−1 xn−1((8))}(2)(1)• 前のステップで求まった xk を用いて同じことをすると xk が求まる。これをどんど(i)ん繰り返して xk を計算していく。ここもヤコビ法と同じ。
(i)

(i+1)

があらかじめ設定した誤差の範囲内に収まったら、繰り返しのループを• xk と xk抜け、その値を答えとする。ここもヤコビ法と同じ。

• ガウス―ザイデル法のプログラムはヤコビ法のコードをほんのわずかに変えるだけでできる。どこを変更すればよいであろうか。

1.5

ラプラス方程式の差分化

以下の 2 次元ラプラス方程式を与えられた境界条件でとくことを考えよう。

∂ 2ψ ∂ 2ψ
+ 2 =0∂x2∂y(9)次のような境界条件の問題を考える。x = 0,y = 0, L を接地し、x = L のみに ψ0 の電圧をかけた時の電位を計算する。
数値的にこの方程式を解くとは、ψ を xy 平面上の格子点での値として求めることである（格子間隔を x 方向は ∆x,y 方向は ∆y ）。
さて今ある格子点 (xi+1 , yi ) でのポテンシャルの値 ψ(xi + ∆x, yj ) を x 方向にテーラー展開すると、

ψ(xi + ∆x, yj ) = ψ(xi , yj ) +∂ψ
1 ∂ 2ψ
1 ∂ 3ψ
2∆x +∆x+∆x3 · · ·∂x2! ∂x23! ∂x3(10)同様に格子点 (xi−1 , yi ) では

ψ(xi − ∆x, yj ) = ψ(xi , yj ) −1 ∂ 2ψ
1 ∂ 3ψ
∂ψ2∆x +∆x−∆x3 · · ·∂x
2! ∂x23! ∂x3(11)

このふたつの式を足すと ∆x の奇数次の項が消えてψ(xi − ∆x, yj ) + ψ(xi + ∆x, yj ) = 2ψ(xi , yj ) +∂ 2ψ
∆x2 + · · ·∂x2

(12)∆x の 4 次以上を無視すると、2 階微分が格子点の ψ の値から評価でき、
1∂ 2ψ
={ψ(xi − ∆x, yj ) + ψ(xi + ∆x, yj ) − 2ψ(xi , yj )}
2∂x
∆x2

(13)

y 方向についてもまったく同様の議論をすると、∂ 2ψ
1={ψ(xi , yj − ∆yj ) + ψ(xi , yj + ∆y) − 2ψ(xi , yj )}2∂y∆y 2(14)上の 2 つの式を加えるとラプラス方程式の左辺が格子点のポテンシャルの値で表現できる。
ラプラス方程式を ψij ≡ ψ(xi , yj ) という記法を使って書き直すと

1
1{ψi−1 j + ψi+1 j − 2ψi j } +{ψi j−1 + ψi j+1 − 2ψi j } = 02∆x∆y 2(15)∆x = ∆y とすることにし、ψi j に関して解くと、
ψi j =1(ψi−1 j + ψi+1 j + ψi j−1 + ψi j+1 )4(i, j = 1, 2, · · · , n)(16)となり、ある格子点 i, j でのポテンシャルの値は単にその周りの 4 格子点でのポテンシャルの値の平均である、という単純な一次式である差分方程式が得られる。

1.6

差分方程式＝連立一次方程式

さて格子点が端を入れて全部で (n + 2)2 個あるとし (図 1 では n = 4)、端の値は境界条件で ﬁx されているとすると、未知の ψi j の値は全部で n2 個あることになる。一方、方程式は式 16 でやはり未知の ψi j の分で n2 個あるので、結局式 16 は n2 元連立の一次方程式を与えている。式 16 の右辺の中には境界の値 ψ0 i , ψn+1 i , ψi 0 , ψi n+1 などが含まれており、これらの方程式を解くには境界条件が必ず必要であることに留意しよう。
さてこの方程式はすでに未知数 ψi j に関して解かれた形になっているので、前回の連立一次方程式を解くときに、反復法で用いた式(1)x1 =..
.

(1)

xk =
..
.
x(1)
n =

{

({(a11 −1 b1 −(0)(0)(0))}(0))}+a12 x2 · · · + a1k xk · · · + a1n−1 xn−1 + a1n x(0)n..
.

akk −1 bk − a11 x1 + a12 x2 · · ·(0)(0)· · · + a1n−1 xn−1 + a1n x(0)n..
.
{((0)(0)(0)(0)−1ann
bn − a11 x1 + a12 x2 · · · + a1k xk · · · + a1n−1 xn−1)}

(17)

6

ψ(x, L) ≡ ψ0y(L, L)

ψi j+1-

ψ(0, y) ≡ 0

ψi−1 j?

6ψi jψi+1 jψ(L, y) ≡ 0

ψi j−10ψ(x, 0) ≡ 0
図 1: 格子点x -をすでに与えていることになる。よって行列 A やベクトル b の形を特に気にする必要は無
く、式 16 を、反復法の中身の式として使ってしまえばよい。次に実際のプログラムを見てみることにする。→トップページ
プログラムの出力がでたら、• gnuplot の splot コマンドや python を使って立体的に可視化してみよう。
• y = 0 に電圧をかけた場合の解や、他の場所に電圧をかけた場合の解も見てみよう。

1.7

Laplace 方程式の変数分離による解法

2 次元ラプラス方程式の数値解を求めたが、変数分離による方法で解いた解と比べてみよう。
∂ 2ψ ∂ 2ψ
+ 2 =0(18)∂x2∂y同様の境界条件の解を求めてみる（図 2）。物理的には x = 0,y = 0, L を接地し、x = L のみに ψ0 の電圧をかけた時の電位を計算することに対応する。

このような問題を解析的に解く場合には普通、変数分離法を用いる。変数分離法では解の 2 つの変数（この場合には x と y ）に対する依存性が掛け算で分かれてしまうという形を仮定するところから出発する。すなわち、

ψ(x, y) = ψx (x)ψy (y)これを元のラプラス方程式 (28)に代入して整理すると、

yL

6

ψ(x, L) ≡ ψ0

(L, L)

ψ(0, y) ≡ 00

ψ(x, 0) ≡ 0

ψ(L, y) ≡ 0L-x

図 2: 境界条件

1 d 2 ψx
1 d 2 ψy=−ψx dx2
ψy dy 2(19)するとこの左辺は x のみの関数、右辺は y のみの関数となるので、結局ある定数 A に等しくならざるを得ない。以下の A の正負によって場合わけをし、境界条件に適合する解を探していく。

(1) A = 0 のとき
1 d 2 ψx
1 d 2 ψy=−=0
ψx dx2ψy dy 2

であるので、それぞれ簡単に積分が実行でき、a, b, c, d を定数として、

ψx = ax + b(20)ψy = cy + d(21)ψ(x, y) = (ax + b)(cy + d)

(22)

よってであるが、これに境界条件 ψ(0, y) ≡ 0, ψ(L, y) ≡ 0 を課すとb(cy + d) ≡ 0(aL + b)(cy + d) ≡ 0となる。1 番目の式から

b = 0 or c = d = 0c = d = 0 の時には ψ ≡ 0 となる。また b = 0 でかつ cy + d ≡/ 0 であると、2 番目の式から
aL = 0 であるが、L ̸= 0 であるので、a = 0 となる。したがって ax + b ≡ 0 であり、やはり ψ ≡ 0 となる。
したがっていずれにせよこの解が与えられた境界条件を満たすためには ψ(x, y) ≡ 0 しかない。これは今興味のない trivial な解なので捨てることにする。
(2) A = k 2 (> 0) のとき (ここで k は正の実数)
1 d 2 ψy
1 d 2 ψx=−= k2
ψx dx2ψy dy 2であるので、まず ψx の方を考えると、解はやはり簡単に積分できて、

ψx = aekx + be−kxしたがってψ(x, y) = (aekx + be−kx )ψy (y)であるが、境界条件 ψ(0, y) = ψ(L, y) ≡ 0 から、ψy ≡/ 0 であれば、a + b = 0 かつ aekL +be−kL = 0 でなければならない。この 2 式から、b を消去すると、a(ekL − e−kL ) = 0となるが、kL = 0 という意味のない場合を除くと必ず a = 0 となり、したがって b = 0 も結論される。したがって ψ(x, y) ≡ 0 となる。ところでこれは ψy ≡/ 0 という仮定の下に得られた結果であるが、ψy ≡ 0 であればそもそも ψ(x, y) ≡ 0 であるので、やはりこの場合も興味のない解となることがわかったので捨てる。

(3) A = −k 2 (< 0) のとき (ここで k は正の実数)
まず ψx に関して解くと、
ψx = a sin (kx) + b cos (kx)

先ほどと同様 ψy ≡ 0 は意味がないので、ψy ≡/ 0 とすると、境界条件 ψ(0, y) ≡ 0 より、まず b = 0 が得られる。続いて ψ(L, y) ≡ 0 より、a sin(kL) = 0 であるが、a = 0 の無意味な解を除くと、k に関して条件がつくことがわかる。すなわちk = kn ≡nπL(n = 1, 2, 3, · · ·)(23)

次に ψy を解くと

ψy = cekn y + de−kn yが得られるが、境界条件 ψ(x, 0) ≡ 0 より、c = −d となる。よって

ψy = c(ekn y − e−kn y )あわせると、

ψ(x, y) = ac sin(kn x)(ekn y − e−kn y )であるが、n は自然数ならば何でもよく、かつ元のラプラス方程式が線形方程式なので、解の和（線形結合）もまた解となることを考慮すると、もっとも一般的な解の形は(∞∑)

nπnπnπan sin
ψ(x, y) =x (e L y − e− L y )Ln=1(24)ac はひとつの係数として新たに an とし、kn は定義どおりに書き下した。この形の解に最後の境界条件 ψ(x, L) = ψ0 を課す。すると(∞∑)

nπan sinx (enπ − e−nπ )ψ0 =Ln=1(25)であるが、これから未知の an を決定することができる。実際には三角関数の直交性、∫ L()

nπLsinx dx =L20()()∫ Lnπmπsinx sinx dx = 0LL02(26)(27)を用いる（上の式は Mathematica か手計算で確認する、ただし m ̸= n）。境界条件 (25) の両辺に sin (mπx/L) をかけて積分すると、m ̸= n の項はすべて消え、

ψ0LL(1 − cos (nπ)) = an (enπ − e−nπ )nπ2となり、an に関してとくことができる。解いた an を用いると、結局

ψ(x, y) =∞∑()nπnπnπ2 1 − (−1)nx (e L y − e− L y )sinnπ−nπ−eLn=1 nπ eが得られる。理解できたらこの級数解を計算し、可視化して先ほどの数値解と比較してみよう。

1.8

ポアソン方程式

ポアソン方程式はラプラス方程式の右辺に電荷密度がある場合であるがこの場合も同様にとくことができる。2 次元の場合の差分方程式は

ϕij =11 ρij(ϕi−1j + ϕi+1j + ϕij−1 + ϕij+1 ) +(∆x)2
44 ϵ0

ここで ρij は (xi , yj ) での電荷 (面) 密度を表している。いま ij の点に電荷 Q を置くのであれば、電荷の保存の式から Q = ρij (∆x)2 である。理解できたら適当な境界条件と電荷分布から電場を求めて見よ。

1.9

提出課題

3 次元のラプラス方程式∂ 2ψ ∂ 2ψ ∂ 2ψ
+ 2 + 2 =0∂x2∂y∂z(28)を立方体内部で解くコードを作れ。自分で適当な境界条件を置いてみて電位がどうなるかを可視化してみよ。

提出期限：次回授業開始時