/tex/num/recc1.tex

漸化式あれこれ
整数又は自然数 n に対してある点列 xn を考えた時、小さな n に対する xn の知識を用いて順次大きな n に対する xn の値を計算する手順を書き下す事が出来る場合がある。
点列の値そのものに興味がある場合に、この手順を漸化式と呼び、xn 又はその和

P

n xn

が収束すると仮定して、収束する値を計算する手順を逐次近似と名付けて区別するのが一般的な言葉の使用法である。しかし点列の個別の値に興味を持つか又は収束値に興味を持つかは、その時々で変わっても良い。面倒だから、ここでは両者を区別せずに、漸化式という言
葉で括ってしまう。ここでは、漸化式に纏わる数値計算の幾らかの話題を提供する。
ある数学者は自然数を以下の様に定義した。
１) １は自然数である
２) 自然数に１を加えたものは自然数である
これで無限個の自然数が出来上がった。
簡単な方程式の解法
ある種の方程式は、以下の様に変形できる。

x = f (x)利用できるならば、ニュートン法を使えば良いが、ニュートン法は微分 f 0 (x) を計算せねばならないので、手間がかかる。そこで、手抜き計算として、以下の方法を推奨する。先ず、仮の解 x0 を持って来る。そして、x1 = f (x0 ) を計算する。|x1 − x0 | が小さければ収束した。収束していなければ、x0 を x1 で置き換えて、計算を続行する。
この方法は、収束の保証は無いが、先ず簡単に試みるべき方程式の解法であるから、覚えておいて損は無い手法の一つである。

(x、y) 平面上に y1 = x と y2 = f (x) というグラフを描いてみると、交点があるか無いか
はすぐにわかるだろう。これで収束の様子も簡単に調べられる。
収束しなければ、方程式を変形して別の表現に変換してみるの良いし、変数変換をしてみるのも良い。ニュートン法等の他の手法に切り替えるのも良いだろう。
正の実数 A の平方根の計算にこの手法を直接適用しても、無駄である。x0 × x1 = A という二つの数を行き来するだけである。しかし、0 < x0 < x1 ならば、x0 <(x0 + x1 )/2 が√√A < x1 という関係にある事は明らかだから、

A の改良された近似になっている事は明らかだろう。この改良法をニュートン法との関係をしらべてみよ。
プランクの輻射の公式を用いて、輻射が最大となる波長の温度依存性を示すウイーンの変移則を導いてみよ。
最後に数値計算が必要である。多分、x = ~ ω/kT と置くと、次の様な式だったと思う。

x = 3(1 − e−x )1ニュートン法にも言及しておこう。F (x) = 0 という方程式を解く時に、微分係数 F 0 (x)が利用でき、初期値 x0 が与えれているならば、x1 = x0 − F (x0 )/F 0 (x0 ) により解を改良する。修正量 |F (x0 )/F 0 (x0 )| が小さくなれば解けた事にし、解けていなければ、x1 を x0 として修正を繰り返す。
ニュートン法は収束域で使用すれば、２次収束するからかなり有効な手法である。収束域に如何にしてたどり着くか？この問題に最も有効な助言は、グラフを描く事である。これは、誰にでも言えるが。
。
次の助言は、解く事は２の次にして解の存在範囲を確定する事に最初の努力をせよ。これは通常は、F (x) の符号が変化する範囲 [xl , xr ] を検出せよという事である。勿論、y = F (x)が x 軸に接するという例外もあるが。
。この時は修正量が発散するが。。
次の助言は、修正量 δ = |F (x0 )/F 0 (x0 )| が解の存在範囲幅 xr − xl よりも充分小さい時には、ニュートン法を適用せずに２分法を利用する、即ち x1 = (xr − xl )/2 を採用し、解に近付くはずの x を採用せずに、区間を狭くする事である。この理由は、区間内で関数値の変化が激しすぎる為に ニュートン法では正解に遅々として近付かないからである。
級数展開の例
２項漸化式の例として、三角関数の級数展開を取り上げよう。

sin x =Xsn =X (−)n x2n+1nn(2n + 1)!!

s1 = x であり、sn+1 = −x2 sn /(2n)(2n + 1) である。
この級数は絶対収束であるから、原理的には任意の x に対して利用できる。実用的には、ある n = N という項 sN があり、そこで絶対値が最大になり、その両側で絶対値は小さくなるという現象が経験される。| sin(x)| ≤ 1 という制限があるので、各項の絶対値が１よりも大きな部分は、和をとると相殺されているはずである。有限桁の計算では、この相殺された桁分だけは有効数字の桁数が少なくり、有効数字の桁数が減る。従って、この漸化式 (展開) には許容誤差をきめると有効範囲が設定されることになる。
この計算法は、ホーナー法に持ち込めるだろう。
指数関数の計算

