私たちの目にする空間は3次元の空間です．より正確にいうならば，局所的には3次元に見えると言った方がいいのかもしれません．私たちの目には地球は平面に見えますし，第一，宇宙の果てがどんなところかは想像さえつきません．しかし，こんな煩わしい議論をここでしようというのではありません．ともあれ，私たちの目にする空間は3次元の空間としましょう．物理現象を捕らえるには，この空間で考えるのは必ずしも適切ではありません．というのは，物質は動いているからです．この動きも込めて空間を設定する方がより自然に物理現象を捕らえることができるようになります．

1つの質点の運動なら，その位置が3次元，そして速度が3次元で計6次元の空間の1点であるとみなし，2つの質点の運動なら，それぞれの位置が 3次元，そしてそれぞれの速度が3次元と考えて，合計12次元の空間の1点であるとみなします．この空間を相空間といいます．

こんなばかでかい次元の空間を考えるとどんなメリットがあるのでしょうか．ニュートンの運動方程式F=maは，力が加速度と質量に比例することを表しています．速度は位置を微分したもの，加速度は速度を微分したものですが，加速度は力によって定まることをこの式は示しています．力は質点の位置と速度で定まりますから，相空間における現在の位置がわかれば，加速度のみたす微分方程式がたてられるので，この微分方程式を解けば(解けるかどうかはわかりませんが，少なくとも原理的には)質点の運動がすべて定まることがわかります．私たちの目にする空間での運動も相空間の点から求めることができます．

1.2 Newton の運動方程式

1.2.1 天動説と地動説

天動説はエジプト人のプトレマイオスにより始まると言われ，太陽を中心にすべての惑星が回っているというものです．これでは惑星(古くは遊星ともいった)の逆行が説明できないので，太陽を中心とする従円があり，従円の点を中心とする周転円の上に惑星が乗っているというという理論です．天

(4)

体の動きは綿密な観測の元に膨大なデータの積み重ねがあり，これを説明するために，周転円の上にさらに周転円を積み上げる必要がでてきて，理論はどんどん複雑さを増していきました．16世紀になったとき，それまでにカエサルの作ったユリウス歴(−45年制定)，は誤差が積み重なり，実際の春分と暦の春分とは10日余りの誤差がでていました．ユリウス暦はアウグストスにより訂正が行われ，それが用いられていて(cf. JulyとAugust)，1年を365 日，4年に1度閏年をもうけるという1年が365.25日とするものです．この誤差を訂正するため，多くの英知が集められ暦の訂正が行われました．これがグレゴリオ暦です．この中にコペルニクス(1473-1543)がいました．グレゴリオ暦は閏年を100年に一度は止め，さらに400年に1度は行うというもので，これで1年は365.2425日になり，実際の1年の365.2422日に限りなく近づきました．こうなれば，3000年に1日の誤差でこれぐらいになれば地球の自転の変化などの方が大きくなりそうです．そして，この暦は1582年に発表されました．カソリックの指示に従うのを快く思わないプロテスタントの国では適用がかなり遅れたようです(イギリスなどでは1752年，日本では 1873年，さらに東方教会では1923年)．

コペルニクスはおそらくこの研究の中で，地動説を発見したものと思われます．では，なぜこの地動説が受け入れられなかったのでしょうか．宗教的な理由だというのは，ガリレオ(1564-1642)の宗教裁判を思い浮かべる後世の勝手な解釈でしょう．第一，コペルニクスは司祭でした．ガリレオの裁判にしても最近の解釈では権力争いの]結果ではないかと言われています．実際に彼は1回目の裁判で無罪を勝ち取りましたが，捏造によって2回目の裁判で有罪を言い渡されたらしいのです．

理由はもっと単純です．地動説の方が誤差が大きかったのです．コペルニクスはその誤差を訂正するために，天動説の周転円を採り入れましたが，その結果，必要な周転円の数は天動説よりも多くなってしまったという笑えない話だそうです．これが訂正された功績はケプラー(1571-1630)にあります．

しかし，彼も師であるチコブラーエ(1546-1601)の膨大な観測データなしには理論を作ることができなかったことでしょう．その当時の観測装置を想像してください．その粗末な装置でのデータですから誤差が多く含まれているはずです．それにもかかわらず，ケプラーはわずかな差に着目して，惑星の起動が楕円であることを見出しました．楕円といってもほぼ完璧な円なのです．信じられない洞察力です．

(5)

1.2. Newtonの運動方程式 5

表1.1: 惑星の軌道惑星離心率公転半径10⁹km

水星 0.2056 0.057 金星 0.0068 0.108 地球 0.0167 0.150 火星 0.0934 0.228 木星 0.0484 0.778 土星 0.0542 1.427 天王星 0.0461 2.871 海王星 0.0085 4.498 冥王星 0.2488 5.914

そして，ケプラーは3つの原理を見出しました．

1. 惑星の軌道は楕円軌道である．

2. 一定の時間に通過する扇型の面積は一定である(面積速度一定)．

3. 公転半径の3乗は周期の2乗に比例する．

ニュートン(1642-1727)はこの原理を統一する原理として次の3つの原理を見出しました．

(6)

表 1.2: 水星，火星，金星，地球，木星 1. 力が働いていない物体は等速直線運動をする．

2. 力は加速度に比例する(F =ma)． 3. 作用反作用の原理

