Re: [Help-gsl] Numerical Integration

[Sam Mason]

Hi John,

Hope you don't mind my naivety, but why are you using GSL for
integration?  Won't you get the right answer with a very simple "for"
loop?  i.e. something like:

  double sum = 0;
  for (int i = 1; i < ny; i++) {
    sum += (y[i-1] + y[i]) / 2;
  }
  return sum * deltat;

It's exact, fast and probably similar code complexity as your definition of F.

I've never used or played with the GSL integration routines, but it
wouldn't surprise me if it was getting confused by the non-smooth
behaviour at the times associated with points.  It may help debug
what's going if you print out values of x that F is getting asked for.
 Cubic spline interpolation may work better, but could introduce
artefacts of its own.

Then again, I may be barking up the wrong tree--I'm not much of a mathematician.


  Sam

On 13 April 2012 00:20, John Chludzinski <address@hidden> wrote:
> I've been trying to use gsl_integration_qags(&F, 0, 20, 0, 1e-2, 1000,
> w, &result, &error) to integrate a function defined by a large set of
> sample values taken at a constant rate.
>
> I provided a routine for F which interpolates (linearly) between
> sample values.  I've tried to play with epsabs and epsrel but get
> either:
>
> "cannot reach tolerance because of roundoff error"
>
> or
>
> "integral is divergent, or slowly convergent".
>
> BTW, simply using Simpson's rule I get: 0.148525.
> The sample values range from 3.5320 and -1.0966 with delta-t being a
> constant 0.003125.
>
> Any ideas?
>
> ---John
>