Bessel 関数の数値計算
桂田 祐史

2005 年 3 月 11 日 (小修正 06/11/21)
概要
円盤領域や円柱領域で熱方程式、波動方程式などの微分方程式を考えるとき必要になる Bessel
関数を数値計算する手段について説明する。

大規模工事中
ミスプリもたくさんあるはず、内容を信じないこと
計算法

1(構想を練っているところ…)

1.1

次数 0, 1, . . . , n に対する値を一斉に計算する方法与えられた n に対して、J0 (x), J1 (x), . . . , Jn (x) を一斉に、あるいは I0 (x), I1 (x), . . . , In (x) を一斉に計算するための、Miller の方法 (Miller algorithm) と呼ばれるアルゴリズムがある。これはBessel 関数の漸化式に基づくものである。
このアルゴリズムの誤差解析、サンプル・プログラムについては、最近 (2006 年) 出版された二宮編 (担当部分は吉田著) [1] を見よ。
n
n
∑
∑
同著には
dj Jn (x),dj In (x) を計算する Deuflhard の方法も載っている。
j=0

j=0

2

色々な環境での計算

2.1

とにかく C コンパイラーがあればNumerical Recipes に整数次の Bessel 関数, 変形 Bessel 関数のプログラムが載っているが、多分単精度だと思われる。
0, 1 次でよければ、大浦卓哉氏の「ベッセル関数 — 整数次」(http://www.kurims.kyoto-u.ac.
jp/~ooura/bessel-j.html) がある。これは C 以外に Fortran 版もある。

2.2

UNIX

UNIX の数学関数ライブラリィには double j0(double), double j1(double), double jn(int,double), double y0(double), double y1(double), double yn(int, double) がある。FreeBSD
や Linux, Cygwin には当然ソース・プログラムがついている。C, C++ や Fortran で気軽に使えるわけである。しかし…これだけである。変形 Bessel 関数などはない。
1

2.3

GSL

素晴らしいことにソースコードが公開されている GNU Scientific Library1 (GSL) の中にも Bessel関数がある。C, C++, Fortran から容易に利用できるわけで嬉しい。ライセンスが気にならない人
ならばこれで決まりか。




/*
*
* gcc -I/usr/local/include*/-L/usr/local/lib -lgsl -lgslcblas testbessel.c -o testbessel#include <stdio.h>#include <gsl/gsl_specfunc.h>/** gsl_sf_bessel_J0()* gsl_sf_bessel_J1()* gsl_sf_bessel_Jn()* gsl_sf_bessel_I0()* gsl_sf_bessel_I1()* gsl_sf_bessel_In()*/
int main(){
int i,n,nu,MaxNu = 3;
double x,dx,xmin,xmax;
n = 100;
xmin = 0.0; xmax = 10.0;
dx = (xmax - xmin) / n;for (i = 0; i <= n; i++) {x = xmin + i * dx;
printf("%f ", x);for (nu = 0; nu <= MaxNu; nu++)
printf("%f ", gsl_sf_bessel_In(nu, x));printf("\n");}return 0;
}


2.4



Java

