4
$\begingroup$

I’m trying to evaluate an improper integral of a 0th order Bessel function of the first kind using Matlab:

v = integral(@(x)besselj(0, x), 0, Inf)

which returns v = 3.7573e+09. However, this should be 1, in theory. What should I do to get 1 in the Matlab?

Improve this question
$\endgroup$
2
  • $\begingroup$ I assume that you also receive a warning, as I do in R2015a, that the "limit on the maximum number of intervals" has been reached and that "it may be difficult to approximate numerically to the requested accuracy"? If so, edit your question, to provide the full warning message. Do you want to do this numerically, or are symbolic solutions acceptable? Since you know what the answer is, this might seem pointless, but I'm guessing it relates to this Math.StackExchange question? $\endgroup$
    – horchler
    Commented Mar 9, 2015 at 19:03
  • $\begingroup$ yes, actually this related to that question that I asked. I am the same person ;). I know the answer just for specific inputs not for every inputs. this is the analytical solution of one specific problem that I've already solve numerically and by solving this equation numerically I want to compare with my numeric method. but unfortunately I'm sure that the answer in wrong here. as I said in the question $\endgroup$
    – Rasa
    Commented Mar 10, 2015 at 20:09

1 Answer 1

Reset to default
4
$\begingroup$

The problem is that you are integrating an oscillatory function over an infinite interval. The MATLAB website doesn't give specifics on the algorithm behind their integral function (it just says 'globally adaptive'). The older quad function used to be adaptive Simpson, so I can assume that integral is the same. What probably is happening is that the code transforms the interval $[0,\infty]$ to $[0,1]$ or $[-1,1]$ or so and hence it is squeezing all the oscillations in this interval. As a result it is adding a lot of positive contributions and subtracting a lot of negative contributions. Hence a totally wrong answer.

For oscillatory integrals like these a couple of tricks exist. First of all, if the function is rapidly decaying, you could truncate the interval (and as such the classical routines will not transform the interval). If you have ideas on the zero crossings, you could loop over all these and integrate between the zeros, taking care to sum positive parts and negative parts first and then subtract in the end.

A second approach might be the double exponential quadrature formula adapted for oscillatory integrals. Please refer to Ooura and Mori, "The double exponential formula for oscillatory functions over the half infinite interval", J. Comp. App. Math. , 38, p353-360, 1991. Ooura has FORTRAN and C implementations on his website.

Finally, since you already know the correct answer to your integral, why do it numerically? Just to understand why integral fails or is there a bigger picture?

Improve this answer
$\endgroup$
1
  • $\begingroup$ thank you so much for your answer, for your question I have to solve this (( math.stackexchange.com/questions/1165170/… )) equation for a list of data and for simplicity I asked this question here. Also I know the answer just for a specific inputs which are boundary or initial conditions that I used to solve the equation not for every input. $\endgroup$
    – Rasa
    Commented Mar 10, 2015 at 20:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.