この3つの原理と2つの物体(それぞれ質量mとM)に働く力は距離rの2 乗に反比例するという，万有引力の法則

F =γmM r²

を見出し，これにより，ケプラーの原理を説明することができました．

1.2.2 微分

ニュートンは微分，すなわち微小な時間の変化という概念に行きつくことで，万有引力の法則からケプラーの原理を導きました．おそらく，微小な時間変化を記述する法則は長い時間をみるより単純であろうと考えたのでしょ

(7)

1.2. Newtonの運動方程式 7 う．x(t)で時刻tの位置を表すとすると，tからt+ ∆の間の速度(平均速度) は

x(t+ ∆)−x(t)

∆

であることは日常的に使われている概念でしょう．一定の速度で運動している場合には，この値が瞬間の速度と同一であるが，時間とともに変化する場合にも瞬間の速度というものを定めたければ，上の値の∆ →0としたときの極限とすればよいであろう，というのがニュートンの考えたことでしょう．

すなわち，時刻tにおける速度v(t)は v(t) = lim

∆→0

x(t+ ∆)−x(t)

∆

とおけばよいことになります．さらに，速度の変化量，加速度も同様に a(t) = lim

∆→0

v(t+ ∆)−v(t)

∆

と考えればよいことになります．さて，第2法則より，力は加速度に比例するのだからF =m×a(t)となります．

1.2.3 落下の法則

地表面での物体の落下を考えましょう．地球の半径に比べれば，数メートルなどとるにたらないものであるので，落下の場合にはrは一定とみなしてよいでしょう．地表面を平なものとみなし，

g=γM r²

とおくと，力は地球の中心に向かって働くことを考慮に入れると，第2法則はy軸方向の加速度はmay(t) =−mg，すなわちay(t) =−gを得ます．一方，x軸方向には力が働いていないのだから，ax(t) = 0を得ます．ずっと加速度が一定なので，v_x(t) =vx(0)，v_y(t) = vy(0)−gtであることがわかります．v_x(0), v_y(0)はx軸，y軸方向への初速度です．そこで最初の位置を (x(0), y(0))とすると，

x(t) =x(0) +v_x(0)t y(t) =y(0) +vy(0)t−1

2gt² を得ます．これから，tを消去すると

y(t) =y(0) +v_y(0)

(x(t)−x(0) v_x(0)

)

−1 2g

(x(t)−x(0) v_x(0)

)2

と放物線を描いて物体が飛ぶことがわかります．

(8)

1.2.4 惑星の運動

太陽を中心の極座標で考えましょう．

( x(t) y(t) )

=r(t) (

cosθ(t) sinθ(t) )

. 微分して

d dt

( x(t) y(t) )

= dr dt(t)

( cosθ(t) sinθ(t) )

+r(t)dθ dt(t)

(−sinθ(t) cosθ(t)

)

d² dt²

( x(t) y(t) )

= d²r dt²(t)

( cosθ(t) sinθ(t) )

+ 2dr dt(t)dθ

dt(t)

)

+r(t)d²θ dt²(t)

)

−r(t) (dθ

dt(t) )2(

cosθ(t) sinθ(t) )

= (

d²r

dt²(t)−r(t) (dθ

dt(t) )2) (

cosθ(t) sinθ(t) )

+ (

2dr dt(t)dθ

dt(t) +r(t)d²θ dt²(t)

) (−sinθ(t) cosθ(t)

)

万有引力γ^mM_r2 は− (

cosθ(t) sinθ(t) )

方向に働くことと (

cosθ(t) sinθ(t) )

と

)

が直交することから，k=γM とおくと第2法則は d²r

dt²(t)−r(t) (dθ

dt(t) )2

=− k

r(t)² (1.1)

2dr dt(t)dθ

dt(t) +r(t)d²θ

dt²(t) = 0 (1.2)

(1.2)は

d dt

( r(t)²dθ

dt(t) )

= 0 すなわち，ある定数hを用いて

r(t)²dθ dt(t) =h

と表せることがわかります．この右辺は微小時間変化したときの扇型の面積の2倍を表しているので，これがケプラーの第2原理です．(1.1)の微分方程式はp= ¹_rと置き換えると，r²^dθ_dt =hに注意すると

dp dθ =−1

r² dr dθ =−1

r² dr dt

dt dθ =−1

h dr dt d²p

dθ² = 1 h

d²r dt²

dt

dθ =−d²r dt²

1 h²p²

(9)

1.2. Newtonの運動方程式 9 により，比較的簡単な線形の微分方程式

d²p

dθ²+p= k h² を得ます．この解の1つは

p= k h² であり，

d²p

dθ² +p= 0

の解は，θ= 0のときp= 0となるように角度を定めれば p=Acosθ

と表されるので，2つをあわせて p= k

h² +Acosθ したがって

r= h²/k 1 +Ah²/kcosθ

を得ます．Ah²/k=eとおけば，これは離心率eの2次曲線を表しています．

実際の惑星では離心率は0とみなしてよいのでr=h²/k，これを面積速度一 定の法則に代入すると

dθ dt = h

r² =k^1/2 r^3/2

を得ます．両辺を1周期積分すると回転の角度は2πになるので，周期をT とすると

2π= k^1/2 r^3/2T を得ます．すなわち

