Gauss-Legendre and Gauss-Jacobi quadrature

Contents

You need to download two m-files

Download GaussLegendre.m and GaussJacobi.m

Gauss-Legendre quadrature on [-1,1]

[x,w]=GaussLegendre(n) computes quadrature nodes $x_1,\ldots,x_n$ and weights $w_1,\ldots,w_n$ such that

$$\int_{-1}^1 f(x) dx \approx w_1 f(x_1)+\cdots+ w_n f(x_n)$$

Example: Use n=2,4,6,8 nodes to approximate

$$I=\int_{-1}^1 (x+1)^{-1/2}e^x dx$$

Note that the function goes to infinity at the left endpoint -1, so it is not surprising that the errors are fairly large. They converge to zero as $n\to\infty$.

f = @(x) (x+1).^(-1/2).*exp(x); % define f(x) so that it works for vectors

syms X                          % compute "exact" value of integral with vpaintegral
exact = double(vpaintegral(f(X),-1,1,'AbsTol',1e-20,'RelTol',1e-20));

for n=2:2:8
  [x,w] = GaussLegendre(n);     % find nodes and weights of Gauss-Legendre rule with n nodes
  Q = sum(w.*f(x));             % Q = w(1)*f(x(1)) + ... + w(n)*f(x(n))
  error = Q - exact;
  fprintf('n=%2g, error=%g\n',n,error)
end

n= 2, error=-0.178422
n= 4, error=-0.0995924
n= 6, error=-0.0693315
n= 8, error=-0.0531315

Gauss-Jacobi quadrature on [-1,1]

[x,w]=GaussJacobi(n,alpha,beta) computes quadrature nodes $x_1,\ldots,x_n$ and weights $w_1,\ldots,w_n$ such that

$$\int_{-1}^1 g(x) (x+1)^\alpha (1-x)^\beta dx \approx w_1 g(x_1)+\cdots +w_n g(x_n)$$

Note that GaussJacobi(n,0,0) is the same as GaussLegendre(n).

Example: Use n=2,4,6,8 nodes to approximate

$$I=\int_{-1}^1 (x+1)^{-1/2}e^x dx$$

Since the integrand contains the singular term $(x+1)^{-1/2}$ at the left endpoint -1 we choose $\alpha=-1/2$, $\beta=0$ and $g(x)=e^x$.

Since the function g(x)=exp(x) is smooth, the errors are very small and rapidly converge to 0.

g = @(x) exp(x);                % define g(x). Note: (x+1)^(-1/2) term is part of quadrature rule
for n=2:2:8
  [x,w] = GaussJacobi(n,-.5,0);    % find nodes and weights of Gauss-Jacobi rule, alpha=-1/2, beta=0
  Q = sum(w.*g(x));             % Q = w(1)*g(x(1)) + ... + w(n)*g(x(n))
  error = Q - exact;
  fprintf('n=%2g, error=%g\n',n,error)
end

n= 2, error=-0.0102284
n= 4, error=-4.04491e-07
n= 6, error=-2.17115e-12
n= 8, error=0

Gauss-Legendre quadrature on [a,b]

We can use nodes and weights from [x,w]=gaussleg(n) to approximate an integral on an interval [a,b]: Let $r:=(b-a)/2$ and $t_j:=(a+b)/2+rx_j$, then

$$\int_{a}^b f(x) dx \approx r \Bigl[ w_1 f(t_1)+\cdots+w_n f(t_n) \Bigr]$$

Example: Use n=2,4,6,8 nodes to approximate

$$I=\int_0^1 x^{-1/3}e^x dx$$

Note that the function goes to infinity at the left endpoint 0, so it is not surprising that the errors are fairly large. They converge to zero as $n\to\infty$.

f = @(x) x.^(-1/3).*exp(x);      % define f(x) so that it works for vectors

syms X                           % compute "exact" value of integral with vpaintegral
exact = double(vpaintegral(f(X),0,1,'AbsTol',1e-20,'RelTol',1e-20));

a = 0; b = 1;
r = (b-a)/2;

for n=2:2:8
  [x,w] = GaussLegendre(n);      % find nodes and weights of Gauss-Legendre rule with n nodes
  t = (a+b)/2 + r*x;             % nodes for [a,b]
  Q = r*sum(w.*f(t));            % Q = r*( w(1)*f(t(1)) + ... + w(n)*f(t(n)) )
  error = Q - exact;
  fprintf('n=%2g, error=%g\n',n,error)
end

n= 2, error=-0.115784
n= 4, error=-0.0536563
n= 6, error=-0.0329847
n= 8, error=-0.0230987

Gauss-Jacobi quadrature on [a,b]

We can use nodes and weights from [x,w]=GaussJacobi(n,alpha,beta) to approximate an integral on an interval [a,b]: Let $r:=(b-a)/2$ and $t_j:=(a+b)/2+rx_j$, then

$$\int_{a}^b g(x) (x-a)^\alpha (b-x)^\beta dx \approx r^{\alpha+\beta+1} \Bigl[ w_1 g(t_1)+\cdots+ w_n g(t_n) \Bigr]$$

Example: Use n=2,4,6,8 nodes to approximate

$$I=\int_0^1 x^{-1/3}e^x dx$$

Since the integrand contains the singular term $x^{-1/3}$ at the left endpoint 0 we choose $\alpha=-1/3$, $\beta=0$ and $g(x)=e^x$.

Since the function g(x)=exp(x) is smooth, the errors are very small and rapidly converge to 0.

g = @(x) exp(x);                 % define g(x). Note: x^(-1/3) term is part of quadrature rule
a = 0; b = 1;
r = (b-a)/2;

for n=2:2:8
  [x,w] = GaussJacobi(n,-1/3,0); % find nodes and weights of Gauss-Jacobi rule, alpha=-1/3, beta=0
  t = (a+b)/2 + r*x;             % nodes for [a,b]
  Q = r^(1-1/3+0)*sum(w.*g(t));  % Q = r^(1+alpha+beta)*( w(1)*g(t(1)) + ... + w(n)*g(t(n) )
  error = Q - exact;
  fprintf('n=%2g, error=%g\n',n,error)
end

n= 2, error=-0.000600555
n= 4, error=-1.46867e-09
n= 6, error=-2.22045e-15
n= 8, error=-1.33227e-15


Published with MATLAB® R2018a