quadratura-gaussiana.c

/* Importando bibliotecas necessárias */
#include <stdio.h> 
#include <math.h>

/* Definindo uma constante N que vai ser o valor para o polinômio de Legendre LN*/
#define N 3

double Pi;
double lroots[N]; /* Array com raízes do polinômio de Legendre */
double weight[N]; /* Pesos da Integração */

/* Função que retorna o valor de um polinômio de Legendre aplicado sobre um ponto 
   o parâmetro n representa a ordem do polinômio de Legendre e x a icógnita
   
   Lembrando da equação de recorrência de um polinômio de Legendre, que é uma fórmula recursiva
   Ln(x) = (2*n-1)*x*Ln-1(x)/n -(n-1)*Ln-2(x)/n
   
   Sendo que L0(x)=1 e L1(x)=x
*/
double lege_eval(int n, double x)
{
        if(n == 1) {
          //Condição que representa que L1(x)=x
          return x;
        } else if(n == 0) {
          //Condição que representa que L0(x)=1
          return 1;
        }
        
        //Equação recursiva de um polinômio de Legendre de ordem n>2
        double retorno = ((2*n-1)*x*lege_eval(n-1, x) -(n-1)*lege_eval(n-2, x))/n;
  
        return retorno;
}

/* Função para encontrar derivada do polinômio de Legendre
   de ordem n no ponto x
*/
double lege_diff(int n, double x)
{
        /* retorna a derivada de um polinômio de Legendre como
           n(x*Ln(x)-Ln-1(x))/(x^2-1)
        */
        return n * (x * lege_eval(n, x) - lege_eval(n - 1, x)) / (x * x - 1);
}

/* Função para encontrar raízes do polinômio de Legendre */
void lege_roots()
{
        int i;
        double x, x1;
  
        /* Acha as n raízes do polinômio de Legendre*/
        for (i = 1; i <= N; i++) {  
                /*Inicializa um valor de x*/
                x = cos(Pi * (i - .25) / (N + .5));
                
                /* Encontra as raízes através do método de Newton-Rampson
                   que consiste em xk = xk-1 - f(xk-1)/f'(xk-1)
                */
                do {
                       x1 = x;
                       x -= lege_eval(N, x) / lege_diff(N, x);
                } while ( fdim( x, x1) > 2e-16 );

               lroots[i - 1] = x;
 
               /* Adiciona o valor da derivada de Ln(x) à derivada x1*/
               x1 = lege_diff(N, x);
          
               /* Cria o coeficiente de peso Ai-1 sendo que a fórmula geral para os coeficientes é
                  Ai = 2/((1-x^2)*x1^2)
                  onde x1 é a derivada da função no ponto x, que é uma raiz do polinômio de Legendre
               */
               weight[i - 1] = 2 / ((1 - x * x) * x1 * x1);
        }
}
 
/* Função que retorna o valor da integral da função, recebe como parâmetro 
   o ponteiro para uma função e os extremos do intervalo de integração definida
*/
double lege_inte(double (*f)(double), double a, double b)
{
        /* Transformação de coordenadas de (a,b) para (-1, 1) utilizando dois coeficientes auxiliares*/
        double c1 = (b - a) / 2, c2 = (b + a) / 2, sum = 0;
  
        int i;
  
        /* Laço que soma para o valor final da integral a expressao Ai*f(c1x + c2) 
           onde c1 e c2 são coeficientes para transformar o intervalo em (-1, 1)
        */
        for (i = 0; i < N; i++)
               sum += weight[i] * f(c1 * lroots[i] + c2);
        
        /* O valor final retornado da integral*/
        return c1 * sum;
}
 
int main()
{
        int i;
        
        /* Declara o valor como arctg(1/1) * 4 */
        Pi = atan2(1, 1) * 4;
 
        /* Chama as funções de carregar a matriz com as raízes do polinômio de Legendre*/
        lege_roots();
 
        /* Imprime as raízes */
        printf("Roots: ");
        for (i = 0; i < N; i++)
                printf(" %g", lroots[i]);
 
        /* Imprime os valores dos pesos */
        printf("\nWeight:");
        for (i = 0; i < N; i++)
                printf(" %g", weight[i]);
 
        /* Imprime o valor da integral e compara com o valor real */
        printf("\nintegrating Exp(x) over [-3, 3]:\n\t%10.8f,\n"
                "compred to actual\n\t%10.8f\n",
               lege_inte(exp, -3, 3), exp(3) - exp(-3));
        return 0;
}