Gaussian quadrature for arbitrary weight functions

Dear all,

there are several Matlab codes available to compute an integral of the form

with Gaussian quadrature, where is a "standard" weight function, e.g. .

I'm looking for a Gaussian quadrature that works for arbitrary weight functions. For instance, a code of the form

[Quadrature_weights,Quadrature_nodes] = CODE(weight_function)

that returns Gaussian quadrature weights and nodes for a given, but arbitrary weight function.

Thanks, Stephan

Answer 1

[Nipun]

Hi Stephan,

I understand that you intend to perform Gaussian quadrature for arbitrary weight functions in MATLAB. I recommend the following steps to achieve the desired results:

Define the weight function and desired number of nodes:

weight_function = @(x) exp(-x);  % example weight function
n = 5;  % number of nodes

Generate the quadrature nodes and weights:

You can use the Golub-Welsch algorithm to find the nodes and weights for any weight function. This method involves computing the eigenvalues and eigenvectors of a Jacobi matrix constructed from the moments of the weight function.

function [nodes, weights] = custom_gaussian_quadrature(weight_function, n)
    % Define a suitable interval for integration, e.g., [a, b]
    a = -1;  % start of interval
    b = 1;   % end of interval
    
    % Compute the moments of the weight function
    moments = zeros(1, 2*n);
    for k = 1:2*n
        moments(k) = integral(@(x) (x.^(k-1)) .* weight_function(x), a, b);
    end
    
    % Construct the Jacobi matrix
    J = zeros(n);
    for i = 1:n
        J(i,i) = moments(2*i-1) / moments(2*i-2);  % diagonal elements
        if i < n
            J(i,i+1) = sqrt(moments(2*i) * moments(2*i-1)) / moments(2*i-2);  % off-diagonal elements
            J(i+1,i) = J(i,i+1);  % symmetric matrix
        end
    end
    
    % Compute the eigenvalues and eigenvectors of the Jacobi matrix
    [V, D] = eig(J);
    nodes = diag(D);  % nodes are the eigenvalues
    weights = moments(1) * V(1,:)'.^2;  % weights are the square of the first element of eigenvectors
end

% Use the function to compute nodes and weights
[nodes, weights] = custom_gaussian_quadrature(weight_function, n);

Use the nodes and weights for integration:

integral_value = sum(weights .* f(nodes));  % where f is the function to be integrated

This approach gives you the Gaussian quadrature nodes and weights for any arbitrary weight function. Adjust the interval [a, b] as per the specific requirement.

Hope this helps.

Regards,

Nipun

[Torsten]

For i = 1,

J(i,i) = moments(2*i-1) / moments(2*i-2);  % diagonal elements

throws an error because it's attempted to access moments(0).

AI failed ?