exp(x) はおなじみの関数であり、グラフを描く事を想定し、x 軸に沿って等間隔に xn = n∆
をとり、ここでの関数値 en ≡ exp (n∆) を計算したいとする。この時、次の漸化式が成り立つ。

en+1 = e1 × en
従って、e0 = 1, e1 = exp (∆) が与えられていると、原理的には、何処まででも計算を実行出来る。指数関数の直接計算は、計算コストがかなりかかるので不利な場合がある。しかし
ここで例示した場合では少しの記憶域と掛け算さえ実行すれば良いので、実用的にもかなり
2

有利な場合がある。
三角関数
例えば台形公式を用いて角度積分をしたいとする。ある角度範囲で、sin (x), cos (x) を等間
隔の x に対して、計算する必要がある。次の公式を高校時代に習った。

sin(x + ∆) = sin(x) cos(∆) + cos(x) sin(∆)
cos(x + ∆) = cos(x) cos(∆) − sin(x) sin(∆)
従って、角度範囲が０度から始まるならば、次の漸化式を書き下す事が出来る。先ず、

s0 = 0, c0 = 1, s1 = sin(∆), c1 = cos(∆) を準備しておく。
sn+1 = sn c1 + cn s1 ,
又は、

Ã

sn+1
cn+1

!

Ã
=

cn+1 = cn c1 − sn s1

c 1 s1−s1 c1

!Ã

sn
cn

!

下の書き方だと、新しいベクトルは古いベクトルを直交変換すれば良いと教えている。この方法で、任意の sin(n∆) と cos(n∆) が計算出来る。
漸化式は、(新しい量) ＝ (古い量)(に何か変換を施す) という形で書かれる場合が多い。

( 古い量) → (新しい量) という変換を時間発展と捉える事が出来る。この時間発展に対して、ある種の保存則が存在する場合もある。例えばエネルギーや角運動量が運動の恒量になっているかもしれない。漸化式が厳密に計算されるならば、この保存則は満足されるのだ
ろうが、計算の常として誤差を含み得るから、時間発展を追跡している内に保存則は満足されなくなる恐れもある。
ここで提示した三角関数の計算は、直交変換で実現されているから、sin2 (x) + cos2 (x) = 1という保存則を常に満足しながら計算される。この意味で、安定な漸化式である。
次のシンプソンの公式を用いて数値積分を、実行する場合を想定する。

Z b
f (x)dx =ah[f (a) + {4f (a + h) + 2f (a + 2h) + · · · + +4f (b − 3h) + 2f (b − 2h) + 4f (b − h)} + f (b)]3

この場合、係数４と２が交互に繰り返すので、次の様な漸化式が考えられる。最初 W = 4 として作業を開始
し、(new)W = 6 − (old)W とすると、４、２、４、２ · · · というトグル動作をする。両端とそれ以外とを別途に計算し、最後に h/3 を掛ければ良い。
台形公式は、簡単な積分公式であるから精度が悪いと思い込んでいる人もいるだろう。そういう人には、

Euler-Maclaurin の和公式を調べる事を勧める。積分の上下端での被積分関数の (高階) 微分の値次第では、台形公式が非常に高精度を与える場合がある。

三項漸化式
３項漸化式の代表として、素人数学者に人気があるというフィボナッチの数列を取り上3

げよう。元々は、兎の数の増え方として、フィボナッチが書いた数学の教科書の問題として登場するという事であるが、階段の上がり方の場合の数として登場する場合もある。

f0 = 0, f1 = 1 を初期値とし、次の漸化式を満足する数列である。
fn+1 = fn + fn−1
この漸化式から一般項を書き下す手法は簡単である。fn = a r n と仮定して、漸化式に代入すると r に関する次の２次方程式をえる。

r2 − r − 1 = 0
この２次方程式の解を r± = 1 ±√

5/2 と書くと、一般解は二つの定数 A, B を用いて

nnfn = A r++ B r−と書かれ、n = 0, 1 での値を代入して定数をきめると、

(Ã

fn =

√ !n )√ !n Ã
√
1− 5
1+ 5
−
/ 5
2
2

漸化式の定義から明らかであるが、二つの定数が漸化式と初期条件を確定するのに必要である。
次の三項漸化式

