phys_552.bessel#
Some utilities for computing properties of the Bessel functions for the DVR basis.
Functions#
|
Return the d'th derivative of sinc(x) = sin(x)/x. |
|
Return the d'th derivative of the bessel functions |
|
Return the first N positive roots of the bessel function |
|
Return a function that computes the d'th derivative of |
Module Contents#
- sinc(x, d=0)#
Return the d’th derivative of sinc(x) = sin(x)/x.
- Parameters:
x ({float, array}) – Argument
d (int, optional) – Order of derivative.
Examples
>>> N = 3 >>> np.allclose(sinc(2.0), np.sin(2.0)/2.0) True >>> print("skip ...");from mmfutils.math.differentiate import differentiate skip ... >>> x = np.array([-1,0,0.5,1]) >>> np.allclose(differentiate(lambda x:sinc(x), x), ... sinc(x, d=1)) True
- J(nu, d=0)#
Return the d’th derivative of the bessel functions \(J_{\nu}(z)\).
Examples
>>> J0 = J(0.5); J1 = J(1.5); J2 = J(2.5); >>> z = 2.5; nu = 1.5 >>> abs(J0(z) + J2(z) - 2*nu/z*J1(z)) < _EPS True
- j_root(nu, N, rel_tol=2 * _EPS)#
Return the first N positive roots of the bessel function J_nu(x).
- Parameters:
Notes
The general method is to first estimate the roots with a bisection/secant method, and then polish them using Newton’s method.
We start by estimating the lower bound for the first for non-negative \(\nu\)
\[\nu + \nu^{1/3} < j_{\nu,1} < \nu + 1.85575 \nu^{1/3} + \pi\]Then, using the fact that the roots are spaced by at least \(\pi\), we step through the sign changes to bracket all of the desired roots.
lowest root using the following heuristics
\[\begin{split}j_{\nu,s} \approx \begin{cases} 2\sqrt{\nu + 1} & -1 < \nu < -0.8\\ \left(\frac{\nu}{2} + \frac{3}{4}\right)\pi & -0.8 < \nu < 2.5\\ \nu + 1.85575 \nu^{1/3} & 2.5 \leq \nu \end{cases}\end{split}\]Examples
>>> nu = 2.5 >>> j_ = j_root(nu, 2000) >>> J_ = J(nu)(j_)
These are roots!
>>> np.max(abs(J_/j_)) < _EPS True
They are also distinct
>>> pi < min(np.diff(j_)) True
And the spacing is decreasing, meaning we have not skipped any.
>>> np.max(np.diff(np.diff(j_))) < 0 True
- J_sqrt_pole(nu, zn, d=0)#
Return a function that computes the d’th derivative of sqrt(z)*J(nu,z)/(z - zn) where zn is a root: J(nu, zn) = 0. This combination appears in the Bessel function DVR basis.
Notes
\[\frac{\sqrt{z}J_{\nu}(z)}{z - z_{n}}\]As \(z\) approaches \(z_n\), this has the form of 0/0, so one can apply a form of l’Hopital’s rule to reduce the round-off error. The specified form of the function has been chosen for special properties of the Bessel functions. Express the function as
\[\begin{split}F(z) &= \frac{f(z)}{z - z_n} = \frac{\sqrt{z}J_{\nu}(z)}{z - z_n}\\ F'(z) &= \frac{f'(z)}{z - z_n} - \frac{f(z)}{(z - z_n)^2}\end{split}\]Let \(\delta = z - z_n\). Close to the singular point we use the Taylor series:
\[\sum_{m=0}^{\infty}\frac{a_m\delta^{m}}{m!}\]\[\begin{split}F(z) &= f'(z_n) + \sum_{m=3}^{\infty}\frac{f^{(m)}(z_n)\delta^{m-1}}{m!} = \sum_{m=0}^{\infty}\frac{f^{(m+1)}(z_n)\delta^{m}}{(m+1)m!}\\ a_m &= \frac{f^{(m+1)}}{m+1}\\ F'(z) &= \sum_{m=3}^{\infty}\frac{(m-1)f^{(m)}(z_n)\delta^{m-2}}{m!} = \sum_{m=1}^{\infty}\frac{f^{(m+2)}(z_n)\delta^{m}}{(m+2)m!}\\ a_m &= \frac{f^{(m+2)}}{m+2}\\\end{split}\]The first few derivatives are presented here:
\[\begin{split}f(z) &= \sqrt{z}J_{\nu}(z)\\ f'(z) &= \frac{J_{\nu}(z)}{2\sqrt{z}} + \sqrt{z}J'_{\nu}(z) = \frac{f(z)}{2z} + \sqrt{z}J'_{\nu}(z)\\ f''(z) &= \sqrt{z}J_{\nu}(z)\left( \frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right) = f(z)\left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)\\ f'''(z) &= f'(z)\left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right) - 2f(z)\frac{\nu^2 - \tfrac{1}{4}}{z^3}\\ f^{(4)}(z) &= f(z)\left[ \left(\frac{\nu^2 - \tfrac{1}{4}}{z^2} - 1\right)^2 + 6\frac{\nu^2 - \tfrac{1}{4}}{z^4}\right] - 4f'(z)\frac{\nu^2 - \tfrac{1}{4}}{z^3}\end{split}\]Evaluated at the root \(z=z_n\) these become:
\[\begin{split}f(z_{n}) &= 0\\ f'(z_{n}) &= \sqrt{z_{n}}J'_{\nu}(z_{n})\\ f''(z_{n}) &= 0\\ f'''(z_{n}) &= f'(z_{n})\left( \frac{\nu^2 - \tfrac{1}{4}}{z_{n}^2} - 1\right)\\ f^{(4)}(z_{n}) &= - 4f'(z_{n})\frac{\nu^2 - \tfrac{1}{4}}{z_{n}^3}\end{split}\]with both the function and the second derivative vanishing.
To determine where to use this formula, we match the estimate roundoff error with the truncation error. The Bessel functions are of order unity and are typically calculated to an absolute accuracy of \(\epsilon\). The round-off error in the numerator is \(\epsilon f(z)\) and \(\epsilon \sqrt{2} z_n\) in the denominator. The roundoff errors in the denominator dominate both cases:
\[\begin{split}\delta F(z) &\sim \epsilon \frac{\sqrt{2}z_n F(z)}{\delta} \sim \frac{\sqrt{2}\epsilon z_n f(z)}{\delta^2} \sim \frac{\sqrt{2}\epsilon z_n f'(z_n)}{\delta}\\ \delta F'(z) &\sim \frac{2\epsilon z_n f(z)}{\delta^3} \sim \frac{2\epsilon z_n f'(z_n)}{\delta^2}\end{split}\]To choose the appropriate transition point, we equate half of this with the truncation error to transition points:
\[\begin{split}\delta_c &\sim \left( \frac{72\epsilon z_n f'(z_n)}{\sqrt{2}f^{(4)}(z_n)} \right)^{1/4} \sim \left(\frac{72\epsilon z_n}{\sqrt{2}} \right)^{1/4}\\ \delta_c' &\sim \left(120\epsilon z_n\right)^{1/5}\end{split}\]the fact that :math:f(z) behaves asymptotically as a \(\sqrt{2/\pi}\cos(z + \phi)\) and so all derivatives have essentially the same magnitude.
Examples
>>> nu = 5.5 >>> zn = j_root(nu,21)[-1] >>> abs(zn - 73.62361318251753391646) < 1e-16 True >>> float(J_sqrt_pole(nu,zn)(zn)) -0.796778576780013...
-0.796778576780013129760
You can also use a vector of zn, but only if it is commensurate with the argument:
>>> zn = j_root(nu,21) >>> float(J_sqrt_pole(nu, zn)(zn[-1])[20]) -0.796778576780013...