Question

Write a program to compute 2 point Gauss quadrature on a single interval. Your input should be a function with f,a,and b. Verify it by checking that the error in using the rule for f(x)=x^3 - 3x^2 +x - 7 is zero. In matlab please.

Answer 1

%Matlab Code for Gaussian Quadrature integration

clear all
close all

%function for Gaussian Quadrature
f=@(x) x.^3 - 3*x.^2 +x - 7;
%for 2 point Gaussian quadrature N=2
N=2;
%upper and lower limit
a=0; b=1;

%integration using gaussian quadrature 2 point
value=gauss_quad(f,N,a,b);

%exact integration
ext_val=integral(f,a,b);

fprintf('For the function f(x)=')
disp(f)

fprintf('Integration value for a=%f to b=%f using Gaussian Quadrature is %f\n',a,b,value)
fprintf('Integration value for a=%f to b=%f using exact integration is %f\n',a,b,ext_val)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Function for Gaussian quadrature
function value=gauss_quad(f,N,a,b)

%
% This script is for computing definite integrals using Legendre-Gauss
% Quadrature. Computes the Legendre-Gauss nodes and weights on an interval
% [a,b] with truncation order N
%
% Suppose you have a continuous function f(x) which is defined on [a,b]
% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f.

    N=N-1;
    N1=N+1; N2=N+2;

    xu=linspace(-1,1,N1)';

    % Initial guess
    y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);

    % Legendre-Gauss Vandermonde Matrix
    L=zeros(N1,N2);

    % Derivative of LGVM
    Lp=zeros(N1,N2);

    % Compute the zeros of the N+1 Legendre Polynomial
    % using the recursion relation and the Newton-Raphson method

    y0=2;

    % Iterate until new points are uniformly within epsilon of old points
    while max(abs(y-y0))>eps

        L(:,1)=1;
        Lp(:,1)=0;

        L(:,2)=y;
        Lp(:,2)=1;

        for k=2:N1
            L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
        end

        Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);

        y0=y;
        y=y0-L(:,N2)./Lp;

    end

    % Linear map from[-1,1] to [a,b]
    xx=(a*(1-y)+b*(1+y))/2;    

    % Compute the weights
    w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
    value = sum(w.*f(xx));
end

%---------------End of Code---------------%