(2π)²r³=kT² とケプラーの第3法則を得ます．

1.2.5 質点の落下

もっとも単純な例ですが，質点の落下を考えてみましょう．落ちるだけなのですから，位置も3次元で考える必要はないでしょう．高さの方向の1次元だけにしましょう．同様に速度方向も縦方向のみを考えれば，質点の運動は2次元で記述できることになります．質点の質量をmとすれば，下向きに重力mgが働いているから運動方程式は

x^′(t) = v(t) mv^′(t) = −mg

(10)

表 1.3: 地球の軌道，楕円??

(11)

1.2. Newtonの運動方程式 11

v x

図1.1: 相空間における質点の運動

と表されることになります．これは容易に解けて，時刻0での位置と速度を x(0) =x0，v(0) =v0とすると

x(t) = x₀+v₀t−1 2gt² v(t) = v0−gt

となります．したがって，tを消去すると質点の落下の運動は x=x0+v₀²

2g −v² 2g

と質点の軌道は図 1.1のように(x, v)平面の放物線として表されることがわかりました．

1.2.6 ふりこの運動

もう1つ，単純なモデルであるふりこの運動について述べておきましょう．

1つの平面を移動するふりこは鉛直線からの角度θとその角速度を考えれば記述できmす．ふりこのおもりをm，長さをlとすると(図1.2)

lθ^′ = v

mv^′ = −mgsinθ

となります．このままでは解けませんが，θが十分に小さいときにはsinθ∼θ とみなせるので，上の方程式は

lθ^′ = v mv^′ = −mgθ

(12)

θ

mg θ

mgsinθ 図1.2: ふりこの運動

で近似できるはずです．この方程式の解は ξ=lθ+i

√g lv とおけば

ξ^′=−i

√g lξ となるので

θ(t) =Asin (√g

lt+ω0

)

で与えられます．ここでAやω0は最初の位置や速度で決定されます．予想通り，ふりこは図 1.3のように周期運動がすることが導かれました．2次元

図 1.3: sinθ=θとみたときの，ふりこの相空間での運動

の相空間の各点(x, v)からその速度にあたるベクトル(x^′, v^′)を生やしてみま

(13)

1.3. 微分方程式とベクトル場 13 しょう．相空間の運動はこの速度ベクトルが接線になるように移動するわけですから，この図をみるだけで大まかな運動がわかることになります．振り子では図1.3のようになり，回転運動が見えています．これをベクトル場といいます．

図 1.4: 本当のふりこの相空間での運動

もっとも，ふりこの振り幅θが大きくなるとsinθ∼θとみなすことができなくなりますので，その解を簡単には表すことはできません．本当の解は微分方程式を解くことができないので記述はできませんが，図 1.4のようにベクトル場に表すとなんとなく解が見えてくることがわかるでしょう．θが大きいところでは図 1.3と大きく異なる挙動をします．

1.3 微分方程式とベクトル場

以上の例からもわかるように，位置だけではなく速度も含めて考えると見通しがよくなることがわかりました．

2次元の微分方程式を考えましょう．

dx

dt = f(x, y) dy

dt = g(x, y)

微分方程式を解くとは，この式をみたす(x(t), y(t))を求めること，すなわち x^′(t) = f(x(t), y(t))

y^′(t) = g(x(t), y(t))

を求めることです．しかし，通常の方程式は解けることはほとんどありません．そこで，せめて与えられた点(x0, y0)の近所の行動だけでもわかる努力をしてみましょう．それにはテイラー展開を用います．まず，(x₀, y₀)の近く

(14)

では

dx

dt = f(x0, y0) dy

dt = g(x0, y0)

と考えられるはずです．これは容易に解けて，t = 0で(x0, y0)にいるとすると

x(t)−x0) = f(x0, y0)t y(t)−y0 = g(x0, y0)t

と直線的に動くことがわかります．問題はf(x0, y0) =g(x0, y0) = 0となってしまった場合です．この場合には今の近似は成り立たなくなります．この場合を特異点といいます．こうなったら，テイラー展開を一歩すすめて

dx

dt = f(x0, y0) +ax(t) +by(t) =ax(t) +by(t) dy

dt = g(x0, y0) +cx(t) +dy(t) =cx(t) +dy(t) まで考えればいいでしょう．ここで

a=∂f

∂x(x0, y0), b= ∂f

∂y(x0, y0), c= ∂g

∂x(x0, y0), d= ∂g

∂y(x0, y0) です．これはベクトルx

¯= (

x y )

と行列

A= (

a b c d

)

を用いると

d dtx

¯=Ax

¯ と線形の微分方程式になります．

まず，もっともやさしい例から始めましょう．それは行列Aが対角行列 (

λ 0 0 µ

)

の場合です．この場合は，もとの形に戻せば dx

dt =λx dy dt =µy

(15)

1.3. 微分方程式とベクトル場 15 となって，2つの独立な方程式が得られることになる．この時刻t= 0における初期値が(x₀, y₀)である解は

x(t) =x₀e^λt y(t) =y0e^µt と解くことができます．次に簡単なのはA=

( λ −µ

µ λ

)

の場合です．この場合には極座標r=√

x²+y²とtanθ=^y_xに変換します．微分方程式は dx

dt =λx−µy dy

dt =µx+λy ですので，

dr dt = d

dt(x²+y²)^1/2)