a fn+1 + b fn + c fn−1 = 0
は、以下のような二項漸化式に変形できる。

a (fn+1 − p fn ) = q (fn − p fn−1 )ここで、p, q は２次方程式を解けば確定する。従って、gn ≡ fn − p fn−1 とおくと、

gn+1 = (q/a) gn = · · · = (q/a)n g1 と書き下せる。
この式を、相続く n に対して書き下すと、
fn+1 − p fn = Rn g1
fn − p fn−1 = Rn−1 g1···
f3 − p f2 = R2 g1
f2 − p f1 = R1 g1
ここで、R ≡ q/a とおいた。上から２番目に式には p、次の式には p2 を掛け、全部の式を加えると、

fn+1 − pn f1 =R
(Rn − pn ) g1
R4 − pこのようにして、定数係数の３項漸化式は簡単に一般項が書き下せる。係数が定数でない場
合にもこの手法は利用できるが、場合に依っては作業が煩雑かも知れない。
三項漸化式、

a fn+1 + b fn + c fn−1 = 0
を fn で割ると、比 rn ≡ (fn /fn+1 ) に関する２項漸化式になる。この２項漸化式を繰り返して利用すると、以下の様な連分数表現が可能である。

rn = d +

e
rn+1

=d+

e
d+

=d+

e

rn+2

e
d+

e

d+

e
..

.

一般には、d や e は n に依存する。連分数の収束する方向に利用すると、面白い世界が広がる。興味があれば、僕のホームページにある連分数の話題も読んでみて下さい。例えば、(旧主要項) ＋ (補正項) ＝ (新主要項) という漸化式に持ち込むと安定した連分数の数値計算ができる。
この様に、連続する３個のパラメータを有する量の間の関係式として表現される例は沢山ある。キーワードとしては、直交多項式の基本漸化式という概念がある。物理の世界では、統一的な話題よりも個別的な話題の方が好まれるから、調和振動子の問題に対する固有関数
としてのエルミートの多項式を例にとろう。
岩波書店発行の数学公式集 (森口・宇田川・一松著) には色々な定義がある旨書かれている。そこで、Schiff の量子力学から引用しよう。シュレーディンガー方程式は、

Hn00 − 2ξ Hn0 + 2 n Hn = 0
であり、Hn (ξ) = (−1)n exp(ξ 2 )dnexp(−ξ 2 )dξ nと表せ、H0 = 1, H1 = 2ξ, H2 = 4 ξ 2 − 2 と書かれている。
更に、次の漸化式が成立するとも記されている。

Hn+1 = 2 ξ Hn − 2 n Hn−1
これで、任意のエルミート多項式が計算できる。気が向いたら、級数展開を利用する方法と、漸化式を利用する方法を利用して計算速度、精度、簡便性等を少し大きめの n に対して
比較してみよ。漸化式の有効性が体感出来るだろう。
三項漸化式の例は非常に多いが、散乱問題の量子力学にしばしば登場する、球ベッセル関
数 jn (x) を取り上げてみよう。数学公式集を見ると、以下の初期値と漸化式が与えられている。

j0 (x) =sin x − x cos x
sin x, j1 (x) =xx2
5jn−1 (x) + jn+1 (x) =2n + 1jn (x)x

一見すると、これで情報は必要充分であるように見える。ところが、数値的にはこれだけの情報では不十分である。良く知られているように、遠心力があるために 0 < x < n では

jn (x) ∝ xn /(2n + 1)!! であるから、小さな x では n と共に関数値が小さくなる。従って、n の小さい方から大きな n へ漸化式を用いて計算すると、有効数字がどんどん減って行く。
更に悪い事が起こる。実は、球ノイマン関数も同じ漸化式を満足し、こちらは n が大きくな
ると関数値が大きくなる。従って、幾らかでも誤差が混じっている j0 (x), j1 (x) を用いて nの増す方向に漸化式を用いると、この誤差は急激に増大して真の値を圧倒してしまう。
球ベッセル関数を、上に与えた漸化式を用いて計算する場合には、以下の手法が採用される。N >> n なる N を想定し、jN +1 (x) = 0, jN (x) = ² を初期値として、漸化式を大きなN から n を下げながら使用する。n = 0 にたどりついたら、j0 (x) = sin(x)/x を用いて全体を規格化する。
この場合の様に、漸化式はどちら向きに使用するのが安全であるかという問題がある。
ある種の経験的な記述が Abramowitz と Stegun が編集した Handbook of Mathematical