検索してみたら色々なものが出回っているようである。そのうちColt Project (http://dsd.lbl.gov/~hoschek/colt/)
を試してみたところ、快適だったので紹介する。
使うための準備は、ダウンロードしたファイルを展開して出て来る colt.jar を適当なところにおいて、そのパス名を環境変数 CLASSPATH に含めればよい。例えば私は Cygwin 上でCygwin 上の tcsh でsetenv CLASSPATH ’c:\j2sdk1.4.2_06\lib\colt.jar;.;’としている。
1

検索すればすぐに見つかる。

2

cern.jet.math というパッケージの中に Bessel というクラスがあるので、例えばこんな感じで
利用する。


TestBessel.java// TestBessel.java
import cern.jet.math.*;public class TestBessel {public static void main(String []args) {System.out.println(""+Bessel.j0(1.0));
}
}




残念ながら変形ベッセル関数は 0, 1 次のものしかないようである。
なお、Colt ライブラリィを全部インストールするのは大げさだと思う場合は、

jar でクラス・ファイルを抜き出すjar xvf colt.jar cern/jet/math/Bessel.class
jar xvf colt.jar cern/jet/math/Constants.classとして cern/jet/math/{Bessel,Constants}.class を抜き出しておけば、CLASSPATH を設定する必要はない2 。

2.5

Mathematica

Mathematica では BesselJ[n,x] で、n 次の Bessel 関数 Jn が計算できる。

グラフを描く

g1=Plot[BesselJ[0,x],{x,0,10}]g2=Plot[Evaluate[Table[BesselJ[n, x], {n, 0, 10}]],{x,0,10}]のようにしてグラフが描ける 。
3

1
0.8
0.6
0.4
0.2
2

4

6

8

10

-0.2
-0.4

2.5.1

零点を求めるB.7 も参照せよ。
以前は、次のようにしていた (これ自体は不要になったが、参考になるところが残っていると思う
ので、削除せずに残す)。
2
3

UNIX ユーザーには説明の要もないでしょう。jar は tar コマンドと使い方が似ていて、気軽に使える。
Mathematica でグラフを保存するには、描画した直後に Export["Bessel.eps", g1] のようにすれば OK.

3


これは不要になった



以下のようにして零点を必要な精度だけ計算することもできる。


FindRoot[BesselJ[0,x]==0,{x,2},WorkingPrecision->50]
Table[FindRoot[BesselJ[0,x]==0,{x,m*Pi+2},WorkingPrecision->50],{m,0,20}]

たくさんの桁数が必要なときは、AccuracyGoal->桁数 のような指定が必要かもしれない (通常は “Automatic”)。また n が大きくなると、MaxIterations->200 のような指定が必要かもしれない (通常は “100”)。
Version 6 以降の Mathematica では、Bessel 関数の零点を計算するために BesselJZero[n,k] (Jnの k 番目の零点), BesselJZero[n,k,x0 ] (Jn の x0 より大きい k 番目の零点), BesselYZero[n,k](Yn の k 番目の零点), BesselYZero[n,k,x0 ] (Yn の x0 より大きい k 番目の零点) が用意されている。
数値計算するには、bjz[n_,m_]:=bjz[n,m]=N[BesselJZero[n,m],30]

のようにして、bjz[n][m] を使うと良いだろう (一度計算したら覚えておいて、それを再利用する)。
Version 5.2 まで、NumericalMath‘BesselZeros‘ というパッケージがあった (http://documents.
wolfram.com/mathematica/Add-onsLinks/StandardPackages/NumericalMath/BesselZeros.html)。
この中には、例えば BesselJZeros[], BesselJPrimeZeros[] が入っていて、前者は今でこそ BesselJZero[]があるので不要になったが (N[BesselJZero[1,2],30] のようなことが出来るから)、後者は今でも必要になる。Jn′ (z) (n = 0, 1, . . . , 10) の正の零点を小さい方から 20 個集めた表が欲しければ、

<<NumericalMath‘BesselZeros‘For[n=0,n<=10,n++,bjpz[n]=BesselJPrimeZeros[n,20,WorkingPrecision->30]]

のようにする。bjpz[n][[m]] は Jn′ の m 番目の零点の 30 桁の数値となる。

2.6

Maple

Maple では、BesselJ(n,x) で、n 次の Bessel 関数 Jn が計算できる。

グラフを描くplot(BesselJ(0,x),x=0..10)plot([seq(BesselJ(n,x),n=0..5)],x=0..10)plot({seq(BesselJ(n,x),n=0..5)},x=0..10)

あるいは





でグラフを描いたり 。
44Mathematica でグラフをセーブするには、描画した直後に Export["Bessel.eps", %] とすれば OK.

4

Jn の正の零点の小さい方から m 番目 νn,m を表す BesselJZeros(n,m) があり、必要な精度だ
け計算することもできる。
evalf[50](BesselZeros(0,1))
evalf[50](BesselZeros(0,1..20))



(Maple は便利だ…)

2.7

gnuplot

これはプログラミング言語とは言えないけれど、グラフ描きで有名な gnuplot では、besj0(),
besj1() が用意されている。従って単に plot besj0(x),besj1(x) とするだけで J0 , J1 のグラフが描ける。
残念ながら、一般の自然数 n に対して Jn (x) を計算する関数は用意されていない。

2.8

精度保証付き数値計算

整数次 Bessel 関数のテイラー展開は交代級数になるので、効率を気にしなければ素朴なアルゴリ
ズムで精度保証付き数値計算ができる。参考文献をあげておこう…

3



余談
• 庵原謙治『古典的 Bessel 函数入門』http://nalab.mind.kobe-u.ac.jp/HOME/iohara/doc/
Bessel.pdf は面白い。
• 大型計算機がまだ強かった時代にコンピューターとつきあい始めた筆者は、NUMPAC のようなサブルーチン・ライブラリィに一目も二目も置いているところがある。そういうわけで
NetNUMPAC (http://netnumpac.fuis.fukui-u.ac.jp/numpac/) には興味があるのだけ
れど…
5

4

2017/1

Netlib に AMOS5 というライブラリィがあることを知った。

README======== readme for AMOS =======A Portable Package for Bessel Functions of a Complex Argument
and Nonnegative Order
This algorithm is a package of subroutines for computing Bessel
functions and Airy functions. The routines are updated
versions of those routines found in TOMS algorithm 644.


Bessel 関数早分かり

A

やっぱり学生は迷子になりがちなので、少し道案内をするつもりでこの節を書く。
なお、Bessel 関数とある程度長くお付き合いすることになる場合、何はともあれ Watson [2] とボウマン [3] を入手することを勧める (中途半端な本を何冊も買うよりも良い)。[2] は古い洋書だが、Amazon で容易かつ安価に入手できる。
勉強する場合に必要となる予備知識は複素関数論である。

A.1

Bessel 関数とは

A.1.1

ご対面

Bessel 関数とは (親族まで含めると色々あるが)、以下の式で定義される関数である。ただし ν は
複素数の定数である。
∞( z )ν ∑(1)Jν (z) :=(2)Yν (z) = lim

2( z )2k
(−1)k
,
k!Γ(ν
+
k
+
1)
2
k=0

cos µπJµ (x) − J−µ (x).
µ→ν
sin µπJν (z) を ν 次の第 1 種ベッセル関数, Yν (z) を ν 次の第 2 種ベッセル関数 (あるいはノイマン関数)と呼ぶ。
老婆心ながら 複素関数論の多価解析関数に慣れている人には、上の説明で十分のはずだが、そうでない人は A.1.3 を見よ。また Γ は (言うまでもなく) ガンマ関数であるが、これについて付録 B.5を用意した。もっとも m が 0 以上の整数の場合は、Γ(m + 1) = m! ということだけ知っておけば取りあえず意味が取れるであろう。

5

http://www.netlib.org/amos/

6A.1.2

今学ぶ意義(急いでいるときは、ここをすっ飛ばしてもよいでしょう。)
Bessel 関数は特殊関数と呼ばれる一群の関数のうちの一種で、例えば円盤領域や円柱領域における微分方程式の問題で活躍する (だから「円柱関数」と呼ばれることもある)。歴史的には最初これらの領域における波動方程式や熱方程式の問題を解くために出現したが、名前のもとになったベッセル (Friedrich Wilhelm Bessel, 1784–1846, Minden に生まれ、Königsberg にて没する, 天文学に顕著な業績がある) は、Kepler の方程式を解くために Bessel 関数の基礎固めをすることになった。
余談 A.1 久賀 [4] には、二体問題の解も周期関数なのだから Fourier 級数展開出来るはずだ、と言って授業で解き出して (勇気がありますね)、うまく行かなかったという話が書いてある。クーガー先生、そういうのはベッセル関数の本 (例えばボウマン [3]) か、数理物理の本に載っています。
円盤にしろ、円柱領域にしろ、特殊な領域にすぎないことは否めない。特殊関数という「方法論」
は、どれだけ幅広い現象を扱えるか (計算により予知できるか) という汎用性の面からはコンピュー
ター・シミュレーションに及ばない。
しかし特殊関数には、コンピューター・シミュレーションにはない利点もある。この方法で得られる解は解析的6 な関数なので、問題がパラメーターを含む場合はパラメーター依存性も目に見える
し、微分などの操作もでき (例えば電位を求めた場合は容易に電場が求まる)、数値計算では分かりにくい種々の情報が一時に得られる。特に解が固有値、固有関数の情報 (固有値の分布や、固有関
数の節がどうなっているか等) も含んでいるのはとてもありがたい場合が多い。たとえ特殊な場合(領域) に限られると言っても、得られる知見の大きさは無視できない。
とかく縁遠くなりがちな特殊関数の世界への入門として、Bessel 関数を学ぶ意義は大きい、と筆者が考えるようになった次第である。

A.1.3

ご対面の後に (すぐ分かること)

まず Jν (z) の定義式からすぐに分かることを調べてみよう。ここに書いてあることを理解するには複素関数論の (ごくごく初歩的な) 知識があればよい7 。
∑(1) のの部分はいわゆる z のベキ級数であり、その収束半径は (後で示すように) ∞ であるか
ら、複素平面 C 全体で正則な関数 (整関数と呼ぶのであった) を与える。例えば導関数を求めたい
場合は、単に項別微分すればよい。
ところが、一見簡単そうに見える (z/2)ν の部分が要注意である。ν が 0 以上の整数でないとき、
この式は(z/2)ν = exp (ν log(z/2))のように解釈すべきものであり、これは一般には整関数にはならないことを覚えているであろうか？
√
複素関数論の授業では、w = z や、w = log z という関数について学んだはずである。log として
いわゆる主値を取ると、C から負軸を除いた領域8Ω := C \ {x ∈ R; x ≤ 0}6
「解析的」という言葉は、普通は、(1) 式で書かれた、という意味と、(2) 局所的にベキ級数展開できる、という二つの意味で使われるが、今の場合、どちらの意味でも解析的である。
7
どうも複素関数論は「何も覚えていません」と言って簡単にギブアップする学生が多い (とても困る)。知識というものは使うことで初めて定着するものである。復習の良い機会だと考えて、面倒がらずに教科書やノートを開くべきである。
8(複素関数論を理解している人には余計なおせっかいであるが) 0 は関数 log の分岐点 (定義を書くのは省略する) であり、単に一価正則な関数を得るためには、この「要注意点」の回りを回ることができないように障害物を設ければよいということである。取り除くのが特に負軸である必然性はない。主値というときは、こうするのが「普通」というだ
けである。

7

で一価正則な関数を得ることができる。特に ν ∈ R のとき、任意の x > 0 に対して Jν (x) ∈ R となる。また ν ≥ 0 ならば Jν (0) にも普通に意味がつけられる:Jν (0) = 0 (ν > 0),

J0 (0) = 1.

この Bessel 関数をどのように導入し、その性質を調べて行くかについては、実は色々な流儀があるが、この文書では、(桂田研学生にとっての) 後の応用を考慮して、微分方程式を基礎として話を進めることにする。

身もふたもないまとめ

任意の ν ∈ C に対して、Jν (z) は (負軸を除いた領域) C \ {x; x ≤ 0} で正則である。特に ν が
0 以上の整数であるときは、Jν (z) は C 全体で正則である。


A.2



Bessel の微分方程式

ν を複素数の定数とするとき、

(
)
1 ′
ν2
y + y + 1− 2 y =0
x
x
′′

(3)

を (ν 次の) Bessel の微分方程式とよぶ。
Bessel 関数 Jν (x), Yν (x) はこの微分方程式の、解の基本系の一つ (一組) として得られる。つまり、この微分方程式の任意の解は
y = AJν (x) + BYν (x)の形に一意的に書けるし、逆にこの形の任意の関数はこの微分方程式の解になる。以下このことを
説明する。

A.2.1

「解の基本系」というものがあることについて大ざっぱな説明をすると、(3) は 2 階の線形微分方程式であり、解空間 (解の全体) は (通常の意味で) 2 次元の線型空間になる。ゆえに 1 組の関数からなる基底 (解の基本系) が存在する、ということである。
急いでいる人は、以上の説明でとりあえず満足して先に進むことをお奨めする。
寝た子を起こす話

常微分方程式の講義で、2 階の線型常微分方程式

y ′′ + p(x)y ′ + q(x)y = 0の解空間は 2 次元の線型空間であるという定理を学んだかもしれないが、そこでは係数 p(x), q(x)は普通の意味で素直な関数となっていたと思われる。我々が扱っている Bessel の微分方程式は 0 が係数の特異点になっているので、上の定理をそのまま使って一丁あがり、というわけにはいかない。
きちんと議論するには、まず特異点 0 を除いた区間 [a, b] (0 < a < b < ∞) に制限して考える。ここでは上の定理が使えて解空間は 2 次元の線型空間であり、関数の組 {φ1 , φ2 } が存在して、任意の解 y に対して、C1 , C2 が一意的に存在して(4)

y = C 1 φ1 + C 2 φ2が成り立つ。φ1 , φ2 は Ω := C \ {x; x ≤ 0} まで正則に拡張できる。これも証明を要することではあるが、それほど難しいことではない (線形だから解は定義域の端まで延長できる) し、我々の場合の Jν (z), Yν (z) はそのことを直接確認することもできるので、省略する。

8さて、こうして拡張した φ1 (z), φ2 (z) は微分方程式の解になる (一致の定理による)。微分方程式
は線形同次であるから、その線型結合 y もやはり微分方程式の解になる。
逆に Ω で微分方程式をみたす y があったとき、それを [a, b] に制限すると、すでに述べたことか
ら、適当な C1 , C2 が存在して (4) が成り立つ。ここで関数関係の延長原理 (これも一致の定理による) を用いると、Ω 全体で (4) が成り立つことが分かる。

A.2.2

Frobenius の方法で解を求める

おまじないとしか聞こえないかもしれないが、一般論を学んだ、あるいはこれから学ぶ機会のある人のために、「Bessel の微分方程式は、0 を確定特異点に持つ微分方程式であるので、Frobeniusの理論で扱うことができる」。
Frobenius の方法とは、
∞
∑
λ
y=x
bk x kk=0を微分方程式に代入して…両辺の係数比較により、よって...
以下、かなり長い議論になるので (書く暇が出来ると良いですね)、先に結果を述べると、Jν (x),J−ν (x) が解として求まる。


命題 A.1 任意の ν ∈ C に対して、Jν (x), J−ν (x) は Bessel の微分方程式 (3) の解である。



A.2.3



解の基本系

ν ∈ C \ Z であるとき、Jν (x), J−ν (x) は 1 次独立である (したがって (3) の解の基本系になる)。
これは両者の級数表示を見て、x → 0 のときの増大度を考えても容易に証明できるが、次の補題か
らも分かる。


補題 A.2 任意の ν ∈ C \ Z に対してW (Jν (x), J−ν (x)) = −2 sin νπ.
πx

ただし W (f, g) は関数 f , g のロンスキアンを表わす記号とする: W (f, g) :=f g.
f ′ g′



証明



(省略)

余談 A.2 この補題から、ν ∈ C \ Z のとき、Jν (x), J−ν (x) は Bessel の微分方程式の解の基本系になるのだが、我々に取って最も重要な応用では ν ∈ Z なので、このことをあえて命題にまとめておくことはしない。
ロンスキアンの等式から想像できることだが、ν = n ∈ Z のとき Jν (x), J−ν (x) は 1 次従属になる。これは容易に確かめられる関係

J−n (z) = (−1)n Jn (z)から明らかである。
さて、いよいよ第二種ベッセル関数 (ノイマン関数ともいう) Yν (x) を導入しよう。

9まず ν ∈ C \ Z の場合にcos νπJν (x) − J−ν (x)sin νπとおく。これは Jν (x), J−ν (x) の線型結合であるから、(3) の解である。また補題 A.2 より容易にYν (x) :=(5)

W (Jν (x), Yν (x)) =2πx(ν ∈ C \ Z)が分かるので、Jν (x), Yν (x) は (3) の解の基本系である。
さらに n ∈ Z に対して広義一様収束極限

lim Yν (x)ν→nが存在する。そこで

Yn (x) := lim Yν (x)ν→nとおくと、これも (3) の解であるが、(5) より9W (Jn (x), Yn (x)) = lim W (Jν (x), Yν (x)) =
ν→n

2
̸= 0
πx

であるから、Jn (x), Yn (x) は (ν = n のときの) (3) の解の基本系である。
以上をまとめると次の定理が得られる。


定理 A.3 (Bessel の微分方程式の解の基本系) 任意の ν ∈ C に対して、Jν (x), Yν (x) は Besselの微分方程式 (3) の解の基本系である。




A.3

Bessel 関数の性質

A.4

Bessel 関数の零点

A.5

Fourier-Bessel 展開

命題 A.4 (Bessel 関数の直交関係式) ν ≥ −1/2 とするとき、∫ 1rJν (µν,n r)Jν (µν,m r)dr =

90Jν+1 (µν,n )2δnm2(n, m ∈ N)


正則関数列が広義一様収束するとき、各階の導関数も広義一様収束することに注意せよ。

10



∫ c
√
xf (x) dx が絶対収束するとき、定理 A.5 (Fourier-Bessel 展開) f : (0, c) → C が連続で、

0

f が有界変分であるような任意の閉区間 I ⊂ (0, c) の任意の x, 任意の n ∈ N ∪ {0} に対して、
f (x) =∞∑

Aj Jn (λj x)j=1が成り立つ。ただし2
Aj := 2c Jn+1 (λj c)2∫ c
xf (x)Jn (λj x) dx (j ∈ N)
0

で、λj は Jn の正の零点を小さい順に並べたものである:

{x > 0; Jn (x) = 0} = {λj }j∈N ,0 < λ1 < λ2 < · · · < λj < λj+1 < · · · .

√
注意 A.6 (不連続関数の展開) f が連続でなくて、単に xf (x) が (0, c) 上 Lebesgue 下積分とするだけで、結果の等式を∞
f (x + 0) + f (x − 0) ∑=Aj Jn (λj x)2
j=1で置き換えたものが成立する。 — 普通の Fourier 級数と同じだ。


定理 A.7 円盤 {(x, y); x2 + y 2 < 1} の完全系
∞
∞ ∑
∞
∑
1∑
f (r, θ) =A0j J0 (µ0j r) +Jn (µnj r)(Anj cos nθ + Bnj sin nθ)
2 j=1
n=1 j=12
Anj =πJn+1 (µnj )22
Bnj =πJn+1 (µnj )2f (r, θ)Jn (µnj r) cos nθ dθ dr,
0

0

0

0

)

∫ 1 (∫ 2π



証明

)

∫ 1 (∫ 2π

f (r, θ)Jn (µnj r) sin nθ dθ dr.
r ∈ (0, 1) を固定して、θ 7→ f (r, θ) を考えると、Fourier 級数展開∞a0 (r) ∑f (r, θ) =+(an (r) cos nθ + bn (r) sin nθ),2
n=1

1
an (r) :=π∫ 2π
f (r, θ) cos nθ dθ,01
bn (r) :=π∫ 2π
f (r, θ) sin nθ dθ0が得られる。an (r), bn (r) は (0, 1) 上の関数だから、Fourier-Bessel 展開ができる。

an (r) =∞∑Anj Jn (µnj r),

j=1

bn (r) =∞∑Bnj Jn (µnj r),

j=1

2
Anj :=
Jn+1 (µnj )2
2
Bnj :=
Jn+1 (µnj )2

以上まとめると定理を得る。

11∫ 1
ran (r)Jn (µnj r) dr,0∫ 1
rbn (r)Jn (µnj r) dr.
0





系 A.8 円盤 {(x, y); x2 + y 2 < R2 } の完全系
∞
∞ ∑
∞
∑
1∑
f (r, θ) =A0j J0 (µ0j r) +Jn (µnj r)(Anj cos nθ + Bnj sin nθ)
2 j=1
n=1 j=1

2
Anj =2
πR Jn+1 (µnj R)2
2
Bnj =
2
πR Jn+1 (µnj R)2


A.6

)

∫ R (∫ 2π

f (r, θ)Jn (µnj r) cos nθ dθ dr,
0

0

0

0

)

∫ R (∫ 2π

f (r, θ)Jn (µnj r) sin nθ dθ dr.


応用 1: 円盤領域のラプラシアンの固有値問題Ω = {(x, y) ∈ R2 ; x2 + y 2 < 1} におけるラプラシアンの固有値問題(6)− △ u = λu (in Ω),(7)

u = 0 (on ∂Ω),(8)u ̸≡ 0を考える。すなわちこの 3 つの方程式をみたす関数 u と定数 λ を求める。λ を固有値、u を λ に属する固有関数とよぶ。
極座標変換
x = r cos θ, y = r sin θ, U (r, θ) = u(x, y)を導入すると(10)

1
1
Urr + Ur + 2 Uθθ = −λUrr
U (1, θ) = 0 (θ ∈ [0, 2π]),(11)

U ̸≡ 0.

(9)(0 < r < 1, θ ∈ [0, 2π]),

変数分離解を求める。すなわち

U (r, θ) = R(r)Θ(θ)の形をしているものを求める。(9) に代入して R(r)Θ(θ) で割るとR′′ 1 R′
1 Θ′′
+
+ 2
= −λ.
R
rR
r Θ
移項して

Θ′′
r2 R′′ + rR′ + λr2 R=− .
RΘ明らかにこの等式の値は定数である。それを µ とおくと、
Θ′′ = −µΘ.

Θ が周期 2π の関数であることを考えると、(∃n ∈ N ∪ {0}) (∃A, B ∈ C) s.t. µ = n2 ,Θ(θ) = A cos nθ + B sin nθ.

µ = n2 のとき、R(r) については
r2 R′′ + rR′ + (λr2 − n2 )R = 0.
12

すなわち

)
(
1 ′
n2
R + R + λ − 2 R = 0.
rr′′

ここで λ > 0 を仮定して

z=
とおくと、

λf

′′

(√

ゆえに

√

λr,

R(r) = f (z)

√

) (
(
√ )
n2
λ ′ (√ )λr +fλr + λ − 2 f
λr = 0.
r
r
)

(
)
n2
1 ′
f (z) + f (z) + 1 − 2 f (z) = 0.
zz′′

これは Bessel の微分方程式である。f (z) は z = 0 で有界であるから、∃C ∈ C s.t. f (z) = CJn (z).
すなわちR(r) = f (z) = f
ところで R(1) = 0 であるから、C ̸= 0 であるから

√

(√

)
√
λr = CJn ( λr).

√CJn ( λ) = 0.

λ は Jn の零点である。
∃m ∈ N s.t.

√λ = µnm .

これから λ = µ2nm .
注意 A.9 (λ > 0 を仮定したことについて) − △ u = λu の両辺に u をかけて積分すると∫∫∫
2∥∇u∥2 dx > 0 (非自明な固有関数は定数でないから)
λ |u| dx = − △ u u dx =ΩΩΩより λ > 0 が分かる。あるいは最大値原理を使って λ > 0 であることを証明することもできる。
一方、計算だけで λ > 0 を示すこともできる (工事中)。
以上から、λ := µnm は円盤領域におけるラプラシアンの固有値であり、unm (r, θ) := Jn (µnm r)(A′ cos nθ + B ′ sin nθ)は λ に属する固有関数であることが分かる。
実はこれ以外に固有値はなく、固有関数もここに現われるものだけで十分であることが分かる。
(工事中)
円盤領域における熱方程式
単位円盤 Ω における

ut = κ △ u (in Ω × (0, ∞)),u(x, t) = 0 (on ∂Ω × (0, ∞)),u(x, 0) = f (x) (x ∈ Ω)
13

の解は
∞
∞ ∑
∞
∑
)
(
1∑
u(r, θ, t) =A0m J0 (µ0m r) +exp −κµ2nm t Jn (µnm r)(Anm cos nθ + Bnm sin nθ),2 m=1
n=1 m=12
Anj =πJn+1 (µnj )22
Bnj =πJn+1 (µnj )2)∫ 1 (∫ 2π

f (r cos θ, r sin θ)Jn (µnj r) cos nθ dθ dr,
0

0

0

0

)

∫ 1 (∫ 2π

f (r cos θ, r sin θ)Jn (µnj r) sin nθ dθ dr.

問

半径 R の円盤の場合は？

問

Neumann 境界条件の場合は？

円盤領域における波動方程式1
utt = △ u (in Ω × (0, ∞)),c2u(x, t) = 0 (on ∂Ω × (0, ∞)),u(x, 0) = ϕ(x) (x ∈ Ω),
ut (x, 0) = ψ(x) (x ∈ Ω)
の解は
∞
∞ ∑
∞
∑
1∑
u(r, θ, t) =A0m J0 (µ0m r) cos(µ0m t) +Jn (µnm r) cos(µmn t) (Anm cos nθ + Bnm sin nθ)
2 m=1
n=1 m=1+∞ ∞∞sin(µ0m t) ∑ ∑sin(µnm t)
1∑+(Cnm cos nθ + Dnm sin nθ) .
C0m J0 (µ0m r)
Jn (µnm r)
2 m=1
µ0m
µ
nm
n=1 m=1

B

おもちゃ箱

B.1

Kepler の方程式

離心近点角 E, 平均近点角 M , 楕円の離心率 e ∈ [0, 1) についてE − e sin E = Mが成り立つが、これを Kepler の方程式と呼ぶ。この解は Bessel 関数を用いて

E=M+

∞
∑
2

n
n=1

Jn (ne) sin nMと解ける。
(M と e が与えられて E を求めるには、ということだろうけれど、数値的に求めるだけならば、直接方程式を Newton 法で解くのが早いだろうな。まあ、パラメーター依存性とか見るにはよいのだろうけれど。)

14



√



(書きかけ) 楕円軌道 (x/a)2 + (y/b)2 = 1 で、焦点 (c, 0) (c = a2 − b2 ) に「太陽」があるとする。近日点 (a, 0) が π = (1, 0) となるように縮尺を変えると、太陽の位置は S = (e, 0) となる(e := c/a)。単位円上を等角速度で動く仮想の天体を考える。ただし、周期と近日点を過ぎるタ
イミングは問題としている惑星と同じとする。ある時刻に平均近点角が M であるとは、その仮想の天体の座標が (cos M, sin M ) ということである。問題の惑星の位置を P = (x, y) とする。
√
単位円上の点 P′ を P′ = (x, 1 − x2 sign y) で定める。E = ∠πOP′ を離心近点角という。E から真近点角 ∠πSP = ν を求めるには
√
ν
1+e
E
tan =
tan
2
1−e
2
を解けば良い。


余談 B.1 Kepler の方程式は文献に現れた初めての超越方程式ということであり、これまでに何百という解法が考えられてきたそうである (平凡社百科事典の記述 (by 堀 源一郎) による)。

B.2

変形 Bessel 関数

(
)
1 ′
ν2
y + y − 1+ 2 y =0xxを変形 Bessel の微分方程式と呼ぶ。
x = ix′ と変数変換すると通常の Bessel の微分方程式が得られる。これから Jν (ix), Yν (ix) が解の基本系になるが、これは分かりづらいので、∞( x )2m( x )ν ∑1,Iν (x) :=2 m=0 m!Γ (ν + m + 1) 2′′

2
I−α (x) − Iα (x)
lim
α→ν
π
sin απとおくと、Iν (x), Kν (x) が解の基本系になる。このうち x = 0 で連続な解を与えるのは Iν (x) の方だけである。
なお、()
1
Iν (x) = exp − νπi Jν (ix)2という関係がある。
なお、
(
)
ν2
1 ′
2
′′
y + y − k + 2 y=0xxの一般解は (通常の Bessel の微分方程式と同様に)Kν (x) :=y = AIν (kx) + BKν (kx).

B.3

零点

n 次 Bessel 関数の正の零点を小さい順に並べたものを
0 < νn,1 < νn,2 < · · · < νn,m < νn,m+1 < · · ·とおくとき、次の15

(1) Jn と Jn+1 の零点は重なることなく交互に並んでいる。
νn,1 < νn+1,1 < νn,2 < νn+1,2 < · · · < νn,m < νn+1,m < νn,m+1 < νn+1,m+1 < · · ·(2) lim νn,m = ∞.
m→∞(3) lim (νn,m+1 − νn,m ) = π.
m→∞(4) lim (νn+1,1 − νn,1 ) = 1.
n→∞

B.4

台形公式による数値計算

母関数 (generating function)

[

(
)] ∑
1
1
exp z t −=Jn (z)tn .
2tn∈Z
留数定理から、原点を正の向きに一周する閉曲線 C に対して
[ (
)]
I
1
1
1
exp z t −Jn (z) =t−(n+1) dz.
2πi C
2
t
これから

1
Jn (z) =
2π

∫ 2π
e

i(z sin θ−nθ)

0

1
dθ =
2π

∫ 2π
cos(z sin θ − nθ) dθ.
0

最後の積分は、解析的周期関数の一周期に渡る積分なので、台形公式で高精度に計算できる。例
えば、J4 (5) の計算に、区間を 32 等分した台形公式で 10−16 オーダーの精度が出る (桂田 [5])。

B.5

ガンマ関数

(Euler の) ガンマ関数とは、Re z > 0 なる z ∈ C に対して次式で定義される。
∫ ∞Γ(z) =
e−x xz−1 dx (Re z > 0).
0

これは正則関数なので、解析接続によって拡張され、0 以下の整数全体が極全体となる有理型関数となる。
実変数の関数とみたときのグラフを次に掲げる。
20
15
10
5
-4

-2

2
-5
-10
-15

16

4

任意の z ∈ C に対して

Γ(z + 1) = zΓ(z)がなりたち、特に n ∈ N に対して

Γ(n) = (n − 1)!.

B.6

漸近挙動

0 < x ≪ 1 のとき
Jn (x) ∼

1 n
x
n
2 n!

(n ∈ N ∪ {0}),Y0 (x) ∼2
log x,π

Yn (x) ∼ −2n (n − 1)! −nxπ(n ∈ N).

x ≫ n のとき√
Jn (x) ∼

B.7

√

(
2
π)
cos x − (2n + 1),
πx
4

Yn (x) ∼(2π)
sin x − (2n + 1)
.
πx
4

零点の表


いいかげん Mathematica プログラム(* n 次第 1 種 Bessel 関数の x0 の近くの根 *)rootnear0[n_, x0_]:= FindRoot[BesselJ[n, x] == 0, {x, x0},
WorkingPrecision -> 100, AccuracyGoal -> 90](* n 次第 1 種 Bessel 関数の x0 の近くの根 *)rootnear[n_,x0_]:= x /. rootnear0[n,x0][[1]](* n 次第 1 種 Bessel 関数の m 番目の根 *)Jroot[0,1]=rootnear[0,2]Jroot[n_,1]:=Jroot[n,1]=rootnear[n,Jroot[n-1,1]+Pi/2]
Jroot[n_,m_]:=Jroot[n,m]=rootnear[n,Jroot[n,m-1]+Pi](既に述べたように、現在の Mathematica では、もっと楽に零点を求めることが出来る。)


m > 1 に対する νn,m は、Jn (x) = 0 の νn,m−1 + π の近くの解として求めている。それはlim (νn,m+1 − νn,m ) = πm→∞に基づく。根拠はあまり強くない (n が小さいときは、m が小さくても νn,m+1 − νn,m ∼ π は良い近似になっているが、n が大きくなるとずれが大きくなって来る)。
m = 1 の場合、つまり νn,1 をどう求めるか。

lim (νn,1 − νn−1,1 ) = 1n→∞という結果を使って、νn−1,1 + 1 の近くの零点として探すか、それとも

νn−1,1 < νn,1 < νn−1,2を使って、νn−1,1 と νn−1,1 の真中、大体 νn−1,1 + π/2 の近くの零点として探すか。

17

1.16
1.14
1.12
20

40

60

80

100

1.08
1.06
1.04

νn+1,1 − νn,1 はどれくらいか？
ListPlot[Table[Jroot[n,1]-Jroot[n-1,1],{n,100}]]
νn−1,1 + 1 くらいで良いのかな？後は FindRoot[] の仕様を調べないと。
B.7.1
1
2
3
4
5
6
7
8
9
10
B.7.2
1
2
3
4
5
6
7
8
9
10

J0 の零点2.40482 55576 95772 76862 16318 79326 45464 31242 44909 14596 71357
5.52007 81102 86310 64959 66041 12813 02742 52218 65478 78290 98537
8.65372 79129 11012 21695 41987 12660 94668 55657 95231 27535 56188
11.79153 44390 14281 61374 30449 11925 45892 20229 24699 69544 67032
14.93091 77084 87785 94776 25939 97388 68220 79158 50115 63302 81587
18.07106 39679 10922 54314 78829 75618 17656 02489 86747 00132 60864
21.21163 66298 79258 95907 83933 50526 30683 61818 08975 97639 98327
24.35247 15307 49302 73705 79447 63178 90718 45693 72675 14892 70224
27.49347 91320 40254 79587 72882 34607 41454 65295 68860 54962 20109
30.63460 64684 31975 11754 95789 26854 23273 72735 71629 17814 71907
J1 の零点3.83170 59702 07512 31561 44358 86308 16076 65645 45274 28780 19287
7.01558 66698 15618 75353 70499 81476 52474 32763 11502 91131 38960
10.17346 81350 62722 07718 57117 76775 84406 98195 12500 19168 55561
13.32369 19363 14223 03239 36841 26947 87675 12166 44731 35786 57854
16.47063 00508 77632 81255 24604 70989 55144 94381 26822 27312 57699
19.61585 85104 68242 02112 50658 84137 50985 02474 02661 88054 46473
22.76008 43805 92771 89805 30051 52182 25759 29053 70738 07322 68720
25.90367 20876 18382 62549 58554 45979 87428 79054 27031 36724 76413
29.04682 85349 16855 06664 78198 83531 96110 04141 71793 08387 56660
32.18967 99109 74403 62662 29841 04460 36921 90528 67711 01497 69203

18

B.7.3
1
2
3
4
5
6
7
8
9
10
B.7.4
n
0
1
2
3
4
5
6
7
8
9
10
B.7.5

J2 の零点5.13562 23018 40682 55630 14016 90137 76545 69737 72347 50055 09433
8.41724 41403 99864 85778 36136 76153 88467 37669 92080 49839 77125
11.61984 11721 49059 42709 41449 86755 33488 52627 09284 34531 37936
14.79595 17823 51260 74666 14713 20227 96527 25762 10042 80731 40849
17.95981 94949 87826 45511 51420 77253 73317 75646 82223 39859 93682
21.11699 70530 21845 59096 28170 82797 69448 96968 38559 34853 79677
24.27011 23135 73102 60958 31524 28183 40906 89636 70207 83166 67677
27.42057 35499 84557 33057 04881 91395 16056 47756 94492 23345 43440
30.56920 44955 16397 03660 40635 39306 27367 24964 14998 40280 13764
33.71651 95092 22699 92195 91984 25218 07642 50150 14168 14688 32921
第一零点

2.40482 55576 95772 76862 16318 79326 45464 31242 44909 14596 71357 06999
3.83170 59702 07512 31561 44358 86308 16076 65645 45274 28780 19287 62298
5.13562 23018 40682 55630 14016 90137 76545 69737 72347 50055 09433 58257
6.38016 18959 23983 50623 66146 41942 70330 53263 03691 90308 83271 07188
7.58834 24345 03804 38506 96300 07985 61741 73699 77901 31298 12110 15515
8.77148 38159 59954 01912 28671 33409 56056 29810 77014 89739 55086 45007
9.93610 95242 17684 89469 30891 26965 19193 15561 76774 40886 79173 76649
11.08637 00192 45083 84576 27644 35929 99914 02727 17259 47518 88322 10993
12.22509 22640 04655 17561 28047 69107 39895 12086 25498 16142 81632 21761
13.35430 04774 35331 06641 99248 83491 92217 62585 57055 64275 26769 50854
14.47550 06865 54541 23845 16376 55413 15197 63048 12629 97022 97678 75271
小さい方から 10 個の零点

µn,m
n m
2.40482555769577276862163187932645464312424490914596713570699909059676583 0 1
3.83170597020751231561443588630816076656454527428780192876229898991883930 1 1
5.13562230184068255630140169013776545697377234750055094335825725745998196 2 1
5.52007811028631064959660411281302742522186547878290985375755203814429082 0 26.38016189592398350623661464194270330532630369190308832710718828425111957 3 1
7.01558666981561875353704998147652474327631150291131389605537782698549601 1 2
7.58834243450380438506963000798561741736997790131298121101551578705267466 4 1
8.41724414039986485778361367615388467376699208049839771258553399526347304 2 28.65372791291101221695419871266094668556579523127535561889147658302259995 0 38.77148381595995401912286713340956056298107701489739550864500722086250787 5 1µ61 = 9.93 · · · , µ32 = 9.76 · · · , µ13 = 10.17 · · · であるから

19

解析学の復習

C

解析的係数を持つ線型常微分方程式は色々なテキストで取り扱われている。どれを参考にしても良いが、卒研で扱う問題との親和性で選ぶと10 、草野 [6] などが良いと思う。

C.1

ロンスキー行列式

R の区間 I で連続な p(x), q(x) が与えられたとき、
y ′′ + p(x)y ′ + q(x)y = 0
を考える。
二つの解 f , g があるとき、(W (x; f, g) := det)
f (x) g(x)f ′ (x) g ′ (x)とおき、これを f , g のロンスキー行列式と呼ぶ。
任意の x0 , x ∈ I に対して、)( ∫ xp(t) dt
W (x; f, g) = W (x0 ; f, g) exp −x0が成り立つ。



命題 C.1 次の (i), (ii), (iii) は互いに同値である。

(i) f , g は I で 1 次独立である。
(ii) ∀x ∈ I

W (x; f, g) ̸= 0.

(iii) ∃x ∈ I

W (x; f, g) ̸= 0.



C.2



解析的係数を持つ線型常微分方程式 (1) 正則の場合

線型常微分方程式a(x)y ′′ + b(x)y ′ + c(x)y = 0
において、すべての係数は開円盤 |x−x0 | < r0 で正則 (ベキ級数展開可能) とする。このとき a(x0 ) ̸= 0のとき x0 を正則点、a(x0 ) = 0 のとき x0 を特異点と呼ぶ。
x0 が正則点である場合、p(x) := b(x)/a(x), q(x) := c(x)/a(x) とおくことで、

y ′′ + p(x)y ′ + q(x)y = 0と変形できる。十分小さい r > 0 が存在して、p(x), q(x) は |x − x0 | < r で正則である。

10

境界値問題、固有値問題が主なテーマであるためだが、その条件を抜きにしてもサービス精神にあふれた良い本で
ある

20





定理 C.2 p(x), q(x) が |x − x0 | < r で正則とするとき、

y ′′ + p(x)y ′ + q(x)y = 0
について次の (1), (2) が成り立つ。

(1) 任意の解は |x − x0 | < r で正則である。
(2) 解空間は 2 次元の線型空間である。


C.3



解析的係数を持つ線型常微分方程式 (2) 確定特異点の場合|x − x0 | < r0 において解析的係数を持つ 2 階線型常微分方程式
a(x)y ′′ + b(x)y ′ + c(x)y = 0(12)において、x0 が特異点 (i.e. a(x0 ) = 0) である場合、p(x) := b(x)/a(x), q(x) := c(x)/a(x) とおくと、x0 は p(x), q(x) の極となる可能性がある。

y ′′ + p(x)y ′ + q(x)y = 0(13)の係数 p(x), q(x) が x0 をそれぞれ高々 1 位の極、高々 2 位の極となっているとき、x0 を確定特異点とよぶ。
(12) にもどると、x0 が a(x) の高々 2 位の零点であるならば、x0 は確定特異点ということになる。
例 C.3 (Euler の微分方程式) a, b を定数とするとき、

x2 y ′′ + axy ′ + by = 0
を Euler の微分方程式という。0 はこの方程式の確定特異点である。x = et とおくと定数係数常微分方程式に帰着して解くことができる。

C.4

一致の定理

(準備中)

参考文献
たけみつ

やすよ

ひろし

[1] 二宮市三, 吉田年雄, 長谷川武光, 秦野 世, 杉浦 洋 , 櫻井鉄也, 細田陽介：数値計算のわざ, 共立
出版 (2006), 「やす」は￥UTF752F で出力できるはず。
[2] Watson, G. N.: A Treatise on the Theory of Bessel Functions, 2nd edition, Cambridge University
Press (1944, 1995).
[3] フランク・ボウマン：ベッセル関数入門, 日新出版 (1963), 平野鉄太郎訳.
[4] 久賀道郎：ドクトル・クーガーの数学講座 1, 日本評論社 (1992).
[5] 桂 田 祐 史：数 値 積 分,
http://nalab.mind.meiji.ac.jp/~mk/lecture/ouyoukaiseki4/
notebook/numerical-integration.pdf (2000).
たかし[6] 草野 尚 ：境界値問題入門, 朝倉書店 (1971 (復刊 2004)).
21