=1

2(x²+y²)⁻^1/2 (

2xdx

dt + 2ydy dt

)

=1

r(x(λx−µy) +y(µx+λy))

=1

rλ(x²+y²) =λr となります．一方，

dtanθ dt = d

dt y

x= xdy/dt−ydx/dt x²

=x(µx+λy)−y(λx−µy)

x² =µx²+y² x² および

dtanθ

dt = 1

cos²θ dθ

dt = (1 + tan²θ)dθ dt

= x²+y² x²

dθ dt ですので，

dθ

dt =x²+y² x²

dtanθ

dt = x²

x²+y² ×µx²+y² x² =µ となります．以上から

dr

dt =ar, dθ dt =µ

を得ます．これは容易に解けて，0での初期値r=r₀とθ=θ₀のときの解は r(t) =r₀e^λt

θ(t) =θ0+µt

(16)

となります．すなわち，角度は一定の角速度µで回り，rは指数的に変化することがわかります．λ <0のときには原点に回転しながら収束し，λ >0なら回転しながら無限に発散することがわかります．このとき原点を渦状点とよび．λ <0のときには安定，λ >0のときには不安定といいます．λ= 0のときには常に一定の円の上を回転するので，原点を渦心点といいます．

いま1つやっておきたいのが (

λ 1 0 λ

)

の場合である．この行列は，ジョルダン標準形とよばれます．この場合はもとの形に戻せば

dx

dt =λx+y dy

dt =λy

となります．下の式はすぐに解けてy(t) =y0e^λt，これを上の式に代入すると dx

dt =λx+y0e^λt

となり，これは線形微分方程式なので，定数変化法で解を x(t) =C(t)e^λt

とおくと，

dx

dt =C(t)λe^λt+C^′(t)e^λt

=λx(t) +C^′(t)e^λt

により，C^′(t) =y₀を得ます．初期値を考慮すればC(t) =x₀+y₀tなので，

x(t) = (x₀+y₀(t−t₀))e^λt y(t) =y₀e^λt

であることがわかります．ベクトル場を用いて解の概略を考えましょう．

1. λ, µ >0のとき 2. λ, µ <0のとき

3. λ >0, µ <0，またはλ <0, µ >0のとき 4. ジョルダン標準形でλ >0のとき

5. ジョルダン標準形でλ <0のとき

この他にλ, µが0に等しい場合もありますが，面白くないので省略しましょう．これらのベクトル場はそれぞれ，図1.5のようになります．(1)，(2)，(4)， (5)のときは原点を結節点，そのうち時間の経過とともに原点に近づいてくる

(17)

1.3. 微分方程式とベクトル場 17

図 1.5: ベクトル場の図，結節点(λ >0とλ <0)

図1.6: ベクトル場の図，鞍点(λ >0とλ <0)

(3)と(5)の場合は安定結節点，(1)と(4)を不安定結節点といいます．(2)では鞍状点といいます．

渦状点や渦心点のベクトル場は，図1.8のように，λ >0のときには，解は原点のまわりを回転しながら遠ざかり，λ <0のときには近づく．λ= 0のときには，図 1.9のように，原点のまわりをぐるぐるまわります．渦状点でも渦心点でもµ >0のときには右まわり，µ <0のときには左まわりになります．

定理 1 行列の固有値は

det(xE−A) = 0

の解である．この方程式を固有方程式という．とくにAが2次正方行列ならば，その固有方程式は

x²−(Tr A)x+ detA= 0

となる．ここで，Eは単位行列，detAはAの行列式，Tr AはAのトレー

(18)

図 1.7: ベクトル場の図，結節点(Jordan標準形，λ >0とλ <0)

図 1.8: 固有値が複素数の場合のベクトル場(λ >0とλ <0)

ス(対角成分の和)を表す．

この定理をみればわかるように実行列であっても固有値は実数とは限りません．

定理 2 2次正方行列に異なる固有値λ, µがあるとしよう．a

¯= (

a₁ a2

) とb

¯= (

b1

b2

)

を対応する固有ベクトルとする．すなわち，a

¯とb

¯は0ベクトルではなく

Aa¯=λa

¯ Ab

¯=µy

¯

をみたすとする．このとき，この2つのベクトルを並べた行列P = (

a₁ b₁ a2 b2

)

は

P⁻¹AP = (

λ 0 0 µ

)

となる．

(19)

1.3. 微分方程式とベクトル場 19

図1.9: 固有値が複素数の場合のベクトル場(λ= 0，µ >0とµ <0) これを行列の対角化といいます．微分方程式

dx¯ dt =Ax

¯ の左からP⁻¹をかけると

dP⁻¹x

dt ¯ =P⁻¹AP P⁻¹x

¯ となりますので，P⁻¹x

¯= y

¯= (

y1

y2

)

とおけば

dy dt¯ =

( λ 0 0 µ

) y

¯

すなわち，対角行列の場合に帰着できます．固有値が複素数になるときには，

固有多項式は実係数の2次方程式ですから，λ±µの形をしています．この場合には行列U =

( 1/√

2 i/√ 2

−i/√ 2 1/√

2 )

とすると

(P U)⁻¹A(P U) =U⁻¹P⁻¹AP U = (

λ −µ

µ λ

)

と表せます．

この他に，固有値がλしかない場合にはP⁻¹AP = (

λ 0 0 λ

)

