%% Gauss-Laguerre Quadrature Evaluation Points and Weights
% This example shows how to solve polynomial equations and systems of equations, and work with the results using Symbolic Math Toolbox™.
%
% Gaussian quadrature rules approximate an integral by sums \int_a^b f(t) w(t) dt \approx \sum_{i=1}^n f(x_i) \alpha_i. Here, the x_i and \alpha_i are parameters of the method, depending on n but not on f. They follow from the choice of the weight function w(t), as follows. Associated to the weight function is a family of orthogonal polynomials. The polynomials' roots are the evaluation points x_i. Finally, the weights \alpha_i are determined by the condition that the method be correct for polynomials of small degree. Consider the weight function w(t) = \exp(-t) on the interval [0, \infty]. This case is known as Gauss-Laguerre quadrature.

syms t
n = 4;
w(t) = exp(-t);

%
% Assume you know the first n members of the family of orthogonal polynomials. In case of the quadrature rule considered here, they turn out to be the Laguerre polynomials.

F = laguerreL(0:n-1, t)

%
% Let L be the n+1 st polynomial, the coefficients of which are still to be determined.

X = sym('X', [1, n+1])
L = poly2sym(X, t)

%
% Represent the orthogonality relations between the Laguerre polynomials F and L in a system of equations sys.

sys = [int(F.*L.*w(t), t, 0, inf) == 0]

%
% Add the condition that the polynomial have norm 1.

sys = [sys, int(L^2.*w(t), 0, inf) == 1]

%
% Solve for the coefficients of L.

S = solve(sys, X)

%
% solve returns the two solutions in a structure array. Display the solutions.

structfun(@display, S)

%
% Make the solution unique by imposing an extra condition that the first coefficient be positive:

sys = [sys, X(1)>0];
S = solve(sys, X)

%
% Substitute the solution into L.

L = subs(L, S)

%
% As expected, this polynomial is the |n|th Laguerre polynomial:

laguerreL(n, t)

%
% The evaluation points x_i are the roots of the polynomial L. Solve L for the evaluation points. The roots are expressed in terms of the root function.

x = solve(L)

%
% The form of the solutions might suggest that nothing has been achieved, but various operations are available on them. Compute floating-point approximations using vpa:

vpa(x)

%
% Some spurious imaginary parts might occur. Prove symbolically that the roots are real numbers:

isAlways(in(x, 'real'))

%
% For polynomials of degree less than or equal to 4, you can use MaxDegree to obtain the solutions in terms of nested radicals instead in terms of root. However, subsequent operations on results of this form would be slow.

xradical = solve(L, 'MaxDegree', 4)

%
% The weights \alpha_i are given by the condition that for polynomials of degree less than n, the quadrature rule must produce exact results. It is sufficient if this holds for a basis of the vector space of these polynomials. This condition results in a system of four equations in four variables.

y = sym('y', [n, 1]);
sys = sym(zeros(n));
for k=0:n-1 
    sys(k+1) = sum(y.*(x.^k)) == int(t^k * w(t), t, 0, inf); 
end
sys

%
% Solve the system both numerically and symbolically. The solution is the desired vector of weights \alpha_i.

[a1, a2, a3, a4] = vpasolve(sys, y)
[alpha1, alpha2, alpha3, alpha4] = solve(sys, y)

%
% Alternatively, you can also obtain the solution as a structure by giving only one output argument.

S = solve(sys, y)
structfun(@double, S)

%
% Convert the structure S to a symbolic array:

Scell = struct2cell(S);
alpha = transpose([Scell{:}])

%
% The symbolic solution looks complicated. Simplify it, and convert it into a floating point vector:

alpha = simplify(alpha)
vpa(alpha)

%
% Increase the readability by replacing the occurrences of the roots x in alpha by abbreviations:

subs(alpha, x, sym('R', [4, 1]))

%
% Sum the weights to show that their sum equals 1:

simplify(sum(alpha))

%
% A different method to obtain the weights of a quadrature rule is to compute them using the formula \alpha_i = \int_a^b w(t) \prod_{j \neq i} \frac{t - x_j}{x_i - x_j} dt. Do this for i=1. It leads to the same result as the other method:

int(w(t) * prod((t - x(2:4)) ./ (x(1) - x(2:4))), t, 0, inf)

%
% The quadrature rule produces exact results even for all polynomials of degree less than or equal to 2n-1, but not for t^{2n}.

simplify(sum(alpha.*(x.^(2*n-1))) -int(t^(2*n-1)*w(t), t, 0, inf))
simplify(sum(alpha.*(x.^(2*n))) -int(t^(2*n)*w(t), t, 0, inf))

%
% Apply the quadrature rule to the cosine, and compare with the exact result:

vpa(sum(alpha.*(cos(x))))
int(cos(t)*w(t), t, 0, inf)

%
% For powers of the cosine, the error oscillates between odd and even powers:

errors = zeros(1, 20);
for k=1:20 
    errors(k) = double(sum(alpha.*(cos(x).^k)) -int(cos(t)^k*w(t), t, 0, inf));
end
plot(real(errors))

%
% Copyright 2015-2016 The MathWorks, Inc.