Re: [Help-gsl] Integration advice

[Will Farr]

Jonny,

You should look up some stuff on Gaussian quadrature (for example,  
it's talked about in Section 4.5 of my copy of Numerical Recipes).   
Gaussian quadrature tries to efficiently evaluate integrals of the from

\int_a^b W(x) f(x) dx,

where W(x) is chosen once for the method, and then the user supplies  
f.  It does this by choosing a set of abscissas (x_i), and weights  
(w_i) such that

\int_a^b W(x) f(x) \approx \sum_i w_i f(x_i)

A Gaussian quadrature routine of order n (with n abscissas and  
weights) will be exact for polynomial f up to order 2n-1.

If you can integrate *either* f or g against polynomials analytically  
then you can use f or g as your W function---that is, if integrals of  
the form

\int_0^1 f(r,z,t) p(t) dt    OR \int_0^1 g(t) p(t) dt,

where p(t) is a polynomial in t are tractable, then you can get a  
great improvement in accuracy by using Gaussian integration (over  
standard, general quadrature).

The GSL already includes routines for integrals of the form \int_a^b  
W(x) f(x) dx, where W is a specified function, and you supply f.  The  
W's which are supported seem to be (maybe somebody more knowledgeable  
will correct me if I get it wrong):

- W(x) = 1 (this is just usual quadrature---perfectly general, but  
wasteful if you can do better by choosing a different W)
- W(x) = (x-a)^alpha (b-x)^beta log^mu(x-a) log^nu(b-x), where a, b  
are limits of integration, and alpha, beta, mu, and nu are parameters  
you feed the routine (restricted to alpha,beta > -1, and mu, nu = {0,  
1}).
- W(x) = sin(x) or cos(x).

Your f and g do not fit these forms, but a textbook on quadrature  
(i.e. Numerical Recipes, Ch 4) should tell you how to develop your own  
sets of weights and abscissas.

Good luck!

Will
On Jun 12, 2008, at 9:24 AM, Jonny Taylor wrote:

> Hi all,
>
> I'm hoping somebody may be able to advise me on a numerical  
> integration problem I have. This may be as much a maths question as  
> a GSL question, but I'd appreciate any input people can offer.
>
> I have an integral over a variable t, an integral involving two  
> other constants r and z. The integral will be evaluated many, many  
> times for different r, z, but over the same range 0 to 1 every time.  
> i.e. the integral has the form:
>
> Int_0^1 { f(r, z, t) g(t) dt }
>
> f(r, z, t) is a relatively simple (though highly oscillatory)  
> function involving an exponential and a Bessel function. g(t) is a  
> nasty mix of trig and Legendre functions.
>
> Assuming the integral can't be solved analytically (I'm pretty sure  
> _I_ can't solve it, at least!), I'm trying to work out the fastest  
> way of numerically integrating the function. What I am hoping is  
> that there is some clever trick (hopefully implemented in gsl...)  
> that will allow me to do the hard work once at the start, and then  
> be able to repeatedly evaluate the integral for different r, z  
> relatively quickly. I feel that because the integral is "separable"  
> into f(r, z, t) and g(t) then there may be some way of speeding up  
> the computation. I have no idea what that method might be though.
>
> Any advice would be greatly appreciated...
>
> Jonny
>
>
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl