数学・物理通信
6巻6号

編集

2016 年 6 月

新関章三・矢野2016 年 6 月 6 日

忠1

目次 (Contents)
1. Jensen の不等式とエントロピー
岩崎正春

2

嶋喜一郎

7

平

11

2. 「光の電磁波説」の誕生

3. ニュートン法の固有値問題への応用
清人

4. ヌメロフの解法—コメント
矢野

忠

20

世戸憲治

21

矢野

忠

28

Masaharu IWASAKI

2

Kiichiro SHIMA

7

5. Bessel 関数と一次元格子振動

6. 編集後記

1. Jensen’s Inequality and Entropy

2. The Birth ofElectromagnetic Theory of Light

3. An Application of Newton’s Method to Eigenvalue Problems
Kiyoto HIRA

11

4. Numerical Solution by Numerov’s Method — A Comment
Tadashi YANO

20

5. Bessel Function and One-Dimensional Lattice Oscillation
Kenji SETO

21

Tadashi YANO

28

6. Editorial Comments

Jensen の不等式とエントロピー
岩崎正春 1Jensen’s Inequality and Entropy
Masaharu Iwasaki2

1

はじめに
エントロピーという概念はもともと物理学から生まれたが，今では情報その他の種々な分野で使われている．

この小論では情報分野で定義されるエントロピーと物理学で使われるエントロピーとの関係について議論する．
その過程において数学で有名な「Jensen の不等式」が登場するが，その初等的な証明についても述べてみたい．

2

情報科学でのエントロピー
はじめに情報科学でのエントロピーの定義からはじめよう．情報科学と言っても広範な分野だが，ここでは『情報理論入門』（ヤグロム，みすず書房）を参考にする．ある試行αは n 個の試行の結果（これを以後「事象」と呼ぶ）を持っているとする：{A1 , A2 , · · · , Ak , · · · , An }それぞれの事象が実現する確率はつぎのようになっている：{p1 , p2 , · · · , pk , · · · , pn }このとき，この試行に対する「あいまいさ」の程度を考えよう．
「あいまいさ」とは実現の不確かさを表すから，log p−1
k に比例するとしてよい．たとえば pk = 1 のときは，事象 Ak は確実に実現するので「あいまいさ」はない．つまり，「あいまいさ」はゼロになる．また，pk が小さくなるにつれ「あいまいさ」は単調に増大するか
ら我々の直観にも合致し，対数を用いる上記の定義が妥当であることが分かる．したがってその期待値として試行全体についての「あいまいさ」を定義する．[1]ϵ(α) = ⟨log p−1 ⟩ = −n∑

pk log pk

(2.1)

k=1

以後，この式を情報理論でのエントロピーと呼ぼう．特別な場合として，事象 A1 が必ず実現する場合，つまり「あいまいさ」が全然ない場合を考えると，p1 = 1, pk = 0 (k ̸= 1) となり，エントロピー ϵ(α) = 0 となる．
つぎに 2 つの連続する試行 α と β について考えよう．試行 β は事象 {B1 , B2 , · · · , Bm } をもち，実現確率が{q1 , q2 , · · · , qm } とすると，ϵ(β) = −m∑

qk log qk

(2.2)

k=1

これらの連続試行 αβ の事象は Ck,l であり，その確率を pk,l とすると，ϵ(αβ) = −n ∑m∑
k=1 l=1

1 高知大学名誉教授
2 miwasaki@cure.ocn.ne.jp

2pk,l log pk,l(2.3)

もし試行 α と β が独立ならば，pk,l = pk ql となるから

ϵ(αβ) = −n∑

pk log pk

m
∑

ql −

l=1

k=1

n
∑
k=1

pk

m
∑

となりエントロピーの「加法性」が成り立っている．ただし，

3

ql log ql = ϵ(α) + ϵ(β)

(2.4)

l=1

∑

k pk =∑

l ql = 1 を用いた．

Jensen の不等式
「あいまいさ」の程度を表すエントロピー ϵ は非負の量である．それでは上限についてはどうなっているのだろうか．それを調べるための道具が以下でとりあげる Jensen 不等式である．
［定理１］（Jensen の不等式）y = f (x) を a から b までの区間の，上に凸関数 3 とする．この区間内に n 個の互いに異なる点 {a ≤ x1 < x2 < · · · < xn ≤ b} をとる．また，和が 1 になる n 個の正数を {p1 , p2 , · · · , pn } と∑nすれば，（ k=1 pk = 1）
( n
)
n
∑
∑
pk f (xk ) < f
pk xk
(3.1)
k=1

k=1

が成立する．ただし，n ≥ 2 とする．
さて，この不等式を証明するため，つぎの補助定理を証明しよう．

[補助定理] 平面上に凸 n 角形 4 があり (n ≥ 3)，各頂点の位置ベクトルを {r 1 , r 2 , · · · , r n } とする．
（図１参照）和が 1 になる n 個の正数 {p1 , p2 , · · · , pn } よりベクトルR≡n∑pk r kk=1を定義すれば，点 R（以後，重心と呼ぶ）はこの n 角形の内部に存在する．

図 1: 凸 n 角形の重心3 関数のグラフの任意の弧 MN が弦 MN の上側にあれば，この関数を上に凸と言う
4 多角形内部の任意の 2 点 M,N を結ぶ線分が多角形内部に含まれるとき凸多角形と言う

3

(3.2)

[証明]

数学的帰納法を用いる．n = 3 の場合明らかに成り立つ．(n − 1) 角形において命題が成立するとする．

R をつぎのように変形する．

1 − pn ∑pk r k + pn r n
1 − pnk=1∑n−1
p′k ≡ pk /(1 − pn ) (1 ≤ k ≤ n − 1) とすると， k=1 p′k = 1 に注意して，
n−1

R=

R = (1 − pn )n−1∑

p′k r k + pn r n ≡ (1 − pn )R′ + pn r n

(3.3)

(3.4)

k=1

帰納法の仮定により，R′ は凸 (n − 1) 角形 {r 1 , r 2 , · · · , r n−1 } の内部，つまり元の凸 n 角形の内部にある．ところで上式は，重心 R が R′ と r n を両端とする線分を pn : (1 − pn ) に分ける内分点であることを意味しているから，R も n 角形の内部に存在することになる．[証明終]
ところで，力学を学んだ学生ならこの補助定理をつぎのように解釈すると分かりやすい．n 角形を各頂点に
質量 {p1 , p2 , · · · , pn } をもつ質点系（剛体）と見なす．このとき点 R は力学での剛体の重心を表している．良く知られているように頂点の一つ r k を支点としてこの剛体を支えると，重心は鉛直線上にあるから R と r kを結ぶ線分はその頂角の内部にある．他の頂点についても同様に考えると，重心は多角形の内部に位置することになる．
つぎに Jensen の不等式の証明に戻ろう．

[定理１の証明]

定理の内容を図２に示す．まず n ≥ 3 の場合を考えよう．凸関数上の n 個の点をベクトル

図 2: 凸関数 f (x) のグラフ

{r 1 , r 2 , · · · , r n } で表し，この n 角形に対して補助定理を適用する．各頂点の座標を r k = (xk , yk ) と置くと，内部の点 R はつぎのように書ける．

R = (Rx , Ry ) =( n∑

pk xk ,

k=1

n
∑

)
pk yk

関数 f (x) が上に凸であるから，Ry < f (Rx ) が成立する．つまり
( n
)
n
∑
∑
pk f (xk ) < f
pk xk
k=1

k=1

4

(3.5)

k=1

(3.6)

となり Jensen の不等式が導かれる．n = 2 の場合も R が，r 1 と r 2 を両端とする線分の内分点になることに注意すると同様に成立する．[証明終]この節の最後に，Jensen の不等式を n 個の点のうち何点かが一致する場合に拡張しておこう．