Functions with Formulas, Graphs, and Mathmatical Tables という本の XIII に書いてある。この本は、数学の公式集としても非常に参考になる。
関数値が大きくなる方向に、漸化式を使用するというのが、一般的な指導原理である。
規格化には、母関数 (今の球ベッセル関数ならば Rayleigh 展開を思いだすとよい) を利用する場合もある。
ある場合には、二つのパラメータがあり、これらのパラメータの間の関係式として複数の漸化式が与えられている場合がある。この例として、Legendre の陪関数 P`, m (x) を取り上げてみよう。この場合も、複数の定義があるようであるが、上記岩波の数学公式集に登場する、 通常の場合を取り上げよう。第３巻 １２６ー１２７ページには色々な漸化式が与えられている。特に m = 0 の場合を先ず考える。即ち Legendre の多項式 P` (x) である。相隣る ` に関しては、以下の関係式がある。

(` + 1) P`+1 (x) = (2` + 1) x P` (x) − ` P`−1 (x)この漸化式と P0 (x) = 1, P1 (x) = x という知識があれば、これで充分である。
次に m を１だけ増やす事を考える。P0,1 (x) = 0 であるから、P1,1 (x) を計算出来れば良い。既に P0,0 (x), P1,0 (x) の知識があるならば、以下の漸化式が利用できる。

(` − m)P`,m (x) − (` + m)P`−1,m (x) +√1 − x2 P`, m+1 (x) = 0

