Checking if gsl_integration_qng() will work

At the end of part 1, I suggested that gsl_integraton_qng() would not be able to successfully integrate the function we want to find the Fourier series of. Here is the function I will find the Fourier series of.

Because this function is discontinuous, we have to approximate it with a Fourier series instead of a Taylor series. If we recall the formulas shown in the last part, we know that we will have to integrate \( f(x) \) from \( 0 \) to \(2l\), among other things. Before calculating the Fourier series, let's try using gsl_integration_qng() to perform basic integration.

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

  /* f(x) = x for 0 <= x < 4, repeating */
  double f(double x, void *params) {
    const double period = 4;
    return fmod(x, period);
  }

  int main() {
    double result, error;
    double low = 0, high = 8;
    gsl_function F = {.function = &f};
    double err_abs = 0, err_rel = 1;
    size_t num_evals;

    gsl_integration_qng(&F, low, high,
            err_abs, err_rel,
            &result, &error, &num_evals);
    printf("result error num_evals\n");
    printf("%f %f %zu\n", result, error, num_evals);
  }

|    result |    error | num_evals |
| 14.804436 | 8.015135 |        21 |

Well, something happened here, but it's not quite what we wanted. We told GSL to take \( \int_0^8 f(x) dx \). Since we're dealing with a discontinuous function, we can rewrite this as \( \int_0^4 x dx + \int_4^8 x dx \) Now I'm no mathematician, but this answer should be 16, not 14.8. Additionally, our error is very large, and if I attempt to lower err_rel to a smaller value, GSL yells at me saying that it failed to reach the requested error.

This is not a problem on GSL's end. We aren't using the right integration function for the job! According to the GSL manual,

The QNG algorithm is a non-adaptive procedure which uses fixed Gauss-Kronrod-Patterson abscissae to sample the integrand at a maximum of 87 points. It is provided for fast integration of smooth functions.

While that's a lot of words, the important part is the last sentence. gsl_integration_qng() is for smooth functions. Our function is not smooth, it has a discontinuity! We will need to find an alternative function in GSL that can handle discontinuous graphs.[fn:1:Strictly speaking, this isn't actually true. Because the discontinuity occurs at the limits of integration for our function (as long as we only integrate from 0 to 4), we can use gsl_integration_qng(). Not all functions are like this, so it's still best to find an alternative.]

gsl_integration_qag() and discontinuities

Fortunately we do not have to look far. The very next function mentioned in the manual is what we need. gsl_integration_qag() has the following signature.

  int gsl_integration_qag(const gsl_function * f,
              double a, double b,
              double epsabs, double epsrel,
              size_t limit, int key,
              gsl_integration_workspace * workspace,
              double * result, double * abserr)

f, a, b, epsabs, epsrel, result, and abserr are all the same as gsl_integration_qng(). We can see that several new terms are introduced however.

workspace is a pointer to an area of memory used by gsl_integration_qag(). We allocate space by using the gsl_integration_workspace_alloc() function, This function is passed an integer to adjust how much memory we want to allocate. Whatever integer we pass to gsl_integration_workspace_alloc(), we need to ensure that limit is the same. This way, gsl_integration_qag() knows how much memory it has available.

key is a integer between 0 through 6. GSL recommends using higher values when integrating smooth functions, and lower values when functions are discontinuous.

Let's see how well gsl_integration_qag() performs!

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

  double f(double x, void *params) {
    const double period = 4;
    return fmod(x, period);
  }

  int main() {
    size_t limit = 1024;
    gsl_integration_workspace *w
      = gsl_integration_workspace_alloc(limit);
    double result, error;
    double low = 0, high = 8;
    gsl_function F = {.function = &f};
    double err_abs = 0, err_rel = 1e-7;

    gsl_integration_qag(&F, low, high,
            err_abs, err_rel,
            limit, 1, w,
            &result, &error);
    /* get into the habit of freeing memory when done! */
    gsl_integration_workspace_free(w);
    printf("result error\n");
    printf("%f %f\n", result, error);
  }

| result | error |
|   16.0 |   0.0 |

Look at that! We are getting the exact value we expected, even though we're integrating with a discontinuity. We are now at the point where we can calculate the Fourier series for our function.