[定理２]f (x) を a から b までの区間の凸関数とし，この区間内に重複を許して n 個の点 {x1 , x2 , · · · , xn } をとる．また，和が 1 になる n 個の正数を {p1 , p2 , · · · , pn } とすれば，
( n
)
n
∑
∑
pk f (xk ) ≤ f
pk xk
k=1

(3.7)

k=1

ただし，等号は n 個の点がすべて一致する場合に限る．

[証明] n 個の点は m 個の互いの異なる x 座標をもつグループに分かれているとする：{xi,j (1 ≤ i ≤ m, 1 ≤ j ≤ ji )}．
つまり各点は縮退している：x′i ≡ xi,j (1 ≤ j ≤ ji )．
a ≤ x′1 < x′2 < · · · < x′m ≤ b
このとき

(3.8)

∑m

i=1 ji = n を満たす．また，和が 1 になる n 個の正数を pi,j (1 ≤ i ≤ m, 1 ≤ j ≤ ji ) とすれば，[定理２] に現れる級数はつぎのように変形できる．
ji
m ∑
∑

ji
m
m ∑
∑
∑
′
pi,j f (xi,j ) =(pi,j )f (xi ) =p′i f (x′i )i=1 j=1

新しく導入された係数 p′i ≡
ji
m ∑
∑

i=1 j=1∑ji

j=1 pi,j もpi,j f (xi,j ) =i=1 j=1

(3.9)

i=1

∑m

′
i=1 pi = 1 を満たすことに注意して

m
∑

(m
∑

p′i f (x′i ) < f

i=1



)
p′i x′i

=f

ji
m ∑
∑


pi,j xi,j (3.10)

i=1 j=1i=1

中央の不等式には [定理１] の Jensen の不等式を適用した．以上の議論は m ≥ 2 の場合であり，m = 1 の場合
は p′1 = 1 であるから不等式でなく等式になる．[証明終]4

エントロピー増大の法則
この節からこれまでの情報科学のエントロピーを熱平衡状態にある量子多体系に応用してみよう．一定のエ

ネルギー E をもつ孤立した量子多体系を考える．エネルギー E をもつ微視的量子状態を事象 Ak とし，それが実現する確率を xk とする．情報科学のエントロピーは前節の定義により

ϵ(α) = −n∑

xk log(xk )

(4.1)

k=1

ただし，

∑

k xk = 1 である．ここで定理２において f (x) = −x log(x) および pk = 1/n (1 ≤ k ≤ n) と置いた

不等式を用いると，

ϵ(α) =

n
∑
k=1

( ∑n
f (xk ) ≤ nf

k=1 xk

n

)
= nf

( )
1
= log(n)n(4.2)

右辺は定数であるからエントロピーの値には上限があることになる．上限を与える確率分布は xk = 1/n (1 ≤ k ≤ n)であるから α0 と置くと，ϵ(α) ≤ ϵ(α0 ) ≡ log n
5

(4.3)

熱力学のエントロピー増大の法則によると，熱平衡状態においてこの最大値が実現していることになる．最大値を与えるのは，すべての微視的状態の実現する確率が同じになる場合である．つまり統計力学 [2] で言う「等重率の原理」がエントロピー増大の条件より導かれた．しかもエントロピーの表式は，エネルギー E をもつ微視的状態の総数の対数である．したがって式全体にかかる比例定数を除き，統計力学のミクロカノニカル分布のエントロピーに一致する．
つぎに温度 T の熱浴と接しているカノニカル分布のエントロピーについても触れておこう．統計力学 [2] によればこのときの微視的状態の確率は

1 −βEke(4.4)
Zで与えられる．
（エントロピーの加法性を用いてミクロカノニカル分布から導出される）ここで β ≡ T −1 ，Ek =∑n
状態 k のエネルギー，Z ≡ k=1 e−βEk は分配関数である．このときの情報科学のエントロピーは次のように計算される．
nn∑∑
1 −βEkϵ(α) =eβEk + log(Z) = β
Ek xk + log(Z)
(4.5)
Z
xk =

k=1

k=1

右辺第１項はエネルギーの期待値 ⟨E⟩ だから，自由エネルギー F ≡ −T log Z を用いて次の良く知られた関係
式に帰着する．

F = ⟨E⟩ − ϵT

(4.6)

以上述べたことにより，情報科学で定義されたエントロピーは物理学のエントロピーと一致し，かつそれを内包するものであると言える．

参考文献[1] アー・エム・ヤグロム，イー・エム・ヤグロム，「情報理論入門」（みすず書房，1958）.
[2] 久保亮五，「統計力学」（共立出版，1952）.

6

「光の電磁波説」の誕生

The Birth ofElectromagnetic Theory of Light
嶋 喜一郎 ∗
Kiichiro Shima†

1

はじめに

√
1856 年のウェーバーとコールラウシュの実験について, 後藤尚久 [1] は誘電率 ε と透磁率 µ に着目して, 1/ εµの値が求められた実験と考えていた. しかし,『電磁気の単位はこうして作られた』[2] で取り上げられている実験をみると, 電磁単位と静電単位の比が問題になっているだけであった. 前者では電磁単位 qm と静電単位 qsの比についての言及はなく, 後者では誘電率 ε と透磁率 µ についてはふれられていなかった. 2 つはどのように関係しているのだろうか.
√
後藤 [1] は 1/ εµ を静電単位系（ε = 1）と電磁単位系（µ = 1）の値を換算するときの係数（速度の次元を√もつ）として取りあげていたが, この 1/ εµ は多くの本で見られるように電磁波の速度を表わしている 1 . マク
スウェルは電磁波の存在を導くときウェーバーとコールラウシュの実験を参考にしている [2]. 一方に電磁波の√速度を表わす式がある. 他方に, その値と関係がありそうな実験結果がある. 1/ εµ と qm /qs の関係は, このときのモデルになっているとみることができる 2 .
この小論は「光の電磁波説」の誕生を歴史的にではなく論理的に把握することをめざす. マクスウェルはど√のようにして 1/ εµ を c（光速度）と同定したのか. 「光の電磁波説」の誕生を次の 3 つの公式と 1 つの定義
式を基礎にして考察する [2][3]3 .

qq ′, （電荷に対するクーロンの法則）
r2mm′
F2 = k2 2 , （磁荷に対するクーロンの法則）
r2mI
F3 = k3, （ビオ・サバールの法則）
rdq
I=, （電流と電荷の関係・定義式）
dt

F1 = k1

(1.1)
(1.2)
(1.3)
(1.4)

ここで，F, q, m, I, r, t は, それぞれ力，電気量，磁気量，電流，距離，時間を示す. また k1 , k2 , k3 は比例定数である．

∗ 元岐阜大学工学部† kiichiro-s@cg8.so-net.ne.jp

√
k/4πµ
√
である [2]. もちろん, その値は 1/ εµ と同じ c である。
2 1856 年の実験は「同じ電気量」を静電単位と電磁単位で測定したものである. 測定値の比 Q /Qsm = c で実験を示す表現が主流である [2]. しかし, ここでは電磁単位と静電単位の電気量の大きさの比 qm /qs = c で実験を表現する.
3 1880 年ころまでは電磁現象は近接作用ではなく遠隔作用によって理解されていた. これらの式を取りあげるのは電磁波以前の時代背景を考慮するからである. また電磁気の単位系を区別する比例定数 k1 , k2 , k3 が過不足なくそろっているからである.
1 ε は正確にいえば真空の誘電率 ε , µ は真空の透磁率 µ である. この式はマクスウェルのものではない. マクスウェルの式は
0
0

7

2

比例定数の制約関係（1.1）(1.2)（1.3）（1.4）（３つの法則と 1 つの定義式）を連立させ, 力が均衡しているときの条件をみていく. ３つの法則を F32 /(F1 × F2 ) = 4 となるように調整する. 例えば F1 = 1[dyn], F2 = 1[dyn], F3 = 2[dyn] とする. これは係数を簡単にし, 関係を簡潔にするためである.

22 m2 I 2 1 r2
1 r2
k32・ 2 × ・ ′ × ・r
k1 qq
k2 mm′=4

