pub(crate) fn j1f_asympt_alpha(x: f64) -> f64Expand description
Note expansion generation below: this is negative series expressed in Sage as positive,
so before any real evaluation x=1/x should be applied.
Generated by SageMath:
def binomial_like(n, m):
prod = QQ(1)
z = QQ(4)*(n**2)
for k in range(1,m + 1):
prod *= (z - (2*k - 1)**2)
return prod / (QQ(2)**(2*m) * (ZZ(m).factorial()))
R = LaurentSeriesRing(RealField(300), 'x',default_prec=300)
x = R.gen()
def Pn_asymptotic(n, y, terms=10):
# now y = 1/x
return sum( (-1)**m * binomial_like(n, 2*m) / (QQ(2)**(2*m)) * y**(QQ(2)*m) for m in range(terms) )
def Qn_asymptotic(n, y, terms=10):
return sum( (-1)**m * binomial_like(n, 2*m + 1) / (QQ(2)**(2*m + 1)) * y**(QQ(2)*m + 1) for m in range(terms) )
P = Pn_asymptotic(1, x, 50)
Q = Qn_asymptotic(1, x, 50)
R_series = (-Q/P)
# alpha is atan(R_series) so we're doing Taylor series atan expansion on R_series
arctan_series_Z = sum([QQ(-1)**k * x**(QQ(2)*k+1) / RealField(700)(RealField(700)(2)*k+1) for k in range(25)])
alpha_series = arctan_series_Z(R_series)
# see the series
print(alpha_series)See notes/bessel_asympt.ipynb for generation