ここでは、複数の漸化式からどの漸化式を利用するのが得策かという選別眼が要求される様になる。
指導原理は、` と m で作られた２次元の格子点を思い描く事である。即ち二つの ` と一つの m、又は 二つの m と一つの ` に関する漸化式を見比べて、どの二つが既知でどれが未知かという事を考えれば良い。

6Legendre 関数が登場したところで、直交多項式と数値積分との関連する話題を提供しておこう。先ず直交多項式の定義から述べよう。定義域 [α, β] で正の関数 w(x) があり、積分Rβw(x)dx は有界であるとする。無限区間ならば、積分が収束するような振舞が w(x) にはα要求される。

n 次多項式 pn (x) =Pn
k=0

ak xk が n = 0, 1, · · · に対して与えられていて、これらが以下の式を満足する時、規格化された直交多項式と言う。
次の準備として、Lagrange 補間式を与えておこう。N 個の点 (xi , f (xi )), (i = 1, 2, · · · N )を通る (N − 1) 次多項式を書き下せという問に対する答えである。頭ごなしであるが、Lk (x) = Πi6=k (x − xi )/(xk − xi ) は N − 1 次多項式であり、i 6= k ならば Lk (xi ) = 0 であり、Lk (xk ) = 1 である事は、具体的に代入すると了解出来る。即ち、Lk (xi ) = δi k と書ける。従って、次式で与えられる L(x) は N 個の点 (xi , f (xi )) を通る N − 1 次式である。
X
L(x) =f (xk ) Lk (x)k

これも、具体的に代入すればすぐに確認出来るだろう。
次に、以下の数値積分を考える。

Z β

I=

w(x) f (x) dx
αf (x) は必要な回数だけ微分可能であると仮定しておく。
数値積分の公式は、積分を以下の様に近似する事であるとする。

I∼

N
X

ρi f (xi )i=1

即ち、積分区間内のある点 xi での関数値 f (xi ) に荷重 ρi を掛けて加えると、積分の近似値を得る。確かに台形公式やシンプソンの公式はこの範疇に入る。
さて、次の問題を考えよう。
出来るだけ少ない回数 f (x) を評価して、出来るだけ正確な積分値を得る為には、(xi , ρi )の組をどのように選べば良いか？
ガウスは以下の様にしてこの問に答えた。pN (x) は直交多項式の内、次数が N のものを取り出したとする。f (x) を 2N − 1 次式だと考えて、f (x) = pN (x) q(x) + r(x) と割算をする。q(x) と r(x) は高々 N − 1 次多項式である。ここで、xi として pN (x) = 0 の N 個の独立根を採用すると、f (xi ) = r(xi ) が成立する。即ち、一般的には 2N − 1 次の多項式の値が、この特殊な点を採用すると N − 1 次多項式で与えられる！勿論、この事が言えるのは、他の点を選択した場合には有り得ない。結局

Z β
w(x) f (x) ∼

I=
α

7

X
i

ρi r(xi )となった。r(x) を近似する N − 1 次式として L(x) を採用すると、

I∼

Z βX
N
α

Z β
r(xi )

w(x) Li (x)dx =

X

α

i=1

とするから、f (xk ) ρk

k

Z β
ρk =

Lk (x)dxαとすれば良い。これで論理が完結した。
纏めると、区間 [α, β] と荷重関数 w(x) が与えられ、積分Rβdx w(x) f (x) dx を計算したければ、次の手順に従えば、関数を 2N − 1 次多項式近似した時の積分結果が N 回の関α

数値計算で得られる。
１) この積分区間と荷重関数に対する直交多項式を作り、pN (x) = 0 の N 個の独立根を計算し、これを xi , (i = 1, 2, · · · N ) とする。
２) 次の積分を行う

Rβ

ρi =

α

w(x) Li (x) dx
P

３) 関数値との積和を計算する。

i

ρi f (xi ) = I

これで、積分が計算出来た。上の説明で分かると思うが、直交多項式は規格化されているという事は何処にも使用していない。
特に α = −1, β = 1, w(x) = 1 の場合には、pN は Legendre の多項式である。積分の定義式や荷重関数に関して色々な場合を調べてみよ。

(xi , ρi ) の値は、数表を利用しても良いし、自分で計算しても良い。ρi はここでは、原理を示す事に留めたが、実際的な計算には直交多項式の教科書を調べ、直交多項式の基本漸化式や Christoffel-Darboux の定理を勉強すると良いだろう。

pN (x) = 0 をニュートン法で解くのに、p0N (x) を含む漸化式を利用出来る場合がある。更に、ρi の計算にもこの微係数が利用できる場合がある。
ひょっとすると、Pl m は l ≤ m に対してのみ定義されていると了解しているかも知れない。上では、l = 0, m = 1 での Pl m の値を使用した事に注意しておこう。０という初期値は簡単で使い易い。
このアイデアを用いて実用的な漸化式を作ってみよう。
独立な二つの角運動量 ~1 , ~2 を合成して ~3 を作る場合に Clebsch-Gordan 係数が登場する。次の式を定義式とする。

|j3 m3 >=X
m1 +m2 =m3|j1 m1 > |j2 m2 >< j1 m1 j2 m2 | j3 m3 >8

m1 と m2 の和がいつも m3 になるような条件付でとられる。ここで、次の演算子に関する恒等式を上の定義式に作用させる。

~23 = ~21 +~22 + 2~1 ·~2
角運動量の演算子と固有関数の間には、以下の関係式が成立する。

J± |j m >=p(j ∓ m)(j ± m + 1) |j m ± 1 >J 2 |j m >= j(j + 1) |j m >これらの式を用い、更に CG 係数の定義式を用いると以下の漸化式を得る。但し、簡単の為に、< j1 m1 j2 m2 | j3 m3 >= C(m1 , m2 ) と書く。

{j3 (j3 + 1)−j1 (j1 + 1) − j2 (j2 + 1) − 2m1 m2 }C(m1 , m2 )p
= (j1 + m1 )(j1 − m1 + 1)(j2 − m2 )(j2 + m2 + 1) C(m1 − 1, m2 + 1)
p
+ (j1 − m1 )(j1 + m1 + 1)(j2 + m2 )(j2 − m2 + 1) C(m1 + 1, m2 − 1)
これが求める３項漸化式であり、図の様に m1 , m2 を横、縦軸とすると、

m2
j

斜め４５度の斜線上に相隣る３個の格

2

j

1

m1

子点での C(m1 , m2 ) の値の関係を与える。初期値は、上で注意したように

C(m1 (= m3 − m2 ) m2 (= j2 + 1)) = 0 とする。従って、最初は実質的に２項漸化式となる。

即ち C(m3 − j2 + 1, j2 − 1) をある定数 C(m3 − j2 , j2 ) で与える式となる。このある定
数を C ≡ C(m3 − j2 , j2 ) とおく。ある定数をそのまま残したまま、３項漸化式を m1 が増える方向に使用すると、全ての項がある定数 C のみを含む形で確定する。ある定数 C は以下の規格化条件から決定できる。

X

C(m1 , m2 )2 = 1

m1 +m2 =m3ある定数の符号は、今の場合正にとっておく。
この３項漸化式を用いた Clebsch-Gordan 係数の計算法は、量子力学の教科書に書いても良いほどの出来栄えであり、初期値、漸化式、規格化条件と符号の規約で構成されている。
ついでにもう一言。ここで与えた漸化式は磁気量子数に関するものであったが、Clebsch-

Gordan 係数には j に関する漸化式もある。例えば M.E. Rose 著、森田訳の角運動量の基9礎理論 (みすず書房) の付録を参照の事。ここで与えた手法が j に関する漸化式にもほとん
どそのまま使用できる。これらの Clebsch-Gordan 係数の計算法は、実用上もかなり使える方法である。極端な場合には、傾向と対策を施す必要があるが。。。
ついでに 6j 係数の３項漸化式の導き方の指針を記録しておこう。ランダウ・リフシッツ、(好村・井上訳) の量子力学２ pp 483-484 を参考する。先ず、Elliot-Biedenharn の恒等式を利用する。

(

j1 j2 j3
j4 j5 j6)(=

j6 j1 j5
j7 j8 j9X)(

φ(j)(2j + 1)j

但し位相は φ(j) = (−1)j+

P

i ji

j4 j2 j6
j9 j8 j) (j2 j1 j3
j7 j j 9)(

j4 j 3 j 5
j7 j 8 j)である。ここで、j7 = 1, j8 = j5 , j9 = j1 を代入すると、jに付いての和は j = j3 , j3 ± 1 の３項が登場する。j7 = 1 を含む 6j には具体的な表式 (参考書第１０表) を代入する。後は、見やすい様に整頓するだけだ。
次はもう少し複雑な例として、クーロンポテンシャルの下での、束縛状態の波動関数を取り上げよう。この問題に対するシュレーディンガー方程式を極座標で書き、解を球面調和関数で展開する。動径部分の波動関数 R` (r) は、以下の微分方程式で与えられる。

