Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve bessel_j in transition region #2137

Merged
merged 3 commits into from
Jan 7, 2025

Conversation

fredrik-johansson
Copy link
Collaborator

  • Add acb_hypgeom_bessel_j_deriv_bound
  • Mitigate cancellation for $J_n(z)$ in the transition region and on wide intervals, by using increased precision + derivative bounds.
  • Essentially only implemented for small integer $\nu$, though evaluation should also be improved in some other cases at least when both $\nu$ and $z$ exact.

Some numerical examples:

>>> sum(CC.bessel_j(0, k) for k in range(100))
[1.41921 +/- 3.56e-6]                # before
[1.419207859380 +/- 2.35e-13]        # after

>>> w = CC.exp_pi_i(QQ(1)/100); sum(CC.bessel_j(0, w*k) for k in range(101))
([0.7803 +/- 4.35e-5] + [0.6276 +/- 6.52e-5]*I)                   # before
([0.78030446659 +/- 2.44e-12] + [0.62756344686 +/- 5.49e-12]*I)   # after

>>> CC.bessel_j(1, "10 +/- 5")
[+/- 1.75e+5]    # before
[+/- 0.370]      # after

>>> CC.bessel_j(1, "100 +/- 5")
[+/- 0.361]      # before
[+/- 0.0848]     # after

>>> CC.bessel_j(1, "10 +/- 0.01")
[+/- 7.77]          # before
[0.04 +/- 6.00e-3]  # after

>>> CC.bessel_j(0, "10 +/- 100")
FlintUnableError: ...   # before
[+/- 1.01]              # after

>>> CC.bessel_j(0, CC("10 +/- 100") * CC.i())
FlintUnableError: ...   # before
[+/- 5.93e+47]          # after

>>> CC.bessel_j(0.5+0.75j, 26)
([0.20772 +/- 3.94e-6] + [-0.14801 +/- 2.12e-6]*I)                        # before
([0.207717454673942 +/- 6.42e-16] + [-0.148010904937155 +/- 5.17e-16]*I)  # after

Added an example integral:

# before
fredrik@mordor:~/src/flint$ build/examples/integrals -i 38 -verbose
I38 = int_0^{1000} J_0(x) J_1(x) J_2(x) dx ...
depth 6/128, eval 4126/68096, 105 leaf intervals
cpu/wall(s): 0.168 0.167
I38 = [0.057674 +/- 1.03e-7]

# after
fredrik@mordor:~/src/flint$ build/examples/integrals -i 38 -verbose
I38 = int_0^{1000} J_0(x) J_1(x) J_2(x) dx ...
depth 5/128, eval 1688/68096, 16 leaf intervals
cpu/wall(s): 0.058 0.058
I38 = [0.057673961810863 +/- 1.35e-16]

This wants a lot of fine-tuning and something better for non-integer $\nu$, but one has to start somewhere...

@fredrik-johansson fredrik-johansson merged commit 151f514 into flintlib:main Jan 7, 2025
12 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant