SAD部屋 - Ticket-26 Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
! Bessel関数の不具合
:Priority:Normal
:Reporter:Akio Morita
:Status:Assigned
:Assigned to:Akio Morita
:Version:3648
:Milestone:?
:Created:2011-05-23
!! Description
#BesselJ及び BesselYの実軸 $x\sim-20$付近で不連続がある
#実数次のBesselJ/Y/I/Kにて、実軸 $x<0$での虚部が 0である
#次数が非整数の実数であるBesselJ/Y/I/Kで実軸上の値の虚部が 0である
#BesselKにて、実軸 $|x| \gg 0$領域で、NaNを返す
!! Attachement files
{{attach_map('Ticket-26')}}
{{attach_form('','')}}
!! Changelog
!!!Akio Morita (2011-05-23 (月) 11:25:32)
(1)については、$z\sim\infty$近傍での$\frac{1}{z}$展開を適用する際に$Re[z]>0$の公式を$Re[z]<0$の領域にも誤って適用し、誤った位相で計算しているのが原因。
r3649にて、cbesasym()内で$Re[z]<0$の場合は $z$の符号変換公式を適用して修正。
!!!Akio Morita (2011-05-23 (月) 11:36:22)
(2,3)については、$\nu\in\mathbf{R}$かつ$z\in\mathbf{R}$の際に、tfbessel()がReal型を返すのが原因。
r3651にて$z<0$の領域にて$\nu\notin\mathbf{Z}$もしくはBesselY, BesselKの場合は、Complex型で返すよう修正。
また、BesselI/Kについては、cbesseli(), cbesselk()も同様の処理が有るので、
r3652にて実数値を返す条件を$Im[z] = 0$かつ$\nu\in\mathbf{Z}$に変更
!!!Akio Morita (2011-05-23 (月) 11:42:46)
(4)については、実軸上$|z|\gg0$領域で発散関数同士を減算することにより
(INF - INF)による NaNと思われる。中間領域では桁落ちも有るので、
BesselKの展開公式を直接適用すべきと思われる。
!!!Akio Morita (2011-05-23 (月) 11:48:13)
*$BesselJ_\nu(z)$の公式
$BesselJ_\nu(-z) = (-z)^\nu z^{-\nu}BesselJ_\nu(z)$
$BesselJ_{-n}(z) = (-1)^nBesselJ_n(z)\;(n\in\mathbf{Z})$
$BesselJ_\nu(z) = \frac{1}{2\pi}\int^{\pi}_{-\pi}\cos(z\sin{\theta} - n\theta)d\theta\;(\nu\in\mathbf{Z})$
周期関数の一周期積分なので、中点公式、台形公式で極めて効率よく積分出来る
$BesselJ_\nu(z) = \frac{(\frac{z}{2})^\nu}{\Gamma(\nu + \frac{1}{2})\Gamma(\frac{1}{2})}\int_{0}^{\pi}\cos(z\cos{\theta})\sin^{2\nu}{\theta}d\theta\;(Re[\nu] > 0)$
*$BesselY_\nu(z)$の公式
$BesselY_\nu(z) = \frac{\cos(\nu\pi)BesselJ_\nu(z) - BesselJ_{-\nu}(z)}{\sin(\nu\pi)}$
$BesselY_\nu(-z) = (-z)^\nu z^{-\nu}BesselY_\nu(z) + 2 (-z)^{\frac{1}{2}} z^{-\frac{1}{2}}BesselJ_{-\nu}(z)$
$BesselY_{-\nu}(z) = \cos(\nu\pi)BesselY_\nu(z) + \sin(\nu\pi)BesselJ_\nu(z)$
*$BesselI_\nu(z)$の公式
$BesselI_\nu(z) = (iz)^{-\nu}z^{\nu}BesselJ_\nu(iz)$
$BesselI_\nu(-z) = (-z)^\nu z^{-\nu}BesselI_\nu(z)$
*$BesselK_\nu(z)$の公式
$BesselK_\nu(z) = \frac{\pi}{2}\frac{BesselI_{-\nu}(z) - BesselI_{\nu}(z)}{\sin(\nu\pi)}$
$BesselK_\nu(-z) = (-z)^{-\nu} z^{\nu}BesselK_\nu(z) + \pi (-z)^{-\frac{1}{2}} z^{\frac{1}{2}}BesselI_\nu(z)$
$BesselK_{-\nu}(z) = BesselK_\nu(z)$
$BesselK_\nu(z) = \frac{\sqrt{\pi}}{\Gamma(\nu+\frac{1}{2})}(\frac{z}{2})^\nu\int_{0}^{\infty}e^{-z\cosh(t)}\sinh^{2\nu}(t)dt\;(Re[\nu]>-\frac{1}{2} \wedge Re[z] > 0)$
!!!Akio Morita (2011-05-23 (月) 14:50:17)
$z$に依存する$\nu$回転因子について補足
*偏角について次の定義で枝を選択
**$\forall z \in \mathbf{C}, \pi \geq Arg[z] > -\pi$
**$-1 = e^{i\pi}$
*$(-z)^\nu z^{-\nu}$
$(-z)^\nu z^{-\nu} = e^{+i\pi\nu}\;(Im[z] < 0 \vee Im[z] = 0 \wedge Re[z] > 0)$
$(-z)^\nu z^{-\nu} = e^{-i\pi\nu}\;(Im[z] > 0 \vee Im[z] = 0 \wedge Re[z] < 0)$
*$(iz)^{-\nu}z^{\nu}$
$(iz)^{-\nu}z^{\nu} = e^{+i\frac{3\pi\nu}{2}}\;(Re[z] < 0 \wedge Im[z] \geq 0)$
$(iz)^{-\nu}z^{\nu} = e^{-i\frac{\pi\nu}{2}}\;(Re[z] \geq 0 \vee Im[z] < 0)$
!!!Akio Morita (2011-05-24 (火) 17:09:39)
r3657にて、Z=0近傍に関しては、複素平面上の枝を含めて修正完了
残る問題は、
*$Re[BesselK_\nu(z)]$が $|z| \gg 1$で 0を保持しない($z\in\mathbf{R}, \nu\in\mathbf{Z}$)
*$|z| \gg 1$領域ので$BesselK_\nu(z)$の有効数字($BesselI_\nu(z)$が発散関数なので、大数同士の減算により桁落ちが発生する)
*$|z| < 20$領域での近似精度
!!!Akio Morita (2011-05-26 (木) 11:58:44)
$BesselK_\nu(z)$の$|z|\gg1$領域の近似に関しては、発散関数同士の演算になるので、元の近似関数の発散項を式レベルで処理する必要がある。
$|Arg[z]| < \pi$かつ$|z|\gg1$の場合、多項式$P_\nu(z)$, $Q_\nu(z)$を用いた以下の近似が成立する
$BesselJ_\nu(z) \sim \sqrt{\frac{2}{\pi z}} [P_\nu(\frac{8}{z})\cos(z_\nu) - Q_\nu(\frac{8}{z})\sin(z_\nu)]$
$BesselY_\nu(z) \sim \sqrt{\frac{2}{\pi z}} [P_\nu(\frac{8}{z})\sin(z_\nu) + Q_\nu(\frac{8}{z})\cos(z_\nu)]$
$z_\nu = z - \frac{\pi}{2}\nu - \frac{\pi}{4}$
$|Arg[z]| \sim \pi$近傍の収束が悪いので注意。
$BesselK_\nu(iz)$の定義に現れる項のうち
$BesselY_\nu(iz) + i BesselJ_\nu(iz)$は
$\sqrt{\frac{2}{\pi z}}[P_\nu(\frac{8}{z}) + i Q_\nu(\frac{8}{z})]e^{iz_\nu}$と書ける
!!!Akio Morita (2011-05-26 (木) 19:41:45)
r3662にて、$BesselK_\nu(z)$の$|z|\gg1$領域での振る舞いを修正
* $|z|\gg1$領域での実部の有効数字を改善
* $z$の実軸上では、実部は無限遠で減衰関数の形を保つ(発散しない)
数学関数的に以下のケースで問題が発生します
* $z\rightarrow-\infty$でのBesselKの虚部がNaNになる
* $z$の一般の偏角で、$|z|\rightarrow\infty$の際にNaNが発生する
発生原因としては、無限大を含む複素数型に対してのComplex[0, 1]の乗算などがあげられます。
多くは、BesselJ, BesselY, BesselIの内部で発生しているので、大域的な修正を必要とします
!!!Akio Morita (2011-05-28 (土) 04:46:30)
*$BesselY_\nu(z)$にて、$|z|\sim5$付近の近似精度が悪い($10^{-7}$程度)
*$BesselJ_\nu(z)$にて、$|z|<20$での近似精度が、$10^{-15}$程度で改善の余地あり
!!!Akio Morita (2011-05-28 (土) 12:07:37)
BesselJに関しては、逆進漸化式の精度評価(libm jn参照)の修正と、
BesselJの和則を計算する際の桁落ち対策が必要。
$|z|\ll0$での級数展開も必要か？
{{its_edit_ticket_form}}
:Priority:Normal
:Reporter:Akio Morita
:Status:Assigned
:Assigned to:Akio Morita
:Version:3648
:Milestone:?
:Created:2011-05-23
!! Description
#BesselJ及び BesselYの実軸 $x\sim-20$付近で不連続がある
#実数次のBesselJ/Y/I/Kにて、実軸 $x<0$での虚部が 0である
#次数が非整数の実数であるBesselJ/Y/I/Kで実軸上の値の虚部が 0である
#BesselKにて、実軸 $|x| \gg 0$領域で、NaNを返す
!! Attachement files
{{attach_map('Ticket-26')}}
{{attach_form('','')}}
!! Changelog
!!!Akio Morita (2011-05-23 (月) 11:25:32)
(1)については、$z\sim\infty$近傍での$\frac{1}{z}$展開を適用する際に$Re[z]>0$の公式を$Re[z]<0$の領域にも誤って適用し、誤った位相で計算しているのが原因。
r3649にて、cbesasym()内で$Re[z]<0$の場合は $z$の符号変換公式を適用して修正。
!!!Akio Morita (2011-05-23 (月) 11:36:22)
(2,3)については、$\nu\in\mathbf{R}$かつ$z\in\mathbf{R}$の際に、tfbessel()がReal型を返すのが原因。
r3651にて$z<0$の領域にて$\nu\notin\mathbf{Z}$もしくはBesselY, BesselKの場合は、Complex型で返すよう修正。
また、BesselI/Kについては、cbesseli(), cbesselk()も同様の処理が有るので、
r3652にて実数値を返す条件を$Im[z] = 0$かつ$\nu\in\mathbf{Z}$に変更
!!!Akio Morita (2011-05-23 (月) 11:42:46)
(4)については、実軸上$|z|\gg0$領域で発散関数同士を減算することにより
(INF - INF)による NaNと思われる。中間領域では桁落ちも有るので、
BesselKの展開公式を直接適用すべきと思われる。
!!!Akio Morita (2011-05-23 (月) 11:48:13)
*$BesselJ_\nu(z)$の公式
$BesselJ_\nu(-z) = (-z)^\nu z^{-\nu}BesselJ_\nu(z)$
$BesselJ_{-n}(z) = (-1)^nBesselJ_n(z)\;(n\in\mathbf{Z})$
$BesselJ_\nu(z) = \frac{1}{2\pi}\int^{\pi}_{-\pi}\cos(z\sin{\theta} - n\theta)d\theta\;(\nu\in\mathbf{Z})$
周期関数の一周期積分なので、中点公式、台形公式で極めて効率よく積分出来る
$BesselJ_\nu(z) = \frac{(\frac{z}{2})^\nu}{\Gamma(\nu + \frac{1}{2})\Gamma(\frac{1}{2})}\int_{0}^{\pi}\cos(z\cos{\theta})\sin^{2\nu}{\theta}d\theta\;(Re[\nu] > 0)$
*$BesselY_\nu(z)$の公式
$BesselY_\nu(z) = \frac{\cos(\nu\pi)BesselJ_\nu(z) - BesselJ_{-\nu}(z)}{\sin(\nu\pi)}$
$BesselY_\nu(-z) = (-z)^\nu z^{-\nu}BesselY_\nu(z) + 2 (-z)^{\frac{1}{2}} z^{-\frac{1}{2}}BesselJ_{-\nu}(z)$
$BesselY_{-\nu}(z) = \cos(\nu\pi)BesselY_\nu(z) + \sin(\nu\pi)BesselJ_\nu(z)$
*$BesselI_\nu(z)$の公式
$BesselI_\nu(z) = (iz)^{-\nu}z^{\nu}BesselJ_\nu(iz)$
$BesselI_\nu(-z) = (-z)^\nu z^{-\nu}BesselI_\nu(z)$
*$BesselK_\nu(z)$の公式
$BesselK_\nu(z) = \frac{\pi}{2}\frac{BesselI_{-\nu}(z) - BesselI_{\nu}(z)}{\sin(\nu\pi)}$
$BesselK_\nu(-z) = (-z)^{-\nu} z^{\nu}BesselK_\nu(z) + \pi (-z)^{-\frac{1}{2}} z^{\frac{1}{2}}BesselI_\nu(z)$
$BesselK_{-\nu}(z) = BesselK_\nu(z)$
$BesselK_\nu(z) = \frac{\sqrt{\pi}}{\Gamma(\nu+\frac{1}{2})}(\frac{z}{2})^\nu\int_{0}^{\infty}e^{-z\cosh(t)}\sinh^{2\nu}(t)dt\;(Re[\nu]>-\frac{1}{2} \wedge Re[z] > 0)$
!!!Akio Morita (2011-05-23 (月) 14:50:17)
$z$に依存する$\nu$回転因子について補足
*偏角について次の定義で枝を選択
**$\forall z \in \mathbf{C}, \pi \geq Arg[z] > -\pi$
**$-1 = e^{i\pi}$
*$(-z)^\nu z^{-\nu}$
$(-z)^\nu z^{-\nu} = e^{+i\pi\nu}\;(Im[z] < 0 \vee Im[z] = 0 \wedge Re[z] > 0)$
$(-z)^\nu z^{-\nu} = e^{-i\pi\nu}\;(Im[z] > 0 \vee Im[z] = 0 \wedge Re[z] < 0)$
*$(iz)^{-\nu}z^{\nu}$
$(iz)^{-\nu}z^{\nu} = e^{+i\frac{3\pi\nu}{2}}\;(Re[z] < 0 \wedge Im[z] \geq 0)$
$(iz)^{-\nu}z^{\nu} = e^{-i\frac{\pi\nu}{2}}\;(Re[z] \geq 0 \vee Im[z] < 0)$
!!!Akio Morita (2011-05-24 (火) 17:09:39)
r3657にて、Z=0近傍に関しては、複素平面上の枝を含めて修正完了
残る問題は、
*$Re[BesselK_\nu(z)]$が $|z| \gg 1$で 0を保持しない($z\in\mathbf{R}, \nu\in\mathbf{Z}$)
*$|z| \gg 1$領域ので$BesselK_\nu(z)$の有効数字($BesselI_\nu(z)$が発散関数なので、大数同士の減算により桁落ちが発生する)
*$|z| < 20$領域での近似精度
!!!Akio Morita (2011-05-26 (木) 11:58:44)
$BesselK_\nu(z)$の$|z|\gg1$領域の近似に関しては、発散関数同士の演算になるので、元の近似関数の発散項を式レベルで処理する必要がある。
$|Arg[z]| < \pi$かつ$|z|\gg1$の場合、多項式$P_\nu(z)$, $Q_\nu(z)$を用いた以下の近似が成立する
$BesselJ_\nu(z) \sim \sqrt{\frac{2}{\pi z}} [P_\nu(\frac{8}{z})\cos(z_\nu) - Q_\nu(\frac{8}{z})\sin(z_\nu)]$
$BesselY_\nu(z) \sim \sqrt{\frac{2}{\pi z}} [P_\nu(\frac{8}{z})\sin(z_\nu) + Q_\nu(\frac{8}{z})\cos(z_\nu)]$
$z_\nu = z - \frac{\pi}{2}\nu - \frac{\pi}{4}$
$|Arg[z]| \sim \pi$近傍の収束が悪いので注意。
$BesselK_\nu(iz)$の定義に現れる項のうち
$BesselY_\nu(iz) + i BesselJ_\nu(iz)$は
$\sqrt{\frac{2}{\pi z}}[P_\nu(\frac{8}{z}) + i Q_\nu(\frac{8}{z})]e^{iz_\nu}$と書ける
!!!Akio Morita (2011-05-26 (木) 19:41:45)
r3662にて、$BesselK_\nu(z)$の$|z|\gg1$領域での振る舞いを修正
* $|z|\gg1$領域での実部の有効数字を改善
* $z$の実軸上では、実部は無限遠で減衰関数の形を保つ(発散しない)
数学関数的に以下のケースで問題が発生します
* $z\rightarrow-\infty$でのBesselKの虚部がNaNになる
* $z$の一般の偏角で、$|z|\rightarrow\infty$の際にNaNが発生する
発生原因としては、無限大を含む複素数型に対してのComplex[0, 1]の乗算などがあげられます。
多くは、BesselJ, BesselY, BesselIの内部で発生しているので、大域的な修正を必要とします
!!!Akio Morita (2011-05-28 (土) 04:46:30)
*$BesselY_\nu(z)$にて、$|z|\sim5$付近の近似精度が悪い($10^{-7}$程度)
*$BesselJ_\nu(z)$にて、$|z|<20$での近似精度が、$10^{-15}$程度で改善の余地あり
!!!Akio Morita (2011-05-28 (土) 12:07:37)
BesselJに関しては、逆進漸化式の精度評価(libm jn参照)の修正と、
BesselJの和則を計算する際の桁落ち対策が必要。
$|z|\ll0$での級数展開も必要か？
{{its_edit_ticket_form}}