と対角化できる場合とP⁻¹AP =

( λ 1 0 λ

)

となる場合があり，下の場合がジョルダン標準形です．

一般的な例を考えてみましょう．

(20)

問い 1 ボルテラの微分方程式 dx

dt =x−xy dy

dt =−y+xy の相図を描いてみよ．

図1.10: ボルテラの微分方程式とファン・デル・ポールの方程式

この方程式は2つの特異点(0,0)と(1,1)をもち，(0,0)は鞍状点，(1,1)は渦心点になります．ベクトル場は，図1.10のようになります．このモデルはx をバッタ，yをカマキリと思ったときの生物モデルです．バッタは自分の数に比例して増大しますが，カマキリがいると食われてしまって減少します．一方，カマキリは自分達だけでは友食いをして減少しますが，バッタがいれば増えることができます．このモデルの解はバッタとカマキリが増加減少をくりかえすという自然現象に近いものとなっています．

問い 2 次の相図を求め，解の漸近的挙動を推定せよ．

dx

dt =x²+y²−y dy

dt =x−2xy

Aが0行列のときには，この方法だけでは解の全体像を想像することはできないし，この方法がいつもうまくいくわけではないことには注意しておきましょう．

例 1 ファン・デル・ポールの方程式 dx

dt =y dy

dt =−εy(1−y)²−x

(21)

1.4. 万有引力と電磁力 21 はε= 0のときには，x²+y²の値は解の時間発展とともに変わりません，したがって，円の上を回転する解になりますが，ε >0になると，1つの周期解が存在し，周期解の近くの解はその周期解に近づいていくことが知られています．この周期解を極限周期解とよびます．この極限周期解は上のような不動点の近傍での線形化では求められません．

1.4 万有引力と電磁力

2つの質量M とmの質点があるとき，mにはMから万有引力 F¯=−GM m

r² e

¯ が働くことが知られています．ここでe

¯はMからm方向への単位ベクトルです．したがって，2点間には引き合う力が働くことになります．地球の表面に 1kgの物体をおいたとき，地球の中心に地球の質量M が集中しているとみなせますから，r= 6400kmを代入したときに働く力が重力加速度g= 9.8m/s² になるわけです．いろいろな場所に単位重さををおいて力の働く方向と強さをベクトルで表すと，ベクトル場のような図が得られます．2つ質点がある

図 1.11: 中心力場と2質点の作る重力場

ときには図1.11の右図のようになります．位置の空間は相空間では質点の運動を表すわけではありませんが，どのような動きをするかの示唆は得ることができます．これを重力場といい，G

¯で表すことにします．

電気でもQ, qクーロンの電荷をもつ点の間にはクーロン力(Charles Coulomb, 1736–1806)

F¯ =kQq r²e

¯

が働きます．重力場とは異なるのは電荷には符号があることです．重力場と同様に各所に1クーロンの電荷をおいたときにその力のベクトルを描いた物

(22)

を電場といいます．マイナスの電荷をもつ点が1つだけあるときには，プラスの単位電荷は引きつけられるので，図1.11になり，プラスの電荷をもつときには反発するのでベクトルの方向が逆になります．2つの電荷がプラスの電荷をもつと反発するので図1.12左図になり，2つの電荷がマイナスの電荷

図1.12: 電力場，反発と吸引

をもつときには，図1.11右図になります．2つの電荷がプラスとマイナスの電荷をもつときには，図1.12右図になります．

N極とS極がばらばらに存在することが確かめられたらば，磁力でも同じ式が成り立つはずです．クーロン力や磁力でも電場E

¯や磁場B

¯を考えることができます．電場はそこに1クーロンの電荷を置いたときにその電荷にかかる力を表します．磁場ではそこに単位磁力をおいたときにその点にかかる力と言いたいところですが，単位磁力という概念がありませんので，これでは考えられません．クーロン力と磁力の間には密接な関係があります．電流が流れれば磁場が生じ，逆に磁力が動けば電流が生じるという相互に同じような働きをします．これを実験的に確かめたのがFaraday(Michael Faraday, 1791–1867) で，数式にまとめたのがMaxwell(James Maxwell, 1831–1879)です．速度v で移動するqクーロンの電荷に働く力は ¯

F¯ =q(E

¯+ v

¯×B

¯) で表されます．ここで×は外積を表します．

電荷が移動するということは電流が流れることになります．電荷が連続的な場合には単位長さあたりの電荷密度をρとすると1時刻に単位面積を通過する電荷の総量はρv

¯です．これは単位面積を通過する電流j

¯

になりますので，

電流密度です．これで電流が流れている場合にも力を求めることができます．

上の式から電場と磁場がわかれば力がわかり，それから運動方程式を適用すれば運動が記述できることになります．それで，必要なのは電場と磁場の間に成り立つ式を考えればいいことになります．これらは以下の式にまとめられます．

(23)

1.4. 万有引力と電磁力 23 1. div E

¯ = _ε^ρ

0

2. rot E

¯ =−^∂B

∂t¯ 3. div B

¯ = 0 4. c²rot B

¯ =^∂E

∂t¯ + j

ε¯₀

これが電磁気学のすべてといってもいいそうです．ここでε0は定数です．

∇=





∂

∂x

∂

∂y

∂

∂z