比例定数を前にもってきて, 両辺を 4 で割る.

k32
m2 I 2 r 2
r2× 2 × ′×
k1 k2
r
qq
mm′

=1

k32
m2 (dq)2 r2
r2
× 2
× ′×
2
k1 k2
r (dt)qqmm′=1

ここに I = dq/dt を代入して,対応する物理量（電荷の電気量 q, 磁荷の磁気量 m, 1 組の r2 ） で割り算して( dr )2
k32
×
=1
k1 k2
dt
したがって,( dr )2
k1 k2
=
= v2
k32
dt

k1 k2 /k32 = v 2 の関係は電磁気の単位系（静電単位系・電磁単位系・ガウス単位系・MKSA 単位系）を構成する基礎になっている 4 . k1 k2 /k32 = v 2 の平方根をとって,
1/2 1/2

k1 k2
k3

3

=v

(2.1)

マクスウェルの思考実験
マクスウェルは n を 1 つの電磁単位の中にある静電単位の数と仮定して（n = qm /qs )5 , この n が速度になることを洞察している. マクスウェルはクーロンの法則による斥力（左辺）とアンペールの法則による引力（右辺）が均衡する条件から導いている.

CC ′ n2 t2= CC ′r2
等号が成り立つのは r = nt, すなわち n = r/t のときである（Maxwell,“A Treatise on Electricity andMagnetism” Art768 [3]）. いいかえれば, 電磁単位と静電単位の比は速度になる.
qm
=v
qs

(3.1)

（2.1）と (3.1) はどちらも力が均衡する条件の中で得られたものである. （2.1）と (3.1) の v は同じであると考えられる.
1/2 1/2

k1 k2
k3

=

qm
qs

(3.2)

「比例定数の制約関係」と「電磁単位と静電単位の比」は対応しているのである.
4 . 静電単位系は k = 1, k = 1, k = v 2 である. 電磁単位系は k = 1, k = 1, k = v 2 である. ガウス単位系は k = 1, k = 1, k = v
3
1
2
3
2
1
1
2
3
である. また MKSA 単位系は電磁単位系に由来する.
5 “n is the number of electrostatic units in one electromagnetic unit”（Maxwell,“A Treatise on Electricity and Magnetism”
Art768 [3]）84電磁波の速度
1/2 1/2

k1 k2
k3

=

qm
qs

比例定数の制約関係は同時に電磁波の速度の式を表わす. ここで, ビオ・サバールの法則の比例定数（連結因子）k3 を k3 = 1 とする. また電荷に対するクーロンの法則の比例定数 k1 , 磁荷に対するクーロンの法則の比
例定数 k2 を, それぞれ k1 = 1/ε, k2 = 1/µ とおくと,
1/2 1/2

k1 k2
k3

1
=√
εµ

(4.1)

したがって,

√

1
qm
=
εµ
qs

(4.2)

である. 電磁波の速度の大きさを決定するには qm /qs の測定値が必要になる. 他方, 1856 年のウェーバーとコー
ルラウシュの実験は qm /qs = c を示している.
したがって,

1
√ =c
εµ

(4.3)

√
マクスウェルは, 電磁単位と静電単位の比 qm /qs を媒介して, 電磁波の速度 1/ εµ は光速度 c であると同定したのである.
（了）（2016.5.7）

参考文献[1] 後藤尚久,『電磁波とはなにか』（講談社, 1984）[2] 木幡重雄,『電磁気の単位はこうして作られた』（工学社, 2003）[3] 「電磁気学の単位系が難しい理由」，ＦＮの高校物理，（2010）http://fnorio.com/0096Electromagnetic_
unit_system1/Electromagnetic_unit_system1.html

補足
嶋喜一郎氏の『「光の電磁波説」の誕生』の補足である．
同じ電気量を静電単位系と電磁単位系で測定する．このときのそれぞれ単位系での電気量の単位を小文字の

qs , qm で表す．またそれらの電気量の単位で測った測定量を大文字の Qs , Qm で表す．そうすると同じ電気量をそれぞれの単位系で測ったのだから，

Qs qs = Qm qmが成り立つ．

9(4.4)

そのとき当然のことながら，(4.4) から

qm
Qs
=
=c
qs
Qm

(4.5)

成り立つ．
Qs
そのことから脚注 2) で「測定値の比 Q= c で実験を示す表現が主流であるが，ここでは電磁単位系と静電m単位系での電気の単位の大きさの qqms = c で実験を表現する」と述べられている．このことは上に述べたことから当然のことであるが，一言注意をしておく．
たとえば，吉田卯三郎著『物理学』下（三省堂，1941 年初版，1963 年四訂版) 184-185 に Weber と Kohlrauschの実験について書いてあり，電流の静電単位系と電磁単位系の電流の両単位の比が真空中の光速の比に等しいとあるが，電磁単位の大きさ
静電単位の大きさ
is
=
im

c=

(4.6)

とある．これは電気量ではなくて，電流について考えているが，測定量の数値の比と単位の大きさの比は逆数関係になることを示している．
（編集者（矢野忠）記）(2016.5.12)

10

ニュートン法の固有値問題への応用
平

清人 1

An Application of Newton’s Method to Eigenvalue Problems
Kiyoto Hira 2

1

はじめに
世戸憲治氏と中西襄氏が，本誌において，量子力学の一般的教科書では馴染みのない Rosen-Morse のポテン

シャル問題を超幾何関数を用いて，解析的に解いていることは物理数学的に興味深い [1] [2]．1968 年に石津武
彦氏が Woods-Saxon 型のポテンシャル問題を S 状態に限ってだが，解析的に解かれていることも同様に興味深い [3]．これらの問題は量子力学のどの教科書のはじめに出てくる 1 次元や 3 次元の井戸型ポテンシャル問題
をより一般化した問題と捉えることができる．
昨今では高精度な性能をもつ PC は大変に身近になった．従って，解析的に解くのが難しいこれらの物理数学的な問題でも，容易に数値計算できるようになった．その結果，問題解決のための数値解析的手法の修得は量子力学だけでなく物理学一般の学習にも不可欠になった．そのような一例として，Rosen-Morse ポテンシャ
ルの束縛状態の数値解法はこのような現状からみて大いに意味があるだろう．

2

Rosen-Morse ポテンシャルの束縛状態の解析的解法の要約
理論的詳細は彼らの論文を参照していただくことにして，本節では数値計算の方法を示し，その結果を検討するために必要な中身だけを要約する．
世戸氏と中西氏は，次の形で書かれる無次元化した 1 次元のシュレーディンガー方程式の微分方程式を[ 2]d+ E − V (x) Ψ = 0
(2.1)
dx2
を考察した．この時に，彼らは２つの型のポテンシャルを考察した．１つは，V (x) = V1 tanh(x) −V2
cosh2 (x)(−∞ < x < ∞)(2.2)

で与えられる．これを本論では，Rosen-Morse potential type 1 と呼ぶ．もう一つは，V (x) = V1 tan(x) +V2
cos2 (x)(−π/2 < x < π/2)(2.3)

で与えられる．これを本論では，Rosen-Morse potential type 2 と呼ぶ．それぞれの具体的考察対象となるシュ
レーディンガーの微分方程式は次のようになる．
[ 2
]
V2
d
+
E
−
V
tanh(x)+Ψ = 0,
1
dx2
cosh2 (x)1 元広島県立高校教諭
2 da43827@pb4.so-net.ne.jp

11

(2.4)

さらにもうひとつは

[

]
d2
V2
+
E
−
V
tan(x)−
Ψ = 0.
1
dx2
cos2 (x)(2.5)(2.4) 式の微分方程式の境界条件は，束縛条件から，
Ψ(x) → 0forx → ±∞(2.6)

for

x → ±π/2(2.7)(2.5) 式の微分方程式の境界条件は，束縛条件から，
Ψ(x) → 0
を満足しなければならない．
彼らは巧みに量子化されたエネルギー固有値 En とその En の固有関数 Ψn (x) を求めることに成功した．まず (2.4) 式に関して，2 つの定数 µ と ν を次のように定義する．

