GSL integrals, bad counts

Question

Why is this program wrong; count integral tan(x) in the range (0, pi/2] (calculated around ~39) where Wolfram Alpha says it's ~7.

my code:

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f (double x, void * params) {
  double alpha = *(double *) params;
  double f = tan(x);
  return f;
}

int
main (void)
{
  gsl_integration_workspace * w 
    = gsl_integration_workspace_alloc (1000);

  double result, error;
  double expected = -4.0;
  double alpha = 1.0;
  gsl_function F;
  F.function = &f;
  F.params = &alpha;
gsl_set_error_handler_off();

  gsl_integration_qag (&F, 0, M_PI/2, 0, 1e-6, 1000, 1, 
                        w, &result, &error); 

  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("estimated error = % .18f\n", error);
  printf ("actual error    = % .18f\n", result - expected);
  printf ("intervals =  %d\n", w->size);

  gsl_integration_workspace_free (w);

  return 0;
}

if i deleted gsl_set_error_handler_off(); i had error "bad integrand behavior".

Answer 1

That integral does not have a finite value.

tan(x) = sin(x)/cos(x)

so

tan(pi/2) = 1/0 = undefined.

Thus, your numerical integral ought to diverge because your function is infinite at the edge of its range.

You can see this analytically:

∫tan(x)dx
= ∫sin(x)/cos(x)dx
Define u=cos(x).  Then du=-sin(x)dx, so
∫sin(x)/cos(x)dx
=-∫(-sin(x))/cos(x)dx
=-∫1/u*du
=-ln|u| + C 
=-ln|cos(x)| + C

So, the integral from 0 to pi/2

= -ln|cos(pi/2)| - (-ln|cos(0)|)
= -ln|0| + 0

But, -ln(0) is undefined, and approaches positive infinity as x -> +0. Numerical integration algorithms will attempt to approximate this infinite integral by summing areas under known slices, and produce wrong, large finite results with error detection disabled. With error detection enabled a good numerical integration algorithm will correctly report an error such as a failure to converge or error evaluating the integrand - which is exactly what you are seeing with gsl when you enable its error detection.

Wolfram Alpha also reports an infinite value for ∫tan(x)dx at pi/2, so I'm not sure where you got the value of ~7.