でナブラといいます．·は内積を表し，ナブラを用いると発散 divは

div





f(x, y, z) g(x, y, z) h(x, y, z)



=∇ ·





f(x, y, z) g(x, y, z) h(x, y, z)



= (∂f

∂x +∂g

∂y +∂h

∂z )

(x, y, z)

と表せます．また×は外積を表し，ナブラを用いると回転 rotは

rot





f(x, y, z) g(x, y, z) h(x, y, z)



=∇ ×





f(x, y, z) g(x, y, z) h(x, y, z)



=





∂h

∂y −^∂g_∂z

∂f

∂z−^∂g_∂x

∂g

∂x−^∂f_∂y



(x, y, z)

と表せます．

同じように関数f(x, y, z)について勾配

gradf(x, y, z) =∇f(x, y, z) =





∂f

∂x

∂f

∂y

∂f

∂z



(x, y, z)

はグラディエントと読みます．

発散，回転，勾配と言われるには当然理由があるのですが，ここではそれについては述べません．

上の式の解説を試みましょう．電荷がまったく存在しないところにも電場は存在します．電荷がまったく存在しない領域を囲む領域Dを考えると，D の表面∂Dの上全体を考えると入ってくる電場と出て行く電場は釣り合っているはずです．もし，電荷が中に存在していればその値が求まるでしょう．これを表したのがdiv E

¯ = _ε^ρ

0 です．磁場では単位磁化というのが存在しないので，常に0になります．それを表したのがdiv B

¯ = 0です．電場は磁場に影響し，磁場は電場に影響することを表したのが残りの式です．

具体的には，ガウスの定理より

∫

D

div E

¯dxdydz=

∫

∂D

E¯·n

¯dS

(24)

ここでn

¯は∂Dの微小面積dSの外側向きの単位法線ベクトルを表します．左辺は電磁場の方程式より ∫

D

ρ ε0

dxdydz

ですからD内の全電荷をε0で割ったものです．右辺は電場の領域Dへの出入りを表面全体について足したものを表していることがわかります．

同様にストークスの定理より，2次元の領域Sの縁を左回りにする曲線を

Cとすると ∫

S

rot E

¯·n

¯dS=

∫

C

E¯ ·ds

¯ が成り立ちます．ここでs

¯はCを左回りにまわる単位接線ベクトル表します．

電場の方程式より左辺は ∫

S

−∂B

∂t¯ ·n

¯dS

右辺は電場が境界で接線方向にむいている量の総和です．Sに垂直な磁場の時間変化がCの境界の電場を与えていることになります．同様に

∫

S

c²rotB·n

¯dS=

∫

S

∂E¯

∂t ·n

¯dS+

∫

S

j ε¯0·n

¯dS

です．したがって，Cの境界での磁場はSを垂直に通る電場の変化と電流によって与えられることになります．

重要なことは真空中の静電磁場の場合です．電荷がないのですから 1. div E

¯ = 0 2. rot E

¯ =−^∂B

∂t¯ 3. div B

¯ = 0 4. c²rot B

¯ =^∂E

∂t¯

となります．前の考察から動く電場は磁場を作り，動く磁場は電場を作るので真空中を伝わっていくことになります．これが電磁波です．電磁波の方程式を考えましょう．(4)を微分して

c²rot∂B

∂t¯ = ∂²E

∂t²¯ 左辺は

−c²rot rot E

¯ =−c²grad div E

¯+c²∆E

¯ =c²∆E

¯ を得ます．したがって，

∂²E

∂t¯² =c²∆E

¯

は波動方程式と呼ばれる波の方程式です．もちろん磁場B

¯も同じ方程式をみたします．ここで

∆ = ∂²

∂x² + ∂²

∂y² + ∂²

∂z²

(25)

1.4. 万有引力と電磁力 25 でラプラシアンといいます．ここで用いた

rot rot = grad div−∆ div rot = 0

は定義に基づけば容易に確かめられます．

この方程式の解は速度cをもつ波，例えば

E¯ =





sin(x+ct) sin(y+ct) sin(z+ct)





になります．このことはよく考えてみると不思議です．電磁波は真空中では速度が一定なのです．速度は相対的なものですから，すれ違う電磁波からお互いをみれば2倍の速度になるはずです．この考察の矛盾をなくすために生まれたのが相対性理論です．

(26)

第 2 ^章 Markov chain

2.1 Markov ^連鎖

離散型の確率変数X1, X2, . . .について

P(Xn =xn|X1=x1, . . . , Xn−1=xn−1) =P(Xn=xn|xn−1=xn−1) をみたすときマルコフ連鎖という．条件付き確率を用いて表現するならば，

BnでXnを可測にする最小のσ–algebra，B˜nでX1, . . . , Xnを可測にする最小のσ–algebraとすると，m > nならば

P{ω:X_m(ω) =x_m|B˜n}(ω) =P{ω:X_m(ω) =x_m| Bn}(ω)

が成り立つことである．すなわち，時刻n以降の事象はnでの値にのみより，

それまでの履歴によらないと言うことである．

X1, X2, . . .が独立ならば，ω∈ {ω:X1(ω) =x1, . . . , Xn(ω) =xn}のとき

P{ω:Xm(ω) =xm|B˜n}(ω) = P{ω: Xm(ω) =xm, X1(ω) =x1, . . . , Xn(ω) =xn} P{ω:X₁(ω) =x₁, . . . , X_n(ω) =x_n}