µ = V1 /2,√
ν = V2 + 1/4 − 1/2.
このとき，エネルギー固有値 En は

En = −µ2− (ν − n)2(ν − n)2である．ただし，n は

n<ν−

√

µ

(2.8)
(2.9)

(2.10)(2.11)

を満たす自然数である．また，そのときの En の固有関数 Ψn (x) は，規格化定数を除いて，µx

Ψn (x) = 2−(ν−n) e− ν−n cosh −(ν−n) (x)F (−n, 2ν + 1 − n,µ1
+ ν − n + 1; (1 − tanh(x)))ν−n2である．このときの関数 u(x) = F (α, β, γ; z) は次の超幾何関数の微分方程式[]d2dx(1 − x) 2 + {γ − (α + β + 1)x}− αβ u(x) = 0dxdx(2.12)(2.13)

を満たす x = 0 で正則な解である．具体的に表すと，F (α, β, γ; z) = 1 +αβα(α + 1)β(β + 1) 2
x+x + ···1!γ
2!γ(γ + 1)

(2.14)である．
次に (2.5) 式に関して，2 つの定数 µ と ν を次のように定義する．

µ = V1 /2,
ν = 1/2 +√(2.15)
V2 + 1/4.

(2.16)

このとき，エネルギー固有値 En は

En = (ν + n)2 −µ2(ν + n)2(n = 0, 1, 2, · · · ).

(2.17)である．また，そのときの固有関数 Ψn (x) は，規格化定数を除いて，

