ベッセル関数を含む超越方程式の解

質問

お忙しい所，回答の方お願いいたします。
私は現在，以下のようなベッセル関数を含む超越方程式

Xn*J1(Xn)=a*J0(Xn)(a:正定数,n=1,2,・・・)

の任意の正定数aを設定した時の解Xnを得たいと考えています。
おそらく数値的に解くのだと思いますが，様々な数値計算プログラミングに関する本を調べているのですが載っておらず，悩んでいます。

私の今の考えはまず

f(x)=x*J1(x)-a*J0(x)=0

の関数のプログラムを組んでグラフを作り，x軸と交わる位置から解の範囲を定めて，ニュートン法や2分法などを使用して解くというものです。
しかし，この方法に自信がありません。

もし「他に簡単な方法がある！」，「このやり方ではこういう点が難しく解けない！」などございましたら
回答の方，何卒よろしくお願いいたします。

回答No.1

adinat

高校数学でもそうですが、パラメータ付きの関数の零点を求める問題は、変数分離法というのが基本です。コンピュータでの数値計算でもそうですが、パラメータを入れてのグラフ描画などは難しいことが多いこと、また変数を分離させて残った関数が、わりに扱いやすい関数になる可能性があること、があるからです。とはいえ万能とは限らないので、もちろんdaitasuki様の提示されたような方法でも問題はありません。

以下、一例です。変数分離させて、xJ1(x)/J0(x)=aという方程式を考えます。左辺をf(x)=xJ1(x)/J0(x)とおきましょう。これは0次ベッセルの零点に一位の極（関数が定義できない点）を持つ、単調増大な関数です。フリーソフトでもこの関数を描画できるソフトがあります。またJ0とJ1の零点は交互に訪れます。たとえればsinとcosの関係ですね。またf(0)=0です。ここから単調に増大して、最初のJ0の零点=2.40482…で∞になります。ここでいったん-∞になって、また単調増加して、J1の最初の零点=3.83170…で0になります。さらに増大して、J0の2番目の零点=5.52007…で再び∞になります。また-∞になって単調増大して、以下これを繰り返します。ベッセル関数の零点は詳しく調べられていて、特殊関数の数表などいろいろなところに載っていますし、プログラムを組んで簡単に調べることもできます。このことを利用すれば、f(x)=aのn番目の解は、J1のn-1番目の零点からJ0のn番目の零点の間に存在することが分かります。これに留意してニュートン法なり2分法なりを用いて解の探索を行ってください。

この問題に限らず、ニュートン法を使うときに留意すべきことは、初期値のとり方を工夫しないとどの解に収束するかが分からない、ということです。大まかでもいいからグラフをイメージして、初期値をどこにとれば、狙った解に収束させることができるかを考えてからプログラミングすべきです。多少アナログ的な考えですが、そうすれば望んだ解が直ちに得られると思います。

以下J0とJ1の零点を参考までに10番目まで計算させましたので、参考にしてください。

J0
1st zero = 2.4048255576957727686216318793264546431242449091460
2nd zero = 5.5200781102863106495966041128130274252218654787829
3rd zero = 8.6537279129110122169541987126609466855657952312754
4th zero = 11.791534439014281613743044911925458922022924699695
5th zero = 14.930917708487785947762593997388682207915850115633
6th zero = 18.071063967910922543147882975618176560248986747001
7th zero = 21.211636629879258959078393350526306836181808975976
8th zero = 24.352471530749302737057944763178907184569372675149
9th zero = 27.493479132040254795877288234607414546529568860550
10th zero = 30.634606468431975117549578926854232737273571629178

J1
1st zero = 3.8317059702075123156144358863081607665645452742878
2nd zero = 7.0155866698156187535370499814765247432763115029010
3rd zero = 10.173468135062722077185711776775844069819512500192
4th zero = 13.323691936314223032393684126947876751216644731358
5th zero = 16.470630050877632812552460470989551449438126822273
6th zero = 19.615858510468242021125065884137509850247402661881
7th zero = 22.760084380592771898053005152182257592905370738073
8th zero = 25.903672087618382625495855445979874287905427031367
9th zero = 29.046828534916855066647819883531961100414171793084
10th zero = 32.189679910974403626622984104460369219052867711015

お礼

回答ありがとうございます。
また，わざわざJ0とJ1の零点の計算結果まで載せていただき，ありがとうございます。

非常に丁寧な説明をしていただき，関数f(x)=xJ1(x)/J0(x)のイメージがよく分かりました。

念のためプログラムを作ってグラフを書いてみましたが，イメージ通りでした。

これを元に，y=a(ためしにa=1.0)との交点から解の範囲を定めて二分法もしくはRegula-falsi法で数値的に解いた所，文献にある5番目(n=5)までの解(この解の解き方がこの文献に載ってなかったのです。)と良好な一致を示す数値解が得られました。

ずっと悩んでいた部分が解決され，すっきりしました。

改めて御礼を申し上げます。
ありがとうございました。