Numerical Integration of sin(1/x)

[Urkel]

Hi everyone,
I am writing a simple code using Numerical Recipes (that bible of numerical method) to

integrate using trapezoid rule the following integral
int_pi/2_inf {sin(x)/x^2} dx
I first make variable change y = 1/x to change limit of integration so that now the integral

becomes int_0_2/pi sin(1/y) dy
I ever took analysis and this function has very special properties because of its extremely rapid

oscillation near x=0. Also, it has a special name, I forgot what it is. Anybody knows what

function sin(1/x) is commonly called in analysis?
Because of this rapid oscillation, one must use much finer sampling density when employing

trapezid rule to get accurate answer.What do you think further simple analytic transformation

that I can make to improve the convergence of this numerical integration using trapezoid rule?

Suggestions would be highly appreciated.

Urkel

 

[rsq_a]

Hi everyone,
I am writing a simple code using Numerical Recipes (that bible of numerical method) to

integrate using trapezoid rule the following integral
int_pi/2_inf {sin(x)/x^2} dx
I first make variable change y = 1/x to change limit of integration so that now the integral

becomes int_0_2/pi sin(1/y) dy
I ever took analysis and this function has very special properties because of its extremely rapid

oscillation near x=0. Also, it has a special name, I forgot what it is. Anybody knows what

function sin(1/x) is commonly called in analysis?
Because of this rapid oscillation, one must use much finer sampling density when employing

trapezid rule to get accurate answer.What do you think further simple analytic transformation

that I can make to improve the convergence of this numerical integration using trapezoid rule?

Suggestions would be highly appreciated.

Urkel

Yes, it's related to the Cosine Integral.

It's not an easy integral to do. Essentially, your best bet is to use high-order series expansions of the solution.

Another route is to go with a spectral method like Fourier or Chebyshev functions. Just make sure that you verify convergence as the number of basis functions increases. I'm fairly sure you won't get anywhere fast with typical finite difference integration.

 

[Urkel]

Thanks for the advice but I particulalry intend to work with finite difference methods.
I use 3 methods at once; trapezoid, Romberg and Gaussian quadrature to integrate out int_pi/2_inf sin(x)/x^2 dx followed first by defining y = 1/x and transforming the integral to int_0_2/pi sin(1/y) dy.
All the methods fail in doing the integral in the latter form when I set the integration limits to something like 0.001 for lower limit and 2/pi for the upper one. When I set the lower to, say, 0.5, all three methods give exactly the same answer so at least I know my codes for 3 algorithms are correct. So, what should I do, what further analytic transformation I need to do, if any, to make the integral converges (the code works) when I set lower limit back to small value close to 0?

Regards

 

[maze]

As x gets smaller and smaller, each of the "bumps" on the curve basically looks like a little spike, with zeros at 1/(n*Pi) and 1/((n+1)*Pi), and maximum 1 at 1/((n+1/2)*Pi).

You could do an integration scheme where you approximate the function with spikes for x<a (for small 'a' determined by error estimates), and with better approximation for x>=a.

Here is an image:
http://img132.imageshack.us/img132/259/spikyintegrationpk1.png

Then you can just sum the alternating series rather than actually evaluating points.

This can be made completely rigorous in terms of getting within an error tolerance, since you also have the bound on the derivative on that interval, which is the derivative on the left hand side. For x in [1/nPi, 1/(n+1)Pi], you have |f'(x)| < |f'(1/(n+1)Pi)| = |cos(1/(n+1)Pi)/((n+1)Pi)^2|. Then you can use this to bound the error in approximation of each spike with the error estimate for polynomial interpolation.