Ψn (x) = e−pn x cos(ν+n) (x)F (−n, 1 − 2ν − n, 1 − (ν + n) + pn i; 1/2(1 + i tan (x))
である．ここで，i は虚数単位であり，また

pn =

µ
ν+1

12

(2.18)(2.19)

と定義される．
本論では，このポテンシャル問題を数値解析的に扱う．ポテンシャルの強さを与えるパラメター V1 と V2に具体的な数値を与えて，微分方程式 (2.4) 式と (2.5) 式を数値的に積分してその解を求めると同時に，エルギー固有値 En とその固有関数 Ψn (x) (n=1, 2, · · · ) を併せて数値的に決定する．微分方程式を積分するために Numerov 法を使い，方程式の解を求める数値的解法として超有名な Newton 法を併せて使うことにより，エネルギー固有値 En とその固有関数 Ψn (x)(n = 1, 2, · · · ) を同時に求めることができる [4]．それらの具体的処方箋は次節で要約して説明をする．

エネルギー固有値とその固有関数の数値解法

3

我々が数値的に積分する微分方程式は次の形をしている．

d2 y= g(x)y,
dx2

(3.1)

ただし，g(x) = −E + V (x) である．(3.1) 式の解 y(x) を x = a の周りでテイラー展開する．

y(a + h) + y(a − h) − 2y(a) = y (2) (a)h2 +22 (4)y (a)h4 + y (6) (a)h6 + · · · ,
4!
6!

(3.2)

ただし，y の肩の添え数は微分の階数を表す．同じように，

y (2) (a + h) + y (2) (a − h) − 2y (2) (a) = y (4) (a)h2 +2 (6)y (a)h4 + · · ·
4!

(3.3)

を得る．(3.1) を使うと，g(a + h)y(a + h) + g(a − h)y(a − h) − 2g(a)y(a) = y (4) (a)h2 +2 (6)y (a)h4 + · · ·
4!

(3.4)

2 2
2 (4)
を得る．そこで，(3.2) − 4!
h × (3.4) から 4!
y (a)h4 を消去すると，

(

)
(
)
(
)
h2
10h2
h2
g(a) y(a)1 − g(a + h) y(a + h) + 1 − g(a − h) y(a − h) − 2 +
12
12
12
1 (6)
y (a)h6 + · · · . (3.5)=−
240

ここで，h の 6 次以上を無視すると,
(
)
(
)
(
)
h2
h2
10h2
1 − g(a + h) y(a + h) + 1 − g(a − h) y(a − h) − 2 +g(a) y(a) = 0
12
12
12

(3.6)

を得る．
この式を使うと，はじめに y(a − h) と y(a) を与えると y(a + h) を決めることができる．この手順を続けると次々と y(a + 2h), y(a + 3h), y(a + 4h), · · · が決まる．こうして (3.1) 式を数値積分して微分方程式を解く方
法が Numerov 法である．この方法を利用するときには初めに 2 つの y の値を与えることが必要である．
いま (3.1) の微分方程式の x の定義域を A ≤ x ≤ B とする．この区間 [A, B] を N 等分に分割し各分割区間の幅を h とすると，

B−A,N
ここで，各分割点 xi は xi = A + ih (i = 0, 1, 2, · · · , N ) と表せる．
h=

2

2

そこで，vi = V (xi ), ci = h12 g(xi ) = h12 (−E + vi ), yi = y(xi )(3.7)(i = 0, 1, · · · , N ) とおけば，(3.6) 式から，次の漸化式を得る．

(1 − ci+1 )yi+1 + (1 − ci−1 )yi−1 − (2 + 10ci )yi = 0.
13

(3.8)

y0 , y1 の初期値を与えて，この漸化式の解 yi (i = 0, 1, 2, · · · , N ) を数値的に求めれば，{(xi , yi )|(i = 0, 1, 2, · · · , N )}(3.9)

が，定義域 A ≤ x ≤ B における微分方程式 (3.1) の数値解を与える．
次にエネルギー固有値を決定する．そのために，(2.6) または (2.7) の境界条件を使って考える．(2.4) 式では V (x) の変化が無視できるほど十分に大きく A と B を取ることにする．そのとき，十分に大きな x = B では (2.6) の境界条件下では y(B, E) ≈ 0 となる．また (2.5) 式では，[A, B] = [−π/2, π/2] と取る．(2.6) の境界条件下では y(B, E) = 0 となる．結局，y(B, E) ＝ 0 となる E を数値的に決定すれば，その時の E がエネ
ルギー固有値 E になる．この方程式の解をもとめるために Newton 法を利用する．
その計算手順を説明する．まず，積分区間 [A, B] を十分に大きい自然数 N で分割する．そして，(3.8) 式の漸化式から，y(B, E) ＝ yN を数値的に求める．そこで，f (E) = yN とおく．f (E) = 0 の解 E は, Newton 法によれば，次の漸化式

En+1 = En + δEn = En −f (En )f ′ (En )(n = 0, 1, 2, · · · ),(3.10)を次々に繰り返していくと En (n = 0, 1, 2…) が解 E に収束することが期待できる．ただし，E0 は f (E0 ) ≈ 0′を満たす試行値である．その際に f (En ) は次のように計算する．
まず，微小量 ∆En を次のように定義する．n = 0 のときは，に取り，n = 1, 2, · · · のときは，′と取る．そして，f (En )∆E0 = 0.001E0(3.11))(f (En−1 )∆En = 0.001δEn−1 = 0.001 − ′f (En−1 )(3.12)(n = 0, 1, 2 · · · ) は，f (En + ∆En ) と f (En − ∆En ) を数値計算し，′f (En ) =f (En + ∆En ) − f (En − ∆En )2∆En(3.13)

nから数値計算する．この手順を繰り返して En (n = 0, 1, 2, · · · ) を計算していく．仮にこの手順を | δE
En | < 1.0e−10を満たす n で適当にやめれば，その時の En が固有値を与え，また (3.9) がその固有関数を与える．

Newton 法の漸化式を使うときに，エネルギー固有値の近似値である初項の E0 が必要になる．それをどのように見当するか考える．そのためには，エネルギー固有値 E が，f (E) = yN とおいたとき，f (E) = 0 を満たす性質を利用すれば見当がつく．そのために，E を適当な範囲 (Ea ≤ E ≤ Eb ) に設定する．そこで，[Ea , Eb ]を小分割して，各分割点 Ek (k = 0, 1, 2, · · · ) ごとに yN (Ek ) を (3.8) 式使って計算をする．そして各 Ek に対して log10 |yN (Ek )| の値をプロットしていけば，エネルギー固有値の近傍では yN (Ek ) の値が 0 に近付くので,log10 |yN (Ek )| の値が急激な変化をするので，エネルギー固有値の近似値の見当がつく．

4

具体的な数値計算例
はじめに (2.4) のシュレーディンガー方程式の微分方程式を取り上げる．V1 ＝ 2.0 と V2 ＝ 8.75 とする．そして，積分区間は [A, B] ＝ [−15.0, 3.8] とし，分割数は N = 500 とする．そのときは (2.8) と (2.9) 式より µ＝ 1.0, ν ＝ 2.5 となる．条件式 (2.11) より，n < 1.5 となる．この時のエルギー固有値 En は，n ＝ 0 と n ＝ 1に限られる．(2.10) 式から E0 ＝ − 6.41 と E1 ＝ − 2.69444444444444 となる．またそれらの固有関数は (2.12)
より，

Ψ0 (x) = N0 2−2.5 e−2x/5 cosh −2.5 (x),()
15
Ψ1 (x) = N1 2−1.5 e−2x/3 cosh −1.5 (x) 1 − (1 − tanh(x) ,
19
14

(4.1)
(4.2)

ただし，N0 と N1 は規格化定数であり数値積分を実行して決める．これらの理論値と我々の計算値を比較す
る．その厳密解との比較は我々の数値解法の信頼性を確かめるためには大変に有意義である．
まず，Rosen-Morse Potential type 1 の E に対する log10 |yN (E)| のグラフを示す (図.1)．このグラフから，
4.5

log10|yN(E)|43.5

log10|yN(E)|

3

2.5

2

1.5

1

0.5

0

-0.5
-10

-9

-8

-7

-6

-5

-4

-3

-2

E

図 1: E に対する log10 |yN (E)| のグラフ. 急峻な極小を取る E の値がエネルギー固有値に対応している ．

E = −6.4 と E = −2.7 あたりに鋭い極小値の谷が見られる．これらの谷の底の E の値が，固有値の近似値である．これら 2 つの値の近傍値を Newton 法の漸化式の初項にとって固有値と固有関数を数値計算した結果を示す (表．1)．
表 1: Newton 法 [Rosen-Morse potential type 1]

E0 （試行値）-6.5

E1-6.415698696654360

E2-6.410024686305030

E3-6.410000163502920

E4-6.410000163046980

収束値-6.410000163046980

理論値

-6.41

E0 （試行値）-2.4

E1-2.689102091455590

E2-2.694440976936360

E3-2.694437981284390

E4-2.694437981283400

収束値-2.694437981283400

理論値-2.694444444444444

この表の結果から，数列 Ej (j = 0, 1, 2, 3, · · · ) は，E0 ＝ − 6.5 では，4 回の繰り返しで収束する．また E0＝ − 2.4 では 4 回の繰り返しで収束する．ともに収束は速いと云える．収束値と理論値を比較してみるとほぼ
有効数字 6 桁の精度である．精度は，積分区間 [A, B] と分割数 N に依存する．
つぎに併せて得られる固有関数の数値解と正確な解析解を重ねて視覚的に示す (図.2)．

15

1

0.8

calculated with E0=-6.4
calculated with E0=-2.7
Eq.(4.1)
Eq.(4.2)

0.6

0.4

Ψ

0.2

0

-0.2

-0.4

-0.6
-12

-10

-8

-6

-4

-2

0

2

4

x

図 2: 横軸は x 座標を表し，縦軸は固有関数値を表す．連続曲線は数値解である．厳密解は x の離散値での値をプロット

正確な解析解は，離散的な x 点で値をクーロス印でプロットした．数値解と解析解の２つを比べると非常に高い精度で一致していることが分かる．計算の結果から固有関数の値は有効数字が 6 桁の精度で合致している
ことが分かる．勿論， 積分区間 [A, B] を広げたり，また区間 [A, B] の分割数 N を大きく取ることよって精度は更に上がる可能性がある．しかし，B を大きく取りすぎると波動関数が発散するというこの微分方程式の固有の性質が現れることにも注意しておく．Numerov 法の漸化式 (3.8) を計算するために初期値として，

y0 = exp(
y1 = exp(√√|E + V1 |x0 ),(4.3)|E + V1 |x1 )(4.4)とした．ただし，x0 ＝ − 15.0 である．漸化式で求めた波動関数は規格化されていないので数値積分を実行して規格化定数を求め規格化した．また，数値計算の結果の固有関数のグラフについて，節がない関数は収束値

E ＝ − 6.410000163047 に対応している．また，節が一つの関数は収束値 E ＝ − 2.69443798128342 に対応している．
つぎに (2.5) のシュレーディンガー方程式の微分方程式を取り上げよう．そして，数値計算の結果と理論的結果を比較する．V1 = 2.0, V2 = 0.75 を採用する．そして，積分区間 [A, B] ＝ [−π/2 + 1.0e − 6, π

/2 − 1.0e − 6]，そして分割数は N = 1200 とした．± 1.0e-6 はポテンシャルの特異点を避けるために導入した．すると (2.15) 式と (2.16) 式から µ＝ 1.0, ν ＝ 1.5 となる．(2.3) 式は x ＝ ± π/2 で ∞ になる井戸型ポテ
ンシャルである．エネルギー固有値 En (n = 0, 1, 2, · · · ) は,(2.17) 式から E0 = 1.80555555555556, E1 = 6.09 ,そして E2 = 12.1683673469388 であることがわかる．またそのときの固有関数は (2.18) より，

Ψ0 (x) = N0 e−2x/3 cos 1.5 (x),(4.5)

Ψ1 (x) = N1 e−2x/5 cos 1.5 (x) (0.4 cos(x) + 1.5 sin(x)) ,
(
)
4
8
457
Ψ2 (x) = N2 e−2x/7 cos 1.5 (x) (3.75 − ) cos(2x) − sin(2x) −
,
49
7
196

(4.6)

ただし，N0 と N1 そして N2 は規格化定数であり，数値積分を実行して求める．
つぎに，E に対する,log10 |yN (E)| のグラフを示す (図.3)．

16

(4.7)

5

log10|yN(E)|4log10|yN(E)|

3

2

1

0

-1

-2

-3
0

2

4

6

8

10

12

14

E

図 3: E に対する log10 |yN (E)| のグラフ. 急峻な極小を取る E の値がエネルギー固有値に対応している ．

このグラフから，E = 1.8 と E = 6.0，そして E = 12.1 あたりに鋭い極小値の谷が見られる．これらの谷の底の E の値が，固有値の近似値である．これら 3 つの値の近傍値を Newton 法の漸化式の初項にとって固有値と固有関数を数値計算した結果を示す (表．2)．
表 2: Newton 法 [Rosen-Morse potential type 2]

E0 （試行値）

1.4

E1

1.730223018280786

E2

1.802379039328814

E3

1.805544320478877

E4

1.805550226321220

E51.805550226349218

収束値

1.805550226349218

理論値

1.805555555555560

E0 （試行値）

6.09

E1

6.089986523500320

E2

6.089986523537520

収束値

6.089986523537520

理論値

6.090000000000000

E0 （試行値）

12.2

E1

12.16822631365359

E2

12.168337379130056

E3

12.168337380460052

E4

12.168337380460848

収束値

12.168337380460848

理論値

12.168367346938800

この表の結果から数列 Ej (j = 0, 1, 2, 3, · · · ) は試行値 E0 ＝ 1.4 では 5 回の繰り返しで収束する．また試行17

値 E0 ＝ 6.09 では 2 回の繰り返しで収束する．そして, 試行値 E0 ＝ 12.2 では 4 回の繰り返しで収束する．正確な解析解と比べると，有効数字が 6 桁の精度をもつことが分かる．さらに，我々の方法の信頼性を確かめる
ために，固有値と併せて得られる固有関数の数値解と厳密解を重ねて視覚的に示す (図.4)．厳密解 (4.5) 式と(4.6) 式，そして (4.7) 式は，離散的な x 点で固有関数値をクーロス印でプロットした．節の数が２つの関数は収束値 E ＝ 1.8055502263492178 に対応する．節の数が３つの関数は 収束値 E ＝ 6.08998652353752 に対応する．そして，節の数が４つの関数は 収束値 E ＝ 12.168337380460848 に対応する．視覚的にも判断される通り，数値解と厳密解の 2 つを比べると非常に高い精度で一致していることが分かる．数値解と厳密解とは計算の結果から有効数字が 6 桁の精度で合致している．
1
calculated with E0=1.4
calculated with E0=6.09
calculated with E0=12.2
Eq.(4.5)
Eq.(4.6)
Eq.(4.7)

0.8
0.6
0.4
0.2

Ψ

0
-0.2
-0.4
-0.6
-0.8
-1
-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

x

図 4: 縦軸の Ψ は計算値と理論値に対応した固有関数を表す

5

おわりに以上見てきたように，ヌメロフ法と ニュートン法の計算アルゴリズムは非常に簡単かつ明解で十分な精度を持ち，初めの試行値を適切に選択すれば，後は自動的に最適値を探してくれる．また視覚的にも満足のいく信
頼性を持っている．多くのポテンシャルの固有値問題に手軽に使えて，演習もできるから量子力学を初めて学ぶ者にとっても教育的であろう．
より具体的に興味をもたれる読者は，筆者の Website を参照していただきたい [5]．プログラムは JAVA 言語で書かれていて Java Applet が使用されているのでセキュリティーの問題で，ブラウザーが動作しない場合がある．その時は問題の解決の仕方が Web 上で公開されているので，検索されればすぐに解決すると思う．

参考文献[1] 世戸憲治，量子力学における Rosen-Morse ポテンシャル問題， 物理 · 数学通信 5 巻 4 号 (2015) 9-19
[2] 世戸憲治，中西襄，量子力学における Rosen-Morse ポテンシャル問題 (2)， 物理 · 数学通信 5 巻 6 号 (2015)16-26[3] T. Ishidzu, Analytical Solution for S-State with the Woods-Saxon Potentials, Prog.Theor. Phys. 40(1968) 796-807
18

[4] M. Atarashi, K. Hira and H. Narumi, On the Kaon-Nucleus Optical Potential at Low Energy, Prog.
Theor. Phys. 60 (1978) 209-219
平 清人， パソコンで量子力学を解こう， パリティ 8 巻 4 号 (1993) 61-63[5] 平 清人,

http://www013.upp.so-net.ne.jp/pathintegral/Demonstration/index.html

19ヌメロフの解法—コメント
Numerical Solution by Numerov’s Method — A Comennt
忠1

矢野Tadashi YANO2
平氏のヌメロフの解法の説明に問題はないのだが，とうとつにテイラー展開が現れてくるような気がするのでちょっとコメントをしたい．
平氏にしたがって y (n) で n 階の微係数を表す記号をつかう．いま関数 y(x) を x = a で数値微分することを考える．

y(a + h) − y(a)hで１階の数値微分を定義する．同様に x = a − h での数値微分はy (1) (a) =(1)y(a) − y(a − h)h(2)y (1) (a) − y (1) (a − h)h(3)y (1) (a − h) =である．このとき２階微分を

y (2) (a) =で定義すれば，y(a + h) + y(a − h) − 2y(a)h2となる．これが (3.2) の左辺で y(a + h) + y(a − h) − 2y(a) を考えた理由であろう．このとき (4) からy (2) (a) =(4)y(a + h) + y(a − h) − 2y(a) = h2 y (2) (a)(5)

が近似的に成り立つが，テイラー展開を考えれば，右辺には高次の微分項が現れる．それが (3.2) である．後は平氏の説明の通りである．

(2016.5.9)

1 元愛媛大学工学部
2 yanotad@earth.ocn.ne.jp20Bessel 関数と一次元格子振動世戸 憲治 ∗

Bessel Function and One-Dimensional Lattice Oscillation
Kenji Seto∗

1 はじめに一次元格子振動は，通常，三角関数を用いて解かれる問題であるが，ここでは，Bessel 関数を用いて解くとどうなるかということに挑戦してみる．

2 方程式の導入

図1

1 次元格子モデル

図 1 に示すように，質量 m の錘と，ばね定数 k のばねを交互に接続した系を考える．ここでは，これら錘とばねは無限個存在するものと仮定し，全体として半無限の形状をなすものとする．錘には左から右に向かって番号を付け，n = 1, 2, · · · とする．振動しているときの時刻 t における n 番目の錘の変位を Un (t) とし，この錘の運動方程式を考える．この錘には，ばねによって，右から k(Un+1 − Un ) の力を受け，左から −k(Un − Un−1 ) の力を受けるので，その運動方程式は，md2
Un = k(Un+1 − 2Un + Un−1 ),dt2

n = 1, 2, · · ·(2.1)となる．この式で，n = 1 としたときの U0 は一番左にあるばねの左端の位置を表し，この位置は強制的に動かされるものとする．より具体的には，時間 t の関数として，

U0 (t) = F (t)(2.2)と与えられるものとする．これが境界条件となる．

3 Bessel 関数による解法
Bessel 関数 Jν (z) の微分公式は，
1
Jν′ = − (Jν+1 − Jν−1 )2∗ 北海学園大学名誉教授

E-mail: seto@pony.ocn.ne.jp

21

(3.1)

と書ける．ここで，プライムは z 微分を表す．これを 2 度繰り返すと，

Jν′′ =

1(Jν+2 − 2Jν + Jν−2 )
4

(3.2)

となる．この式は方程式 (2.1) を解くために使える．そこで，(2.1) 式の Un (t) を，Un (t) = J2n (ωt)
とおいて方程式に代入すると，

√
ω=2

(3.3)

k def
= ω0
m

(3.4)

であれば満たされる．したがって，J2n (ω0 t) は解となる．
以下，いくつかの初期条件，境界条件を与えたときの解を求めてみよう．

[ U0 (t) = ℓJα (ω0 t) とした場合]
α を任意の無次元で非負の定数, ℓ を長さの次元を持つ定数として，境界条件，初期条件を，時間微分を上付ドットで表すことにして，

U0 (t) = ℓJα (ω0 t),

Un (0) = 0,U̇n (0) = 0,n = 1, 2, · · ·(3.5)としたとき，その解はあきらかに，Un (t) = ℓJ2n+α (ω0 t)

(3.6)

となる．
より一般的な結果をだすために，つぎの積分公式を引用する *1 ．

∫ ∞
0

(
)
2 sin (ν − ν ′ )π/2dz=.
Jν (z)Jν ′ (z)zπ(ν 2 − ν ′2 )(3.7)

この式の証明は，Jν (z), Jν ′ (z) が満たす Bessel の微分方程式を

[ d2
dz

+
2

(
1 d
ν 2 )]+ 1 − 2 Jν (z) = 0,
z dz
z

[ d2
dz

+
2

(
1 d
ν ′2 )]+ 1 − 2 Jν ′ (z) = 0
z dz
z

(3.8)

と書いておき，この第 1 式に zJν ′ (z) を，また，第 2 式に zJν (z) を掛けて，辺々を引き算すると，

d [ ( dJνdJν ′ )]Jν Jν ′z Jν ′− Jν
= (ν 2 − ν ′2 )
dz
dz
dz
z

(3.9)

となる．これを z について 0 から ∞ まで積分し，∫ ∞Jν (z)Jν ′ (z)0[ ( dJ
1dJν ′ )]∞
dz
ν
′
= 2
z
J
−
J
ν
ν
z
ν − ν ′2dz
dz
0

(3.10)となる．あとは，(3.1) の微分公式と，Bessel 関数の z が十分大きいときの漸近形 *2√Jν (z) ≃(2
2ν + 1 )cos z −ππz4を用いる．さらに，三角関数の積和公式を用いて計算していくと，結果の右辺の式が得られる．
*1 「数学公式 3」(岩波全書) P.193. 上から 6 番目の公式*2 同書，P.154. 上から 4 番目の公式

22(3.11)

(3.7) 式において，ν − ν ′ がゼロでない偶数となるときはこの右辺がゼロとなる．また，ν − ν ′ = 0 のときは，分母分子がゼロとなるので，ν → ν ′ の極限をとることで，右辺 = 1/(2ν) となる．したがって，r, r ′ を非負整数，α を任意の非負実数としたとき，∫ ∞J2r+α (z)J2r′ +α (z)
0

dz
1
=
δr,r′
z
2(2r + α)(3.12)という直交式が得られる．

[U0 (t) = ℓ(ω0 t)α とした場合]
境界条件，初期条件を

U0 (t) = ℓ(ω0 t)α ,Un (0) = 0,

とおく．このときの解を，Un (t) = ℓ∞∑U̇n (0) = 0,n = 1, 2, · · ·(3.13)Cr J2(n+r)+α (ω0 t)(3.14)

r=0とする．ここに，Cr は r に依存する定数とする．これに境界条件である (3.13) の第１式を適用すると∞∑Cr J2r+α (ω0 t) = (ω0 t)α(3.15)

r=0となる．この式から Cr を求めるため，両辺に J2r′ +α (ω0 t)/t を掛けてから，t で 0 から ∞ まで積分すると，
∞
∑

∫ ∞
Cr

J2r+α (ω0 t)J2r′ +α (ω0 t)
0

r=0

dt
=
t

∫ ∞

(ω0 t)α J2r′ +α (ω0 t)

0

dt
t

(3.16)となる．この左辺には (3.12) の直交式が使え，右辺は Weber の積分 *3∫ ∞dz
2α−1 Γ ((ν + α)/2)=zΓ ((ν − α)/2 + 1)

(3.17)2α (2r + α)Γ (r + α)r!

(3.18)z α Jν (z)0が使える．結果は，

Cr =と求められる．これを (3.14) 式に代入し，Un (t) = 2α ℓ∞∑
(2r + α)Γ (r + α)r=0r!

J2(n+r)+α (ω0 t)(3.19)と，Bessel 関数の重ね合せの形で，解が得られたことになる．この式で．特に．n = 0 としたときは，

U0 (t) = 2α ℓ∞∑
(2r + α)Γ (r + α)

r!

r=0

J2r+α (ω0 t)(3.20)

J2r+α (z)(3.21)となるが，これは，z α を Bessel 関数で展開した式 *4z α = 2α∞∑
(2r + α)Γ (r + α)r=0r!

*3 同書，P.192. 一番下の式*4 同書，P.211. 下から 4 番目の式

23

を用いると，(3.20) 式右辺は ℓ(ω0 t)α となって，確かに境界条件の (3.13) 第 1 式を満たしている．

[U0 (t) = ℓ sin (ωa t) の場合]
境界条件，初期条件を

U0 (t) = ℓ sin(ωa t),

Un (0) = 0,n = 1, 2, · · ·U̇n (0) = 0,(3.22)とした場合を考える．これは，ばねの左端を振動数 ωa で強制的に振動させた場合である．その解をUn (t) = ℓ∞∑

Cr J2(n+r)+1 (ω0 t)(3.23)Cr J2r+1 (ω0 t) = sin(ωa t)(3.24)

r=0とする．これに境界条件を適用すると，∞∑
r=0となる．これには，公式 *5

∫ ∞
0


[
]
1

 sin ν sin−1 (a/b) ,dz  ν
sin(az)Jν (bz)
=
)ν
(
b
1
z

√
 sin(νπ/2),ν
a + a2 − b20<a≤b(3.25)
0<b≤aが使える．(3.24) 式の両辺に J2r′ +1 (ω0 t)/t を掛け，t について 0 から ∞ まで積分する．このとき，項別積分が許されるものとし *6 ，この公式，および，(3.12) の直交式を利用すると，


[
]
−1

2
sin
(2r
+
1)
sin
(ω
/ω
)
,
a
0


)2r+1
(
Cr =
ω0
r


√
,
2(−1)

ωa + ωa2 − ω02ωa ≤ ω0(3.26)

ω0 ≤ ωaと Cr が求められる．これを (3.23) 式に代入すると，[] sin (2r + 1) sin−1 (ωa /ω0 ) ()
Un (t) = 2ℓJ2(n+r)+1 (ω0 t),
2r+1
ω
r


√ 0

r=0 (−1)ωa + ωa2 − ω02


∞ 

∑


ωa ≤ ω0
ω0 ≤ ωa(3.27)と求められたことになる．
この式で，n = 0 としたときは，境界条件である (3.22) の第 1 式 U0 (t) = ℓ sin(ωa t) が成立するはずである．
すなわち，

[
] 

−1
sin
(2r
+
1)
sin
(ω
/ω
)
a
0


∞ 

∑
(
)
2
J2r+1 (ω0 t) = sin(ωa t),
2r+1
ω
r


√ 0

r=0 (−1)
22ωa + ωa − ω0ωa ≤ ω0

(3.28)
ω0 ≤ ωa*5 同書，P.203. 上から 5 番目の式
*6 結果論から言うと，ω0 < ωa の場合はこの項別積分が許されない．したがって，この場合は間違えた結果がでてしまうが，ここでは敢えて，そのままにしておき，この場合については，後で述べることにする．なお，“Higher Transcendental Functions”, Bateman
Manuscript Project, Vol.2, P.64 〜 P.65 にかけて，この (3.12) 式の直交性を用いた任意関数の展開式がでているが，そこでは，“We drive formally the expansion , , ,” と，非常に曖昧模糊とした書き方をしており，項別積分に関しては，“see Hardy” としか書かれていない．これは，この変換式に関し，一般的な確固とした理論が存在しないことを意味するものと理解される．

24

が成立しなければならない．このうち，ωa ≤ ω0 の第 1 式の方で，r = 0, 1, 2 の場合について，三角関数の倍角公式を用い，a = ωa /ω0 とおくことにして，


a,
(
) 
−1
sin (2r + 1) sin (a) = −4a3 + 3a,
16a5 − 20a3 + 5a,

r=0
r=1
r=2

(3.29)となるので，これと Bessel 関数を級数展開した，Jν (z) =∞∑k=0( z )ν+2k(−1)kk!Γ (k + ν + 1) 2

(3.30)

11(ωa t)3 + (ωa t)5 − · · ·3!
5!

(3.31)

を用いて (3.28) の第 1 式左辺を展開すると，左辺 = ωa t −となる．これは確かに sin(ωa t) を展開した式と一致する．
これをよりスマートに処理するために，公式 *7

sin(az) = 2∞∑(−1)r T2r+1 (a)J2r+1 (z)(3.32)

r=0

を引用する．ここに，T2r+1 (a) は奇数次の Tchebycheﬀ 多項式である．この式を証明するには，Bessel 関数の母関数表示 *8

exp

[z (
2

t−

∞
∑
1 )]
=
tk Jk (z)t(3.33)k=−∞から出発するのが最も簡単である．ここで，J−k (z) = (−1)k Jk (z) を用い，さらに，k が偶数の部分と奇数の部分に分けると，

exp

[z (
2

t−

∞
∞
∞
∑
∑
∑
(k
)
( 2r −2r )( 2r+1 −(2r+1) )1 )]= J0 (z)+t +(−t)−k Jk (z) = J0 (z)+t +tJ2r (z)+t−tJ2r+1 (z)
t
r=1
r=0
k=1
(3.34)となる．ここで，t = eiθ とおくと，

eiz sin θ = J0 (z) + 2∞∑

cos(2rθ)J2r (z) + 2i

r=1

∞
∑

(
)
sin (2r + 1)θ J2r+1 (z)(3.35)

r=0となる．この式の実部，虚部をとると，∞∑()
cos z sin θ = J0 (z) + 2
cos(2rθ)J2r (z),

∞
∑
(
)
(
)
sin z sin θ = 2
sin (2r + 1)θ J2r+1 (z)r=1(3.36)

r=0となる．ここで必要なのはこの第 2 式の方で，a = sin θ とおくと，

sin(az) = 2

∞
∑

(
)
sin (2r + 1) sin−1 (a) J2r+1 (z)(3.37)

r=0*7 同書，P.89 には何故か a = π/2 の場合しか載っていないが，「数学辞典」(岩波書店) の付録公式 20 には一般の場合が載っている．
*8 同書，P.149

25

となる．ここで，次数が奇数のときの Tchebycheﬀ 多項式 *9()
T2r+1 (a) = (−1)r sin (2r + 1) sin−1 (a)(3.38)

を利用すると，所望の (3.32) 式が得られる．この (3.32) 式で a = ωa /ω0 , z = ω0 t とおいたときは，

sin(ωa t) = 2∞∑(−1)r T2r+1 (ωa /ω0 )J2r+1 (ω0 t)(3.39)

r=0となって，これは，(3.28) の第 1 式にほかならない．また，このときの解 Un (t) は Tchebycheﬀ 多項式を用いて，Un (t) = 2ℓ∞∑(−1)r T2r+1 (ωa /ω0 )J2(n+r)+1 (ω0 t)(3.40)

r=0と書かれる．なお，この式で，ωa = ω0 すなわち，強制振動の ωa を固有振動の ω0 としたときは，Tchebycheﬀ
多項式は次数に関係なく，変数が 1 のときは，Tn (1) = 1 となるので，∞∑
Un (t) = 2ℓ(−1)r J2(n+r)+1 (ω0 t)(3.41)

r=0と式が簡単化される．ただし，このときは，共振現象が起こり，t が十分に大きいとき，Un (t) ∼= 2ℓ(−1)n

∞
∑
r=0

√

2
cos(ω0 t − 43 π) → ±∞πω0 t

(3.42)と発散してしまう．

(3.28) の第 2 式，ω0 ≤ ωa の場合はどうなるであろうか．もし，これが，第 1 式の ωa ≤ ω0 の領域のものがそのまま解析接続されるのであれば，(3.28) の第 2 式もそのまま成立するはずである．しかし，この (3.28) 式の第

1 式と第 2 式は解析接続では繋がっていない．また，実際に，この ω0 ≤ ωa の場合を r の冪が小さいところで，級数展開しても，sin(ωa t) を展開したものとはならない．これは，(3.24) 式に J2r′ +1 (ω0 t)/t を掛けて t 積分をしたときに，t について一様収束していないために，項別積分が許されない式になっているためと考えられる．
そこで，この ω0 ≤ ωa の場合は別に解を作り直す必要がでてくる．(3.34) 式のところで，t = ieθ とおくと，

eiz cosh θ = J0 (z) + 2

∞
∞
∑
∑
(
)
(−1)r cosh(2rθ)J2r (z) + 2i(−1)r cosh (2r + 1)θ J2r+1 (z)r=1(3.43)

r=0となり，この式の実部，虚部をとると，∞∑()
cos z cosh θ = J0 (z) + 2(−1)r cosh(2rθ)J2r (z),

∞
∑
(
)
(
)
sin z cosh θ = 2(−1)r cosh (2r + 1)θ J2r+1 (z)

r=1

r=0

(3.44)となる．ここで必要なのはこの第 2 式の方で，a = cosh θ とおくと，

sin(az) = 2∞∑()(−1)r cosh (2r + 1) cosh−1 (a) J2r+1 (z)(3.45)

r=0となる．この式は，|z| ≤ 1 における Tchebycheﬀ 多項式の定義 *10

Tn (z) =

√
√
)n (
)n ]
1 [(
z + i 1 − z2 + z − i 1 − z22

*9 同書，P.88. 下から 6 番目の式
*10 同書，P.88, 上から 3 番目の式

26

(3.46)

を |z| > 1 の領域に解析接続したTn (z) =において，cosh−1 (z) = log(z +

√
√
)n (
)n ]
1 [(
z + z2 − 1 + z − z2 − 1
2(3.47)()
Tn (z) = cosh n cosh−1 (z)(3.48)√
z 2 − 1) を用いると，と書けるので，(3.45) 式は，sin(az) = 2∞∑(−1)r T2r+1 (a)J2r+1 (z)(3.49)

r=0となる．これは形の上では (3.32) 式とまったく同じものであるが，|a| ≥ 1 の領域に拡張したものである．この式で，前と同じく，a = ωa /ω0 , z = ω0 t とおくと，(3.39) 式が ω0 ≤ ωa の場合でも成立することになる．
以上のことをまとめると，境界条件，初期条件 (3.22) 式を満たす方程式 (2.1) の解は，ω0 , ωa の大小とは無関
係に，Un (t) = 2ℓ∞∑(−1)r T2r+1 (ωa /ω0 )J2(n+r)+1 (ω0 t)(3.50)

r=0と求められたことになる．
この結果から，t が十分に大きいときの Un (t) を Bessel 関数の漸近形 (3.11) 式で見積もると，Un (t) ∼= 2ℓ(−1)n√2cos(ω0 t − 43 π)
T2r+1 (ωa /ω0 )
πω
t
0
r=0

∞
∑

(3.51)となる．ここで，ω0 < ωa の場合は，Tchebycheﬀ 多項式 (3.48) から，T2r+1 (ωa /ω0 ) > 1 なので，この和は，一様収束はおろか，一般に発散してしまう．これが，(3.24) 式から (3.26) 式に至る過程で，ω0 < ωa の場合に項別積分が許されずに，正しい結果が得られなかった原因となっている．

4 おわりに今回ほど，計算に苦労したことはない．(3.12) 式の直交性を用いて任意関数を展開したときは，正しい結果が出る場合と，間違えた結果が出る場合があるということである．これは，通常の Fourier 級数，Fourier 変換，Fourier-Bessel 変換ではなかったことで，初めは，そのことに気付かずに悩んでしまった．例によって中西先生に相談し，ようやくその理由が見えてきた．これは項別積分ができなくなることがその理由とのことであった．私は，初め，これは，Bessel 関数 J2r+1 (z) が完全性を持たないためと考え，その方向での計算を進めてみたが，これは無駄な結果に終わってしまった．それにつけても，Bessel 関数の奥深さと恐ろしさを存分に味わってしまった．

[ 謝辞 ]
今回も，京都大学名誉教授の中西襄先生にこの原稿を見ていただきました．「おわりに」のところで述べたように，初めは分からないところだらけで，先生とメールを何度も交換しながらようやく完成までこぎつけることができました．その意味で，この論文は先生なくして書くことができなかったものです．先生に，心から深謝いたし
ます．

27

編集後記
3 月に続いて今回も連続して 3 つの号を発行することになった．いつも毎号なかなか多種多様な内容の論文が掲載される．人が違えば考えることもちがうことを表しているということであろう．
数学と物理とはいっても結構な広い範囲をカバーしている．エントロピーは最初は熱力学の分野の概念だったのかもしれないが，いまではもっと広い意味に使われている．詳しく岩崎さんの論文を読んでみたい．
熊本地方の地震はもちろんそこに住んでいる人に大きな被害と影響を与えているが，それだけではな
く，日本の地震学にもいままで経験しなかったような新しい経験を与えているように思われる．
この号を発行しようとして，世戸さんから投稿された論文はどうなったかというメールをいただいてようやく彼の論文を忘れていたことに気がついた．ページ数も余裕があったので今号に掲載できたが，どうも編集者失格である．
こういうことが起こるのはもちろん望ましくはないが，多量の投稿をあつかっているので，論文の著者からの自分の論文がどうなったかという問い合わせの編集者へのメールがときには必要になっている．
編集者へのご協力をお願いしたい．
（矢野28

忠）