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.

Sign up for GitHub

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
Merged

Improve bessel_j in transition region #2137

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0796831
Better bessel_j in transition region and on wide intervals
fredrik-johansson Jan 7, 2025
9ec49c4
more test code
fredrik-johansson Jan 7, 2025
986bb94
restore evaluation at infinity
fredrik-johansson Jan 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions doc/source/acb_hypgeom.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -445,6 +445,11 @@ Error functions and Fresnel integrals
Bessel functions
-------------------------------------------------------------------------------

.. function:: void acb_hypgeom_bessel_j_deriv_bound(mag_t res, const acb_t nu, const acb_t z, ulong d)

Sets *res* to a bound, possibly crude, for `|J^{(d)}_{\nu}(z)|`.
Currently only specialized for small integer `\nu` and small `d`.

.. function:: void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec)

Computes the Bessel function of the first kind
Expand Down
30 changes: 29 additions & 1 deletion examples/integrals.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -743,11 +743,33 @@ f_si(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
return 0;
}

/* f(z) = J_0(z) J_1(z) J_2(z) */
int
f_triple_bessel_j(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
{
acb_t t;

if (order > 1)
flint_abort(); /* Would be needed for Taylor method. */

acb_init(t);
acb_hypgeom_bessel_j(res, t, z, prec);
acb_set_ui(t, 1);
acb_hypgeom_bessel_j(t, t, z, prec);
acb_mul(res, res, t, prec);
acb_set_ui(t, 2);
acb_hypgeom_bessel_j(t, t, z, prec);
acb_mul(res, res, t, prec);
acb_clear(t);

return 0;
}

/* ------------------------------------------------------------------------- */
/* Main test program */
/* ------------------------------------------------------------------------- */

#define NUM_INTEGRALS 38
#define NUM_INTEGRALS 39

const char * descr[NUM_INTEGRALS] =
{
Expand Down Expand Up @@ -789,6 +811,7 @@ const char * descr[NUM_INTEGRALS] =
"int_{-1-i}^{-1+i} 1/sqrt(x) dx",
"int_0^{inf} 1/gamma(x) dx (using domain truncation)",
"int_0^{1000} Si(x) dx",
"int_0^{1000} J_0(x) J_1(x) J_2(x) dx",
};

int main(int argc, char *argv[])
Expand Down Expand Up @@ -1291,6 +1314,11 @@ int main(int argc, char *argv[])
acb_calc_integrate(s, f_si, NULL, a, b, goal, tol, options, prec);
break;

case 38:
acb_zero(a);
acb_set_ui(b, 1000);
acb_calc_integrate(s, f_triple_bessel_j, NULL, a, b, goal, tol, options, prec);
break;

default:
abort();
Expand Down
1 change: 1 addition & 0 deletions src/acb_hypgeom.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,7 @@ void acb_hypgeom_m_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, i
void acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
void acb_hypgeom_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);

void acb_hypgeom_bessel_j_deriv_bound(mag_t res, const acb_t nu, const acb_t z, ulong d);
void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec);
void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec);
void acb_hypgeom_bessel_j(acb_t res, const acb_t nu, const acb_t z, slong prec);
Expand Down
Loading