Is it possible to return an error in bessel.lisp
when, for example, the numerical value of
a bessel function is so large that the code/system
can't cope, and the answer returned for a calculation
is completely wrong?
The example below deals only with integrate
with bessel_j, but the same problem can occur
with other bessel functions.
Here we compute the integral of bessel_j(1,%i*x)
over larger and larger ranges, until we get an
obscure error message and a completely wrong
answer.
------------------------------
(%i29) float(integrate(bessel_j(1,%i*x),x,1,10));
(%o29) 2814.450562588503*%i
(%i30) float(integrate(bessel_j(1,%i*x),x,1,100));
(%o30) 1.0737517071310738E+42*%i
(%i31) float(integrate(bessel_j(1,%i*x),x,1,500));
(%o31) 2.5048094765700785E+215*%i
(%i32) float(integrate(bessel_j(1,%i*x),x,1,700));
(%o32) 1.5295933476718737E+302*%i
(%i33) float(integrate(bessel_j(1,%i*x),x,1,800));
zbesj ierr = 2
(%o33) -1.266065877752008*%i
(%i34) errcatch (float(integrate(bessel_j(1,%i*x),x,1,800)));
zbesj ierr = 2
(%o34) [-1.266065877752008*%i]
---------------------------------------
The actual answer is about %i*3.8e345, whereas
the returned answer is the wrong sign and magnitude.
Ted Woollett