= P{ω:X_m(ω) =x_m}

P{ω:X_m(ω) =x_m| Bn}(ω) = P{ω: Xm(ω) =xm, Xn(ω) =xn} P{ω:X_n(ω) =x_n}

= P{ω:X_m(ω) =x_m} であるから，Markovである．

例 2 X1, X2, . . .を公平な硬貨投げで，

Xn=



 +1 表

−1 裏としよう．このとき，

S_n =X₁+· · ·+X_n

はマルコフ連鎖になる．これを対称random walkという．

確率変数X₁, X₂, . . .が独立同分布なら，条件付き確率を行列で表すことができる．例えば，2つの値a, bしかとらないときには

Π = (

P(X2=a|X1=a) P(X2=b|X1=a) P(X₂=b|X₁=a) P(X₂=b|X₁=b) )

(27)

2.1. Markov連鎖 27 になる．これを推移確率行列という．

π= (P(X₁=a), P(X₁=b)) を初期確率という．また

Π²= (

P(X3=a|X1=a) P(X3=b|X1=a) P(X₃=b|X₁=a) P(X₃=b|X₁=b) )

および

πP = (P(X2=a), P(X2=b)) などをみたすことに注意しよう．

一般の推移確率行列

Π =







p11 p12 · · · p1n

p21 p22 · · · p2n

... ... . .. ... pn1 pn2 · · · pnn







と初期確率

π= (π1, . . . , πn) は

1. pij ≥0 2. ∑n

j=1pij = 1 3. π_i≥0 4. ∑n

i=1πi= 1 をみたす．

∑

i

π_ip_ij=π_j

をすべてのjについてみたす，すなわちπΠ =πをみたすとき，定常であるという．

問題 1 1≤i, j≤nについて，推移確率が

pij =Ciij で与えられるとき，C_iを求めてください．

(28)

解.

∑n j=1

p_ij =c_iin(n+ 1) 2 より

C_i =1

i × 2 n(n+ 1)

□

2.2 ^極限定理

推移確率行列Πについて，Πⁿの(i, j)成分をp⁽ⁿ⁾_i,j で表すことにしよう．これは初めにiにいたときに時刻nにjにいる確率であることは行列のかけ算をしてみればわかるだろう．とくにn= 2ならば

p⁽²⁾_i,j =

∑k l=1

pi,lpl,j

であることからわかる．

推移確率行列Πが非退化であるとは，任意のi, jについて，あるnが存在して，p⁽ⁿ⁾_i,j >0をみたすこととする．非退化でない行列は非退化な行列に分解できることは容易にわかる．

非退化な推移確率行列について，n→ ∞の場合を考えてみよう．

定理 3(Perron–Frobenius) Aを非退化かつ成分が非負な行列とするとき，

その最大固有値は正かつ単純であり，対応する固有ベクトルとして成分がすべて正なものを選ぶことができる．

を用いよう．1が固有値であることは





 1 ... 1





が固有ベクトルであることからわかり，さらに，Πⁿが推移確率行列になっていることから，1が最大固有値であることがわかる．なぜなら，もし絶対値が1以上の固有値があればΠⁿ





 1

... 1





 は発散してしまう．さらに1が最大固有値であることから，定常確率が存在することも示された．

1以外の固有値がすべて絶対値で1より小さいときには，任意の初期状態 πについてπΠⁿは定常確率に収束することがわかり，さらに

1 n

n−1

∑

k=0

e^ikθ→0

(29)

2.3. 吸収壁ランダムウォーク 29 であることを用いると

π1 n

n−1

∑

k=0

Π^k

も定常状態に収束することがわかる．これをエルゴード性という．

2.3 吸収壁ランダムウォーク

自分は1円，相手はN−1円の，2人合わせてN円所持しているとする．

1回あたり1円を賭けてギャンブルをするとき，私が破産する確率を求めてみよう．負ける確率をp，勝つ確率をqとする．

これを解くにはpnで，初めnのところにいたときに破産する確率とする．

欲しいのはp₁だけだが，このように一般の場合を求めるとよい．1回ギャンブルをするとn−1かn+ 1に移動するから

pn=p×pn−1+q×pn+1

が成り立つ．さらに境界条件p₀= 1，p_N = 0が成り立つ．式を書き換えると pn+1=1

qpn−p qpn−1

であるから ( pn+1

p_n )

= (

1/q −p/q

1 0

) ( pn

p_n₋₁ )

ここでP = (

1/q −p/q

1 0

)

の固有値は1と^p_q であるから，対角化する行列を U とおいて，

( α β )

=U⁻¹ (

p1

p0

)

とおくと，

pn=α+β (p

q )n

が成り立つ．境界条件より

α+β= 1, α+β (p

q )N

= 0 であるので，

α= p^N

p^N −q^N, β=− q^N p^N −q^N より，

pn= p^N

p^N −q^N − q^N p^N −q^N

(p q

)n

(30)

を得る．この式はp=q= ¹₂のときには使えない．それには，上の式を整理して

pn =pⁿ(p^N⁻ⁿ−q^N⁻ⁿ) p^N −q^N として，p, q→ ¹₂ ととれば

p_n= 1− n N が導ける．より，正確にはP =

( 2 −1

1 0

)

