Gaussian quadrature with Chebyshev polynomials

We are going to write a code to perform integration using Chebyshev polynomials for Gaussian quadrature. The idea is to approximate the integral of a function over the interval using Chebyshev nodes and weights.

The Chebyshev quadrature rule gives the nodes and weights for integrating with a weight function

For an -point Chebyshev quadrature, the nodes and weights are given by:

Here a to do list that can help:

  1. Prompt the user to choose a function to integrate from the following list:

$$\mbox{[1]}\quad f(x) = e^{-\cos(x)^2}$$

$$\mbox{[2]}\quad f(x) = x^2 $$

$$\mbox{[3]}\quad f(x) = 3x^6 + 2x^3 - 3x^2 + 2 $$

$$\mbox{[4]}\quad f(x) = 3x^{10} + 2x^3 - 3x^2 + 2 $$

$$\mbox{[5]}\quad f(x) = 3x^{100} - 3x^{50} + 2x^3 - 3x^2 + 2 $$

$$\mbox{[6]}\quad f(x) = e^{-x^2} $$

  1. Allow the user to specify the maximum number of Chebyshev nodes, , for the quadrature.
  2. Generate Chebyshev nodes and weights for the integration.
  3. Calculate and print the integral for different values of .
C code
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double f(double x, int lab) {
    switch (lab) {
        case 1: return exp(-cos(x) * cos(x));
        case 2: return x * x;
        case 3: return 3 * pow(x, 6) + 2 * pow(x, 3) - 3 * x * x + 2;
        case 4: return 3 * pow(x, 10) + 2 * pow(x, 3) - 3 * x * x + 2;
        case 5: return 3 * pow(x, 100) - 3 * pow(x, 50) + 2 * pow(x, 3) - 3 * x * x + 2;
        case 6: return exp(-x * x);
        default: return 0.0;
    }
}

void generate_weights_abscissae(double *w, double *x, int N) {
    const double pi = 3.141592653589793;
    *w = pi / N; // Chebyshev weight
    for (int j = 0; j < N; j++) {
        x[j] = cos(pi * (j + 0.5) / N); // Chebyshev node
    }
}

double gauss_chebyshev(double w, double *x, int N, int lab) {
    double integral = 0.0;
    for (int j = 0; j < N; j++) {
        integral += w * f(x[j], lab)*sqrt(1.-x[j]*x[j]);
    }
    return integral;
}

int main() {
    int Nmax, lab;

    printf("Choose the function to be integrated from -1 to 1 (just type the number):\n");
    printf("(1) exp(-cos(x)*cos(x))\n");
    printf("(2) x*x\n");
    printf("(3) 3*x^6 + 2*x^3 - 3*x^2 + 2\n");
    printf("(4) 3*x^10 + 2*x^3 - 3*x^2 + 2\n");
    printf("(5) 3*x^100 - 3*x^50 + 2*x^3 - 3*x^2 + 2\n");
    printf("(6) exp(-x*x)\n");
    scanf("%d", &lab);

    if (lab < 1 || lab > 6) {
        fprintf(stderr, "ERROR: you must insert an integer number between 1 and 6.\n");
        return 1;
    }

    printf("Enter the maximum number of Chebyshev points (Nmax): ");
    scanf("%d", &Nmax);

    if (Nmax < 2) {
        fprintf(stderr, "ERROR: Nmax must be at least 2.\n");
        return 1;
    }

    // Open the output file
    FILE *file = fopen("results.txt", "w");
    if (file == NULL) {
        fprintf(stderr, "ERROR: Could not open file for writing.\n");
        return 1;
    }

    // Write header to the output file
    fprintf(file, "#N       Integral\n");

    for (int N = 2; N <= Nmax; N++) {
        // Allocate memory for the Chebyshev nodes
        double *x = (double *)malloc(N * sizeof(double));
        if (x == NULL) {
            fprintf(stderr, "ERROR: Memory allocation failed.\n");
            fclose(file);
            return 1;
        }

        double w;
        generate_weights_abscissae(&w, x, N);

        double integral = gauss_chebyshev(w, x, N, lab);

        // Write the result to the file
        fprintf(file, "%-3d     %.6f\n", N, integral);

        // Free allocated memory
        free(x);
    }

    // Close the output file
    fclose(file);

    printf("Results written to results.txt\n");

    return 0;
}
FORTRAN code
program chebyshev_integration
  implicit none
  integer :: N,Nmax ! Number of Chebyshev points
  real(8) :: a, b, integral, W
  integer :: j, lab
  real(8), allocatable, dimension(:) :: x  ! Arrays to store nodes and weights
  
  ! Define the limits of integration
  a = -1.0d0
  b = 1.0d0

  ! Ask for input
  print *, "Choose the function to be integrated from -1 to 1 (just type the number):"
  print *, "(1) exp(-cos(x)*cos(x))"
  print *, "(2) x*x"
  print *, "(3) 3*x**6+2*x**3-3*x**2+2"
  print *, "(4) 3*x**10+2*x**3-3*x**2+2"
  print *, "(5) 3*x**100-3*x**50+2*x**3-3*x**2+2"
  print *, "(6) exp(-x*x)"

  !Read input
  read(*, *) lab

  !Check input
  if(lab.lt.0.or.lab.gt.6)then
    print*, "ERROR: you must insert an integer number between 1 and 6."
    stop
  endif
  
  print *, "Choose the function to be integrated from -1 to 1 (just type the number):"

  !Read input
  read(*, *) Nmax

  !Check input
  if(Nmax.lt.0)then
    print*, "ERROR: you must insert a positive integer number."
    stop
  endif
  
  !Header for the output
  write(666,'(A3, 7x, A3)') "#N", "Integral"

  do N = 2, Nmax

    !Allocate the array for the abscissae
    if(allocated(x))deallocate(x)
    allocate(x(N))

    ! Generate Chebyshev nodes and weights
    call generate_weights_abscissae(w,x)

    ! Perform the integration
    integral = gauss_chebyshev(w,x,lab)

    ! Print the result
    write(666,'(I3,5x,F10.6)') N, integral
  enddo
contains

  ! Define the function to integrate here
  real(8) function f(x,lab)
    real(8), intent(in) :: x
    integer,intent(in)  :: lab
    select case(lab)
    case(1)
      f = exp(-cos(x)*cos(x))                 ! Example: Integrate exp(-x^2) over [-1, 1]
    case(2)
      f = x*x
    case(3)
      f = 3*x**6+2*x**3-3*x**2+2
    case(4)
      f = 3*x**10+2*x**3-3*x**2+2
    case(5)
      f = 3*x**100-3*x**50+2*x**3-3*x**2+2
    case(6)
      f = exp(-x*x)
    end select
  end function

  subroutine generate_weights_abscissae(w,x)
    implicit none
    real(8), intent(inout)               :: w
    real(8), dimension(:), intent(inout) :: x
    real(8), parameter                   :: pi = 3.141592653589793d0
    integer :: N
    N = size(x)
    w = pi / real(N) ! Chebyshev weight
    do j = 1, N
      x(j) = cos(pi * (real(j) - 0.5d0) / real(N))  ! Chebyshev node
    end do
  end subroutine
  
  real(8) function gauss_chebyshev(w,x,lab)
    real(8), intent(in)               :: w
    real(8), dimension(:), intent(in) :: x
    integer, intent(in)               :: lab
    integer                           :: j, N
    N = size(x)
    gauss_chebyshev = 0.0d0
    do j = 1, N
      gauss_chebyshev = gauss_chebyshev + W * f(x(j),lab)*sqrt(1-x(j)*x(j))
    end do
  end function


end program chebyshev_integration