sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.44k stars 480 forks source link

Incorrect output(?) from asymptotic_expansions.SingularityAnalysis #33994

Open mezzarobba opened 2 years ago

mezzarobba commented 2 years ago

I believe that the coefficient of n^(-3/2) in

sage: asy = asymptotic_expansions.SingularityAnalysis('n', alpha=1/2, beta=2, normalized=False, precision=6); asy
1/sqrt(pi)*n^(-1/2)*log(n)^2 + (2*(euler_gamma + 2*log(2))/sqrt(pi))*n^(-1/2)*log(n) + (-1/2*pi^(3/2) + (euler_gamma^2 + 4*euler_gamma*log(2) + 4*log(2)^2)/sqrt(pi))*n^(-1/2) - 1/8/sqrt(pi)*n^(-3/2)*log(n)^2 + (1/4*(3*euler_gamma + 6*log(2) - 8)/sqrt(pi) - (euler_gamma + 2*log(2) - 2)/sqrt(pi))*n^(-3/2)*log(n) + (-1/24*(48*euler_gamma - 9*euler_gamma^2 - 12*(3*euler_gamma - 8)*log(2) - 36*log(2)^2 - 64)/sqrt(pi) + 1/2*(4*euler_gamma - euler_gamma^2 - 4*(euler_gamma - 2)*log(2) - 4*log(2)^2 - 4)/sqrt(pi) - 1/48*(9*pi^2 + 64)/sqrt(pi) + 1/4*(pi^2 - 8)/sqrt(pi))*n^(-3/2) + O(n^(-5/2)*log(n)^2)
sage: asy.parent().change_parameter(coefficient_ring=CBF)(asy)
[0.564189583547756 +/- 3.91e-16]*n^(-1/2)*log(n)^2 + [2.21558380774574 +/- 3.71e-15]*n^(-1/2)*log(n) + [-0.60900348841611 +/- 4.74e-15]*n^(-1/2) + [-0.0705236979434695 +/- 6.27e-17]*n^(-3/2)*log(n)^2 + [-0.27694797596822 +/- 3.43e-15]*n^(-3/2)*log(n) + [-1.4283801200753 +/- 4.21e-14]*n^(-3/2) + O(n^(-5/2)*log(n)^2)

is incorrect. Bruno Salvy's equivalent (in Maple) confirms and gives c + 11/(3√π) (≈ 0.64) instead, where c ≈ -1.42 is the value in the above output.

CC: @cheuberg @behackl @dkrenn

Component: asymptotic expansions

Issue created by migration from https://trac.sagemath.org/ticket/33994

mezzarobba commented 2 years ago
comment:1
sage: z = SR.var('z')
sage: expr = 1/(1-z)^(1/2)*log(1/(1-z))^2
sage: ref = list(expr.series(z, 1000).truncate().polynomial(QQ))
sage: asy = asymptotic_expansions.SingularityAnalysis('n', alpha=1/2, beta=2, normalized=False, precision=6)
sage: ini = asy.truncate(5).exact_part().symbolic_expression()
sage: [(RBF(ref[n] - ini(n=n))*RBF(n)**(3/2)) for n in range(980, 1000)]
[[0.640997960 +/- 4.69e-10],
 [0.640997394 +/- 2.43e-10],
 [0.640996829 +/- 1.62e-10],
 [0.640996265 +/- 3.37e-10],
 [0.640995702 +/- 4.93e-10],
 [0.640995139 +/- 5.73e-10],
 [0.640994578 +/- 4.68e-10],
 [0.640994018 +/- 3.68e-10],
 [0.640993459 +/- 2.95e-10],
 [0.640992901 +/- 2.11e-10],
 [0.640992344 +/- 1.54e-10],
 [0.640991788 +/- 9.93e-11],
 [0.640991233 +/- 1.51e-10],
 [0.640990679 +/- 2.13e-10],
 [0.640990126 +/- 2.71e-10],
 [0.640989574 +/- 3.37e-10],
 [0.640989023 +/- 4.21e-10],
 [0.640988473 +/- 5.19e-10],
 [0.640987923 +/- 5.86e-10],
 [0.640987375 +/- 4.60e-10]]
mezzarobba commented 2 years ago
comment:2

oops