が，固有値1のJordan標準形であることから

p_n =α+nβ になることに境界条件を用いればよい．

逆に，相手が破産する確率をq_n とおくと，これはpとqを入れ替えて，

N−nをみればよいので

qn= q^N⁻ⁿ(qⁿ−pⁿ) q^N −p^N

このことから， pn+qn = 1，つまり，そのうち，どちらかが破産する確率が1であることがわかる．

では，どのくらいの時間で破産するのだろうか，p_n,kで最初nにいて，k 時刻後に破産する確率とする．上と同様に，

pn,k =p×pn−1,k−1+q×pn+1,k−1 (k≥1) が成り立つ．両辺にt^kをかけて，k= 1, . . .について加えると

∑∞ k=1

pn,kt^k =pt

∑∞ k=0

pn−1,kt^k+qt

∑∞ k=0

pn+1,kt^k

pn(t) =∑_∞

k=0pn,kt^kとおくと，上の式は

pn(t)−pn,0=ptpn−1(t) +qtpn+1(t)

ここで，n̸= 0ならば，p_n,0= 0，とくにp₀(t)≡1かつp_N(t)≡0である．

2.4 反射壁ランダムウォーク

両側の壁で反射して戻るランダムウォークを反射壁ランダムウォークと言います．1つの例をあげておきましょう．

(31)

2.4. 反射壁ランダムウォーク 31

2.4.1 エーレンフェストの壷

実際の物理現象ではさまざまな要因が絡み合い，複雑な様相を呈するためにメカニズムの原理が見つけにくい場合が少なくありません．そこで原理の部分だけを取り出して自然の仕組みを解明してみようというわけです．エーレンフェストの壷は，多数の粒子の運動を記述する統計力学のモデルとして提唱されたわかりやすいモデルです．

2つの壷を用意します．それと玉がN個あるとします．玉には番号がふってあります．これらの玉を2つの壷に分けていれ，1時刻ごとに1からNまでの番号を1つ選んで，その番号の玉を入っている壷から取り出して，他の壷へと移動します．これをひたすらに繰り返すというのがこのモデルです．物理の対象としては，2つの球を細い管でつないだものを用意します．この装置の中に空気の粒子が多数入っていますが，その合計個数をN個とします．

空気の粒子はでたらめに動きまわりますから，ときどき管を通って他の球へと移動していきます．これを長い時間観察したらどんな風になるだろうかというわけです．

たとえ，初めに片方の球が空っぽでも，時間がたてば多い方から少ない方へと移動するチャンスが多いわけですから，長い時間には2つの球に入っている粒子の数(エーレンフェストの壷なら2つの壷に入っている玉の数)は等しくなって来るだろうという予想がたちます．これを数学的にチェックしたいというわけです．両方の球に半分ずつ粒子が入った状態は平衡状態とみなすことができます．このモデルは，同時に平衡状態とはなんだろうという問題の解決にもなっているはずです．

簡単なモデル

多数の粒子のモデルを扱う前に，もっとシンプルなモデルで記号などを準備しておきましょう．2つの状態AとBというのがあるとします．状態Aからは1時刻たつと確率pで状態Bに移動します．ということは確率1−pで状態Aのままでいるわけです．状態Bにいると1時刻後には必ず状態Aに移動するものとしましょう．これは工場の機械の状態をモデル化したものです．状態Aは故障なく機械が働いている状態，状態Bは機械が壊れてしまった状態とすると，正常なときにはある確率pで壊れて状態Bに移動しますが，

壊れたら必ず次の時刻では修理して正常な状態に戻すという風に考えてください．こうしたものは確率オートマトンなどといわれ，コンピュータサイエンスでも用いられる概念です．

このメカニズムは行列

P = (

1−p p

1 0

)

(32)

で表現できます．行列の(1,1)成分P11は状態Aにいるという条件のもとで1 時刻後にも状態Aにいる確率，(1,2)成分P₁₂は状態Aにいるという条件のもとで状態Bに1時刻後に移動する確率，(2,1)成分P21は状態Bにいるという条件のもとで状態Aに1時刻後に移動する確率，最後に(2,2)成分P₂₂ は状態Bにいるという条件のもとで1時刻後にも状態Bにいる確率を与えていますので，この行列を推移確率行列といいます．

P²を考えてみましょう．この(1,1)成分P₁₁² は行列のかけ算からP11P11+ P12P21に等しいわけですが，この初めの項は初め状態Aにいるという条件のもとで，次もAにいてさらにその次もAにいる確率を表しています．第2 項は初め状態Aにいるという条件のもとで，次は状態Bにいって，またその次に状態Aに戻って来る確率を表しています．ということは，P₁₁² は初めA にいるという条件のもとで2時刻後にもAにいる確率を表しています．一般にP_ijⁿについてもおわかりですね．

それでは，今状態Aにいる確率をπ₁，状態Bにいる確率をπ₂としましょう．1時刻後に状態Aにいる確率は今Aにいて，次もAにいる確率と今B にいて次にAに移動する確率の和ですからπ₁P₁₁+π₂P₂₁と表せます．これは行列を用いると(π1, π2)P の第1成分ですし，この横ベクトルの第2成分は同様に1時刻後に状態Bにいる確率に等しくなっています．横ベクトル π= (π1, π2)を初期確率といいます．