1 d
r2 dr

µ

r

2 d R`dr¶µ2 η κ ` (` + 1)
+
− κ +
r
r2
2

¶

R` = 0

ここで、以下の記号が簡略化の為に導入された。

κ2 = −

2mE
,
~2

η=

Z z e2 m
4 π ² 0 ~2 κ

この方程式には二つの独立な解が存在し、大きな r に対して一方は r と共に指数関数的に増大し、他方は減衰する。束縛状態では、後者の解を採用する。独立変数を、無次元化する
ために、ρ = κ r と変更しておく。各種の ` に対して、大きな r に対してこの波動関数を計算したいという目的を設定してみよう。例えば反陽子原子やリュードベリ原子と呼ばれる非常に大きなエネルギー量子数を持つ原子を想定してみると良いだろう。
大きな r を想定し、1/r で解を展開してみる。但し、無次元量 ρ = k r を独立変数とする。
結果だけを書くと、R` (r) = ρ` e−ρ W1 (` + 1 + η, 2` + 1, ρ)=∞XΓ(` + 1 + η + n) 1
Γ(2` + 2) ` −ρρ e (2ρ)−`−1−ηΓ(` + 1 + η)Γ(` + 1 − η − n) (2ρ)nn=0∞X(2` + 1)!!
e−ρΓ(` + 1 + η + n) 1
=
`+1+η
1+η
Γ(` + 1 + η) 2ρΓ(` + 1 − η − n) (2ρ)nn=010ここで、W1 (α, β, z) は Kummer の合流型の超幾何関数である。教科書に良くでて来る、
水素原子型の束縛状態では、上の式の分母に登場したΓ関数があるところで発散するという条件から、クーロンパラメータ η が整数という条件をつける。
合流型の超幾何関数に慣れていない学生は、抵抗を感ずるかもしれないが、通常の級数展開による解法で、最後の式にたどり着くのは、困難でないだろう。
数学公式集には W1 に対して、パラメータの値が１だけ異なる場合の３項漸化式が登録されているが、今の場合、α = `+1+η, β = 2`+1 であるから、W1 (α, β, ρ), W1 (α±1, β ±2, ρ)というβに対しては一つ飛ばしの漸化式を作らねばならない。
(1)

(2)

β

(3)

β

f

β

g

e
d

α

α

α

b

c

a

問題を右側の図で言うと、a, d, g に関する漸化式を作りたいとする。即ち、連続する各種の軌道角運動量に対する負エネルギークーロン波動関数を ρ を固定して計算したい。参考書には、W (α, β, z) に関して、上の図の左の３個の漸化式が記載されている。これらの漸化式の特徴は W (α, β, z) に対し、(１) では β を１増し、α は固定又は１増した時の関係
式、(２) では α は１増し、β を固定した時と、１だけ増した時、(３) は、α は固定し、β のみ１だけ増減させた時の関係である。
次の手順で６個の未知数 a, b, · · · , g に関する５個の関係式を作る。(１) より、a, b, c、(２) より、b, c, d、(３) より c, d, e、(３) より d, e, f 、(１) より e, f, g 。b, c, e, f を消去し、a, d, g に関する漸化式を得る。
沢山の式を書くのは面倒だから、注意点だけを述べよう。上の手法で、W (α, β, z) に関する漸化式から、R` に関する３点漸化式を作ってみると、η が自然数 (即ち、束縛状態に対応する) の時には、分母が特定の ` に対して０となり、途中で漸化式が途切れてしまう。
この困難を回避するには、R` だけに注目するのでなく、その微分 R`0 (r) も同時に視野に入れておく事を勧める。微分方程式を数値的に解く時は、境界値問題を解いているはずであり、境界での関数値とその微分を同時に考慮するのが望ましいという、一般論に立って考え直すと良い。
実用的には、波動関数の数値や時にはその微分も必要となるから、これらの事を勘案した漸化式を作っておくのが望ましい。今の、負エネルギーのクーロン波動関数がその例になっている。
数学の教科書を見ると、α, β, z は複素数でも Kummer の合流型超幾何関数は定義されている。例えば実数に対するαが必要な時に、α ± i という様な関数を想定し、虚軸に沿っ11た漸化式を作ると、収束が早い場合がある。
折角大学で複素関数に付いての講義を聴いたのだから、こちらの方向へも目を向けてみよう。新天地が開ける場合がある。

もう一つ微分方程式の解法に関する、知っておいて損をしない知識を書いておこう。２階の微分方程式 F (x, y(x), y 0 (x), y 00 (x)) = 0 が与えられていて、点 x0 での値 y(x0 ), y 0 (x0 )が既知だとする。２階微分方程式 F = 0 を y 00 (x0 ) に付いて解くと y 00 (x0 ) も既知となる。
一般に、F = 0 を n 回微分すると、y の高階微分に関する漸化式が作られるが、上で得られた y, y 0 , y 00 を利用すると、任意の高階微分が芋蔓式に計算できる。これらの知識とTayler 展開を利用すると、x0 の近くの x の値が計算できる。
この様にして、微分方程式を数値的に解く事が可能である。
以下の様に纏められる。

0
1
2
3
4
5
6

微分方程式 F (x, y(x), y 0 (x), y 00 (x)) = 0 を与える。

y 00 (x) に付いて解く。
微分方程式を n 回微分し、y (n) (x) に関する漸化式をつくる。
境界条件より、初期値 x, y(x), y 0 (x), y 00 (x) を決定する。
任意の高階微分を計算する。
テイラー展開を用いる、y(x + h), y 0 (x + h) を計算する。

x + h を新しい x だとして [4] に戻る。
これで、微分方程式を高精度に計算する、手順の一般的な指導原理が出来た。
例として、ニュートンの万有引力下での運動方程式を取り上げよう。先ず、以下の様に変形する。r の n 階微分を r(n) 、C を定数として、

r3 r(2) = C r
この式を n 回微分する。

r3 r(n+2) + 3 n r2 r(n+1) + 3 n(n − 1) r r(n) + n(n − 1)(n − 2)r(n−1) = Cr(n)これで微分に関する４項漸化式が出来上がった。初期条件から r = r(0) , r(1) がきまると、運動方程式により r(2) が計算出来る。更なる高次微分は、この漸化式を利用すれば良い。
次の点を決めるには、許容誤差 ² と歩幅 h を与え、hn |r(n) |/n! < ² となる次数 n まで、高次微分を計算する。通俗本に登場する ４段４次の Runge-Kutta 法よりもはるかに有用な計算手法である事が理解出来よう。即ちここで与えた手法では、精度を上げるのに歩幅を減らす自由度と、次数を上げる自由度があるのに対して、Runge-Kutta 法では歩幅を減らす
自由度しかない。
12

１階微分を含まない、２階の微分方程式
少し異なる例として、微分演算子と微分方程式を取り上げよう。次の様な２階微分方程式を取り上げる。

d2 y(x)= f (x) y(x)dx2

この方程式は、１階微分を含まない線形の微分方程式と言う意味で、特殊な例とも考えられるが、ニュートンの運動方程式や、極座標表示のシュレーディンガー方程式にもこの形に変形できるものがある事は良く知られている。
この形の微分方程式の一般的な数値解法としては、 Henrici 著の ‘Discrete variable

Methods in Ordinary Differential Equations’ を参照すると良い。ここでは別の観点から取り掛かろう。微分演算子 (D) と中心差分演算子 (δ) 対にし以下の記号を定義する。
D f (x) ≡df (x), δf (x) ≡ f (x + h/2) − f (x − h/2)dx

ここで h は絶対値が小さなある定数である。中心差分の定義式の右辺をテイラー展開する。

δf (x) =X ½ 1 µ hD ¶n
n!

n

2

1
−
n!

µ

−hD
2

¶n ¾

f (x)

指数関数の定義を用いると、右辺は形式的に微分演算子を指数関数の肩に乗せられるから、

δf (x) = {eh D/2 − e−h D/2 }f (x) = 2 sinh (h D/2)f (x)この関係式は、任意の連続微分可能な f (x) に対して成立するはずである。f (x) をはずしておいて、逆に解くと、

2
1
D = sinh−1 δ/2 =
h
h

¾
½
3 δ5
δ3
+
···
δ−
24 640

最右辺では、逆正接関数を級数展開している。これで、微分演算子が、中心差分演算子で表現出来た。両辺を２乗すると、

1
D = 2
h
2

½

δ6
δ4
+
···
δ −
12 90

¾

2

これで２階微分演算子を、中心差分演算子で置き換える準備ができた。２階微分演算子はパリティーを保存するので、右辺では偶数次項のみが登場している。
ここから先が、数値的な工夫の部分である。両辺に 1 + δ 2 /12 を掛けて、O(δ 4 ) の項を消して、その上で O(δ 6 ) を無視する近似を導入する。その結果、(1 + δ 2 /12)D2 = δ 2 /h2
この式の両辺に右から y(x) を掛け、先に与えた微分方程式を代入すると、左辺は、(1 + δ 2 /12)D2 y(x) = f (x)y(x)1+ {f (x + h)y(x + h) − 2f (x)y(x) + f (x − h)y(x − h)}1213一方、右辺は

δ2
1
y(x)
=
{y(x + h) − 2y(x) + y(x − h)}h2h2であるから、µ

h2 f (x + h)1−
12¶y(x + h) = 12y(x)¶µ¶µh2 f (x)h2 f (x − h)
−10 1 −
y(x) − 1 −y(x − h)1212

この式は、y(x) に関する３点漸化式になっている。従って、上に与えた微分方程式は, 点x, x − h での y の値を与えるとこの漸化式を用いて x + h での y の値を計算する事が出来る。この漸化式は、O(h5 ) 迄は正確であるから、手数が省ける割には高精度の近似式となっている。
別の書き方をすると、yn = y(n h), fn = f (n h), An = (1 − h2 fn /12)yn とおくと、

An+1 = 12yn − 10An − An−1
この式は、大きな数 yn と小さな数 h2 yn /12 を同列に扱っているという意味では、センスが悪い。次の様にすると、情報落ちを幾らか救う事が可能である。

Dn = An − An−1 とおき、次の手順に依る。
1
2
3

Dn+1 = Dn + h2 fn yn
An+1 = An + Dn+1
−1
yn+1 = {1 − h2 fn+1 /12} An+1

最後の式は、割算を含んでいるので、

yn+1 = {An+1 + 10 An + An−1 } /12と書き換える方が良い場合が多いだろう。この式を使うためには、Dn+1 = Dn + h2 fn (1 + h2 f n/12) An という更なる近似を導入する必要がある。
この手法を、連立方程式を解く場合に拡張する事が出来る。その場合、y はベクトルであり、f (x) は行列となるが、ほとんどの式はそのまま成立する。和、差、積は行列の和、差、積と置き換えれば良い。割算が、逆行列を掛ける事に対応する。逆行列の計算は非常に手間 (時間) を要するので、最後に提案した置き換えが意
味を持って来る。 ニュートンの運動方程式には、この手法は適さないかもしれない。微分方程式の数値解法には、予測子修正子法というのがある。この手法が適さない時には、修正子として使用できないか？と考えてみると良いだろう。
導き型をみれば分かるように、通常の４段４次のルンゲークッタ法よりも、計算が簡単であるのに、精度が良い点。
ここでの出発点として取り上げた微分方程式は、偶数階微分のみを含み、パリティーを保存する。従って 差分方程式に変換した時、奇数次階差を含まない様に変換出来た事が一つの重要な点である。

14

もう一つの重要点は、(1 + h2 δ 2 /12) を掛けて O(h4 ) の係数を０とした点である。これらにより、非常に単純ながら高精度な積分公式が作られた。

15