目次
合流型超幾何方程式の1つであるベッセルの微分方程式
d2ydx2+1xdydx+(1−ν2x2)=0
は、x=0にただ1つの確定特異点を持つ。一般解は、A,Bを定数として
y(x)=AJν(x)+BNν(x)
と表される。ここで、Jν(x):(\nu)次の第1種ベッセル関数 Nν(x):ν次の第2種ベッセル関数(ノイマン関数)である。xが実数でνが整数である場合のグラフは以下のようになる。
図1 第1種ベッセル関数のグラフ
図2 第2種ベッセル関数のグラフ
第2種ベッセル関数の値はx→0で発散する。第1種ベッセル関数については、0次のみx=0における値が非零である。
円筒座標系でラプラス方程式やヘルムホルツ方程式を解くときにベッセルの微分方程式( ν=0 )がよく出てくる。(たとえば円筒導体における電流密度の半径方向分布)
今回は、実数変数の整数次第1種ベッセルの微分方程式をアナログコンピュータで解いて、第1種ベッセル関数のグラフを描くことを目的とする。(以後νは整数nとして表記することにする。)
重要な公式をここで挙げておく[1]。
(a) 第1種ベッセル関数の級数表示
Jn(x)=(x2)n∞∑k=0(−1)k(x/2)2kk!Γ(n+k+1)
(b) 微分及び漸化式
J′0(x)=−J1(x)
J′n(x)=nxJn(x)−Jn+1(x)
2{Jn−1(x)+Jn+1(x)}=nxJn(x)
(5), (6)より
J′n(x)=12{Jn−1(x)−Jn+1(x)}
ベッセル微分方程式に限らず、最高微分階数項の係数に独立変数が掛かる方程式はアナログ計算機での陽関数法による実装が難しい。時刻0で式の分母が0になり、割算器が飽和するからである。各々の式の形に合わせて適切な実装法を模索せねばならない。
ただ、0次の第1種ベッセル関数のみ得られれば良いのなら、実はそこまで難しくない。
図1のグラフからわかる通り、0次のベッセル関数は初期値1が最大値であり、1次のベッセル関数も最大値は1以下であることから、電圧換算は大変楽である。演算方程式は原方程式そのものを使えばよい。陽関数法で演算する場合、
d2Jndx2+1xdJndx+(1−n2x2)Jn=0d2Jndx2=−1xdJndx+n2x2Jn−Jn
としてn=0を代入すると
d2J0dx2=−1xdJ0dx−J0
を得る。
しかしながら、このまま割算器で発生させた1/xをdJ0dxの項と掛け合わせると、1階微分項の係数がx=0で発散し、電圧信号としては飽和してしまうので、ひと工夫が必要になる。
式(3)の項別微分により、独立変数原点における微係数を求めると
dJndx=J′n(x)={0(n≠0)12(n=1)
であり、(9)式における信号dJ0dxはx=0で分母と分子どちらも0に近く、電圧信号としては有限値をとる。つまり、割算器でdJ0dx/xを直接生成することで、信号の発散を防ぐことができる。
(4)式から、1次のベッセル関数は0次ベッセル関数の微分値を反転させたものに等しいことがわかるので、式(9)の回路に反転器を追加するのみで求まる。ただし、n=1のみを生成する場合はこの方法は使えない。なぜなら、(8)式にn=1を代入した式
d2J1dx2=−1xdJ1dx+1x2J1−J1
において、独立変数原点における微係数 dJ1(0)dx は0でない値をとり、 dJ1(0)dx/xの項は飽和する可能性が高い。さらにJ1/x2の項もあらたに生じる。
一般のnの場合も同じような理由で困難が生じる。それ以上に深刻なのが、初期条件の問題である。
整数n次ベッセル関数のx=0における値(初期値)として、非零の値が出てくる最初の階数kは、k=nの場合である。それ以降は、偶数次であれば偶数階で、奇数次であれば奇数階で非零の初期値が出現する。
以下に、J0(0)=1としたときのJ(k)n(0)の値を示す。
表1 n次第1種ベッセル関数Jn(x)のx=0におけるk階微分係数
n\k | (0) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
0 | 1 | 0 | -1/2 | 0 | 3/8 | 0 | -5/16 | 0 | 35/128 | 0 | -63/256 |
1 | 0 | 1/2 | 0 | -3/8 | 0 | 5/16 | 0 | -35/128 | 0 | 63/256 | 0 |
2 | 0 | 0 | 1/4 | 0 | -1/4 | 0 | 15/64 | 0 | -7/32 | 0 | 105/512 |
3 | 0 | 0 | 0 | 1/8 | 0 | -5/32 | 0 | 21/128 | 0 | -21/128 | 0 |
4 | 0 | 0 | 0 | 0 | 1/16 | 0 | -3/32 | 0 | 7/64 | 0 | -15/128 |
5 | 0 | 0 | 0 | 0 | 0 | 1/32 | 0 | -7/128 | 0 | 9/128 | 0 |
6 | 0 | 0 | 0 | 0 | 0 | 0 | 1/64 | 0 | -1/32 | 0 | 45/1024 |
7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1/128 | 0 | -9/512 | 0 |
8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1/256 | 0 | -5/512 |
9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1/512 | 0 |
このことから、n=3以上の場合は Jn(x),dJnxdx(0)の2つの初期条件が両方0になることがわかるので、(11)式では有効な演算回路として構成できない。
(4)の漸化式 J′n−1(x)=(n−1)Jn−1(x)x−Jn(x) を用いて、(n-1)次のベッセル関数とその導関数からn次のベッセル関数を次々と求める手法も考えられるが、その場合割算器の被除数(分子)のx=0における値がn≥1で0となってしまい、そのうえn個もの除算器が必要になる。非線形演算要素の数がコスト面や精度面でのボトルネックになりやすいので、その点からもあまりお薦めできない。
n次ベッセル関数Jn(x)のアナログコンピュータによる安定した生成において必要な条件は、
となるだろう。
これらの条件をクリアしてn≥1次ベッセル関数を生成する手法としていくつかの演算回路構成法が考案されているので、ここでまとめてみる。
n次ベッセル関数が満たすベッセル微分方程式
x2d2Jndx2+xdJndx+(x2−n2)Jn=0
の両辺を独立変数xについて1回微分し、最高階について整理する。
−d3Jndx3=3xd2Jndx2+(1+1−n2x2)dJndx+2xJn
さらに、n=1を代入すれば
−d3J1dx3=3xd2J1dx2+dJ1dx+2xJ1
となる。上式は前述の条件を満たすので、これを演算方程式とする回路から J1(x)が得られる。
今度は(13)式を2階微分してn=2を代入し、最高階について整理すると
−d4J2dx4=5xd3J2dx3+d2J2dx2+4xdJ2dx+2x2J2
(16)も条件を満たすので、この式からJ2(x) が得られる。
一般のnの場合も同様にして
−dn+2Jndxn+2=2n+1xdn+1Jndxn+1+dnJndxn+2nxd(n−1)Jndxn−1+n(n−1)x2dn−2Jndxn−2
(16)式をもとに演算方程式をつくると
−aτn+2aJdn+2Jndτn+2=aτn+1aJ(2n+1)aTTdn+1Jndτn+1+aτnaJdnJndτn+2naτn−1aJaTTdn−1Jndτn−1+n(n−1)aτn−2aJaT2T2dn−2Jndτn−2−aτn+2dn+2Jndτn+2=aτn+12n+1T/aTdn+1Jndτn+1+aτndnJndτn+aτn−12nT/aTdn−1Jndτn−1+aτn−2n(n−1)T2/aT2dn−2Jndτn−2−dn+2Jndτn+2=1aτ2dnJndτn+1T[(2n+1)aTaτdn+1Jndτn+1+2naTaτ1a2τdn−1Jndτn−1+1T{n(n−1)aT2aτ21aτ2dn−2Jndτn−2}]
(17)式が一般のnの場合の演算方程式となる。
この方法は非零の初期値を持つn階微分項が独立変数で割られないというメリットがある。さらに、初期条件として非零の値が必ず設定できるのが特徴であり、最初に出てくる非零初期条件
J(k)n=(12)nJn(0)δn,j
(ただし δn,jはクロネッカーのデルタ)を第2の積分器に入力すればよいだけである。
無論欠点を完全に廃せたわけではなく、初期値が0の項の除算 (0÷0) は避けられない。ただし、除算器出力が戻される積分器が1つのみであるため、nが大きいほど多くの積分器を通過し誤差は「ある程度は」平滑低減されることが期待できる。一番最初に示したJ0(x)の演算回路中の時間微分項を用いてJ1(x)を得るような回路では、除算出力から積分器1つ分しか離れてない箇所の信号を引っ張り出すため、誤差が無視できなかった。以上のことから、Binglac と Humo の方法は、その従来手法よりは優れていると考えてよいだろう。
n次の第1種ベッセル関数を発生させるのに必要な積分器はn+3個である。
本方式の詳細については文献 [3] を参照して頂きたい。
微分方程式の独立変数を従属変数(媒介変数)にして、式中の特異点を除去する方法である 。0次および1次のベッセル関数の生成で使える。
原方程式の従属変数をy、独立変数をθとおく。θ≥0の範囲でn次ベッセル関数の解を得ることを考えよう。
θ2d2ydθ2+θdydθ+(θ2−n2)Jn=0
dydθ=xとして、連立1階常微分方程式にすると
{dydθ=xdxdθ=(n2−θ2)y+θxθ2
ここで、計算機上で演算時間τとする変数 t を導入し、真の独立変数とする。θはただtについての関数であるとする。(θとtの関係式はまだ決定されない)
(19)の2式左辺に対し連鎖律を適用し、tについての微分に直すと
{dydt=xdθdtdxdt=(n2−θ2)y+θxθ2dθdt
(20)の関係性から、tとθの関係を
dθdt=−ξθ2
として定めれば特異点が消えて都合がよいことがわかる。ただし、ξは定数。
{dydt=ξθ2xdxdt=(n2−θ2)y+θxθ2ξθ2dθdt=ξθ2
∴ddt(yxθ)=ξ(θ2x(n2−θ2)y+θxθ2)
(21)式の作る三次元空間上の軌道\{y(t)\,,x(t)\,,\theta(t)\}を(x,\theta)平面に投影したものが(19)式の解となる。実定数\xiの正負は解軌道上の点が演算時間t(演算変数としては\tauと表記する)が増大するにあたり軌道上をどちらの方向に動くかを決定する。もし正ならば \theta\rightarrow\infty\,(t\rightarrow\infty)となり、負であるならば \theta\rightarrow 0\,(t\rightarrow\infty)となるよう動く。
この効果は時間換算係数a_{\tau}の符号の正負の効果と同様であり、絶対値は積分時定数に影響するので、\xiは時間換算係数であると解釈することもできる。
\xi=1と定めると
\left\{\begin{array}{ll]}\displaystyle\frac{\mathrm{d}y}{\mathrm{d}t}=\theta^2x\\\\\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}=\left(n^2-\theta^2\right)y+\theta{x}\\\\\displaystyle\frac{\mathrm{d}\theta}{\mathrm{d}t}=\theta^2\end{array}\right.\tag{24}
(22)を演算回路として実現すれば、除算は不要になる。かりに\xiを負にすると、 \theta\rightarrow 0\,(t\rightarrow\infty)という「逆行解」が得られることに注意。逆に考えれば、初期条件ではなく\theta=0以外での境界条件でベッセルの微分方程式を解きたい場合には最適な手法と言えるだろう。
なお、ルジャンドルの微分方程式もこの方法で解ける。
本方式について詳しく知りたい方は文献 [4] [5] を参照して頂きたい。
演算の理論 2.6節で述べた方法の応用である。文献[2]のp.90-p.91でも言及されている。
(1) n\geq{1}のとき
ベッセル微分方程式の右辺の0を(k+1) 階微分項としておく。ただし、非負整数 k は、n=1,2のとき3、n\geq{4}のときnとする。独立変数をtとすると
t^2\frac{\mathrm{d}^2J_n}{\mathrm{d}t^2}+t\frac{\mathrm{d}J_n}{\mathrm{d}t}+\bigl(t^2-n^2\bigr)J_n=-\frac{1}{\mu}\frac{\mathrm{d}^{k+1}J_n}{\mathrm{d}t^{k+1}}\tag{25}
ここで\muはオープンループアンプのゲインであり、非常に大きい値をとるとする。
(2) n=0の場合
原方程式の両辺を独立変数で割った式をもとに、
t\frac{\mathrm{d}^2J_0}{\mathrm{d}t^2}+\frac{\mathrm{d}J_0}{\mathrm{d}t}+t\,J_0=-\frac{1}{\mu}\frac{\mathrm{d}^{3}J_0}{\mathrm{d}t^{3}}\tag{26}
とすればよい。
このように、形式的にk+1階微分項を高利得増幅器のオープンゲインで除した項を導入することで、3次以上の場合でも非零の初期条件を確保しつつ、2階微分項について陽に解くことなく実装できる。演算回路の構成もわかりやすいうえ、除算器が不要であり、なおかつ時間換算係数と電圧換算係数の設定を独立して行えることが利点である。ハイブリッド計算方式に最適な手法と言える。もちろん、ベッセル微分方程式以外にも一般に適用可能な手法である。
最大の利点は、換算係数を変えずに回路から除算器を排することができ、精度のよい演算が行えることにある。係数器の設定の手間も少ない。以上のことから、本方法が2.1および2.2節の手法より優れているといえよう。
この方法の詳細は文献 [6] (US Patent 3652843A) を参照して頂きたい。
※ 以上3つの他にも、独立変数をx=\varepsilon e^{-\beta{t}}\,,(\varepsilon\approx 0,\beta\gg 0)として時間原点近傍でのゼロ除算を回避する近似方法、式をもとに評価関数を設定し、最急降下のプログラムを行う方法などがあるが、ここでは省略する。
まず、Binglacの論文通りに実装した結果を示す。a_{\tau}=1、a_{T}=1/10である。
図3 Binglacの方法によるベッセル関数の演算結果
数値計算による解を緑の点でプロットした。両者を比較すると、n=0,n=1 でかなり精密な解が得られているが、n=2以上では誤差が著しく出ていることがわかる。詳しく調べると、回路中の最高階項に入力される高利得増幅器の波形が、t=0近傍でかなり暴れていることが判明した。対策として、アナログコンピュータにあらかじめ設けて置いたデジタル信号駆動型のアナログスイッチ回路と、比較器の2値出力信号を利用して、演算開始後一定期間は最高階積分器の入力はJ_nのみ帰還させ、高利得増幅器(除算器)の出力はその後に投入するような結線プログラムを試みたが、多少の改善は見られたものの実用には至らなかった。
以上の結果を鑑み、より精度の良い結果を得るため従来の除算器を最急降下法を利用した除算回路に置き換えることとした。
a_{T}=1/10, a_{\tau}=1 としたときの、最急降下除算を利用した3次ベッセル微分方程式の演算回路を図に示す。
図4 n=3の場合の演算回路
除算の初期値は別個に設定する必要がある。私は級数展開して極限をとり初期値を計算した。最急降下のゲインは数学的には大きくとった方が収束がよいのだが、計算機回路的な制約上、大きすぎるのも誤差の元となるから最適値を探る必要がある。
n=0 から n=3 の場合の低速モードにおける演算結果を次に示す。
図5 演算結果 (Binglac + 最急降下法による除算)
n=3 まではかなり精確な解を得られた。しかし、より高次の場合は時間が経つと誤差が無視できなくなった。
回路からわかる通り、3次以上では解関数を開ループによる単純縦続積分で求めるので、その積分で定常誤差が累積して現れてしまうことが原因と思われる。
特許文献[6]の方法に忠実に従った場合の結果を示す。
図6 黒川らの方法によるベッセル関数の演算結果
0次と1次のベッセル関数については正確な解を出しているが、n=3の波形はところどころに歪が生じている。こちらも、高利得増幅器の出力が超低周波で振動していたことが直接の原因である。
そこで文献[8]を参考に、高利得増幅器の利得を係数器で微調整できるようにしてみた。アンプのオープンゲインが有限であることに起因する誤差を補正することが目的である。その場合の \ddot{y(t)}+y(t)=0を解くサークルテスト回路と、高速モードにおける演算結果を図に示す。
図7 サークルテスト
係数器の \gamma の値を0.5前後で微調節する。
以上の補正を踏まえた3次のベッセル微分方程式演算回路は以下のようになる。引き続き a_{\tau}=1,\,a_{T}=0.1 とした。
図8 n=3の場合の演算回路(黒川らによる方法)
n=0からn=3の演算結果を以下に示す。
図9 黒川らの方法による演算結果(利得補正あり)
少なくとも3次までは高精度な解が得られる。
しかし、これ以上高次の解では誤差が目立つ。Binglacの方法と異なり開ループ積分による誤差は生じえないが、次数が大きくなるにしたがい、ループ内の縦続積分器個数が多くなるうえに、 t=0 で非零の値をとる積分項が1つだけという条件であるから、回路中唯一の安定負帰還ループ(図8中でP10を含むメインループ)の帰還量が演算開始後しばらくたったあとも小さいままであり、そのあいだに微小擾乱が拡大することが原因と考える。
[1] 森口繁一, 宇田川銈久, 一松信.「岩波全書 数学公式Ⅲ-特殊関数-」第1版, 岩波書店, 1975.
[2] 藤田広一. 「アナログ計算機のプログラム」第12版, 昭晃堂, 1978.
[3] S. P. Binglac, E. A. Humo. "Analog Computer Generation of Bessel Functions of Arbitrary Order" IEEE TRANSACTIONS ON ELECTRONIC COMPUTERS. VOL. EC-14, NO. 6, 1965.
[4] A. Hausner, "Parametric Techniques for Eliminating Division and Treating Singularities in Computer Solutions of Ordinary Differential Equations," in IRE Transactions on Electronic Computers, vol. EC-11, no. 1, pp. 42-45, Feb. 1962.
[5] T. A. Newton, "Some Parametric Techniques in the Analog Solution of Ordinary Differential Equations," in IEEE Transactions on Computers, vol. C-24, no. 1, pp. 1-8, Jan. 1975.
[6] Kazuo Kurokawa, Ikuo Matsuda. Setup system in analog computer. US Patent 3652843A. Mar. 28, 1972. https://patents.google.com/patent/US3652843
[7] 森 繁雄, 芹沢 正三, 安藤 四郎. 「工業基礎数学 微分方程式」第10版, 東京書籍 ,1986.
[8] 松田 郁夫, 黒川 一夫. 「ハイブリッド計算システムにおける計算方式 ーハイブリッド計算機としてのアナログ演算部の計算方式ー」電気試験所彙報, 1969, 第33巻1号, p.22-p.35