Generating our Fourier series

Let's recall that we can find our \(a_n\) and \(b_n\) coefficients with the following formulas.

\[ a_n = \frac{1}{l} \int_{0}^{2l} f(x) \cos(\frac{n \pi x }{l}) dx \]

\[ b_n = \frac{1}{l} \int_{0}^{2l} f(x) \sin(\frac{n \pi x }{l}) dx \]

GSL expects us to only pass one function to it as an argument. That means that, unlike before, we can't just pass a pointer to our function f(x), as GSL won't be multiplying it by the cosine and sine terms.

There are a couple of solutions to this that I can think of. The first is to make use of the void * argument that GSL includes in the gsl_function structure. Using this, we could pass extra information to f(x), adapting it to our specific n value.

Alternatively, we could write a parent function, like aorb_subn(). This function would then be stored in the gsl_function structure. The advantage of this approach is that the function we want the Fourier series of, f(x), is distinct in our code. This should make it a bit easier to change f(x) to any function that we want. This is the approach that I will take.

To make this code hopefully a bit more modular, I moved the integration out of main, and instead into a function called get_aorb_subn(). This function calculates the nth Fourier coefficient for a function f(x), using the variable get_a to determine if it should calculate \( a_n \) or \( b_n \). It may be better to wrap this in a fourier structure that contains the two terms, but I elected not to do that to hopefully maintain some vague semblance of readability.

This is what the code to calculate the Fourier series looks like.

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

  #define PI 3.14159
  const double period = 4;
  const double l = period / 2;

  struct aorb_params {
    double (*f)(double x, void *params);
    int n;
    bool calc_a;
  };

  double f(double x, void *params) {
    return fmod(x, period);
  }

  double aorb_subn(double x, struct aorb_params *params) {
    int n = params->n;
    double trig = params->calc_a ? cos(PI*n*x/l) : sin(PI*n*x/l);
    return params->f(x, NULL) * trig;
  }

  double get_aorb_subn(double (*f)(double x, void *params), int n, 
               double low, double high, bool get_a) {
    double normalization = 1/l;
    if (get_a && n == 0) normalization /= 2;

    size_t limit = 1024;
    gsl_integration_workspace *w
      = gsl_integration_workspace_alloc(limit);
    double result, error;
    gsl_function F = { .function = &aorb_subn,
      .params = &(struct aorb_params){ .f = f, .n = n,
             .calc_a = get_a } };
    double err_abs = 0, err_rel = 1e-7;

    gsl_integration_qag(&F, low, high,
            err_abs, err_rel,
            limit, 1, w,
            &result, &error);
    gsl_integration_workspace_free(w);

    return normalization * result;
  }

  int main() {
    const int upto = 10;
    double a_subn[upto], b_subn[upto];
    double low = 0, high = period;

    for (int i = 0; i < upto; i++) {
      a_subn[i] = get_aorb_subn(f, i, low, high, true);
      b_subn[i] = get_aorb_subn(f, i, low, high, false);
      printf("%d %f %f\n", i, a_subn[i], b_subn[i]);
    }
  }

| n | a_subn |    b_subn |
|---+--------+-----------|
| 0 |    2.0 |       0.0 |
| 1 | -7e-06 | -1.273242 |
| 2 | -7e-06 | -0.636621 |
| 3 | -7e-06 | -0.424414 |
| 4 | -7e-06 |  -0.31831 |
| 5 | -7e-06 | -0.254648 |
| 6 | -7e-06 | -0.212207 |
| 7 | -7e-06 | -0.181892 |
| 8 | -7e-06 | -0.159155 |
| 9 | -7e-06 | -0.141471 |

Like I mentioned before, aorb_subn is a "parent function" that combines f(x) and the sine/cosine term. This is the function we integrate, and we provide it enough information to know what term we are calculating.

Interestingly, if we use the \( \pi \) constant M_PI, the integration fails due to roundoff error. Fortunately we can get "close enough" values by just using a few less digits.

Once we graph the Fourier series generated by these coefficients, this is what we get.

And there we are! Here is a Fourier series generated for a discontinuous periodic function. The more terms we add, the more accurate the approximation.