Multiple Discretization

Multiple Discretization

This tutorial shows how to compute recurrence coefficients for non-trivial weight functions, and how they are being used for quadrature. The method we use is called multiple discretization, and follows W. Gautschi's book "Orthogonal Polynomials: Computation and Approximation", specifically Section 2.2.4, and Example 2.38.

Suppose we have the weight function

\[\forall t \in [-1,1], \gamma \in [0,1]: \quad w(t;\gamma) = \gamma + (1-\gamma) \frac{1}{\sqrt{1-t^2}},\]

and we would like to solve

\[\int_{-1}^{1} f(t) w(t;c) \mathrm{d}t = \sum_{\nu=1}^{N} f(\tau_\nu) w_\nu\]

by some quadrature rule. We will see that ad-hoc quadrature rules will fail to solve the integral even for the simplest choice $f \equiv 1$. However, finding the recurrence coefficients of the underlying orthogonal polynomials, and then finding the quadrature rule will be the way to go.

Let us first try to solve the integral for $f \equiv 1$ by Fejer's rule.

using PolyChaos, LinearAlgebra
γ = 0.5;
int_exact = 1+pi/2; # exact value of the integral
function my_w(t::Float64,γ::Float64)
    γ + (1-γ)*1/sqrt(1-t^2)
end

N = 1000;
n,w = fejer(N);
int_fejer = dot(w,my_w.(n,γ))
print("Fejer error:\t$(abs(int_exact-int_fejer))\twith $N nodes")
Fejer error:	0.00034489625618583375	with 1000 nodes

Clearly, that is not satisfying. Well, the term $\gamma$ of the weight $w$ makes us think of Gauss-Legendre integration, so let's try it instead.

function quad_gaussleg(N::Int,γ::Float64)
    a,b=rm_legendre(N)
    n,w=golubwelsch(a,b)
end
n,w = quad_gaussleg(N+1,γ)
int_gaussleg = dot(w,γ .+ (1-γ)/sqrt.(1 .- n.^2))
print("Gauss-Legendre error:\t$(abs(int_exact-int_gaussleg))\twith $N nodes")
Gauss-Legendre error:	1.5692263158654423	with 1000 nodes

Even worse! Well, we can factor out $\frac{1}{\sqrt{1-t^2}}$, making the integral amenable to a Gauss-Chebyshev rule. So, let's give it anothery try.

function quad_gausscheb(N::Int64,γ::Float64)
    a,b = rm_chebyshev1(N)
    n,w = golubwelsch(a,b)
end
n,w = quad_gausscheb(N+1,γ)
int_gausscheb = dot(w,γ .+ (1-γ)*sqrt.(1 .- n.^2))
# int=sum(xw(:,2).*(1+sqrt(1-xw(:,1).^2)))
print("Gauss-Chebyshev error:\t$(abs(int_exact-int_gausscheb))\twith $(length(n)) nodes")
Gauss-Chebyshev error:	4.112336333683686e-7	with 1000 nodes

Okay, that's better, but it took us a lot of nodes to get this result. Is there a different way? Indeed, there is. As we have noticed, the weight $w$ has a lot in common with Gauss-Legendre and Gauss-Chebyshev. We can decompose the integral as follows

\[\int_{-1}^1 f(t) w(t) \mathrm{d}t = \sum_{i=1}^{m} \int_{-1}^{1} f(t) w_i(t) \mathrm{d} t,\]

with

\[\begin{align*} w_1(t) &= \gamma \\ w_2(t) &= (1-\gamma) \frac{1}{\sqrt{1-t^2}}. \end{align*}\]

To the weight $w_1$ we can apply Gauss-Legendre quadrature, to the weight $w_2$ we can apply Gauss-Chebyshev quadrature (with tiny modifications). This discretization of the measure can be used in our favor. The function mcdiscretization() takes the $m$ discretization rules as an input

function quad_gaussleg_mod(N::Int,γ::Float64)
    n,w = quad_gaussleg(N+1,γ)
    return n,γ*w
end
function quad_gausscheb_mod(N::Int,γ::Float64)
            n,w = quad_gausscheb(N+1,γ)
    return n,(1-γ)*w
end

N = 8
a,b = mcdiscretization(N,[n->quad_gaussleg_mod(n,γ); n->quad_gausscheb_mod(n,γ)])
n,w = golubwelsch(a,b)
int_mc = sum(w)
print("Discretization error:\t$(abs(int_exact-int_mc))\twith $(length(n)) nodes")
Discretization error:	4.440892098500626e-16	with 7 nodes

Et voilà, no error with fewer nodes. (For this example, we'd need in fact just a single node.)

The function mcdiscretization() is able to construct the recurrence coefficients of the orthogonal polynomials relative to the weight $w$. Let's inspect the values of the recurrence coefficients a little more. For $\gamma = 0$, we are in the world of Chebyshev polynomials, for $\gamma = 1$, we enter the realm of Legendre polynomials. And in between? That's exactly where the weight $w$ comes in: it can be thought of as an interpolatory weight, interpolating Legendre polynomials and Chebyshev polynomials. Let's verify this by plotting the recurrence coefficients for several values of $\gamma$.

Γ = 0:0.1:1;
ab = [ mcdiscretization(N,[n->quad_gaussleg_mod(n,gam); n->quad_gausscheb_mod(n,gam)]) for gam in Γ ];
bb = hcat([ ab[i][2] for i=1:length(Γ)]...);
b_leg = rm_legendre(N)[2];
b_cheb = rm_chebyshev1(N)[2]
bb[:,1]-b_cheb
8-element Array{Float64,1}:
 -1.7763568394002505e-15
 -3.3306690738754696e-16
  3.3306690738754696e-16
 -8.326672684688674e-17 
  2.220446049250313e-16 
  3.3306690738754696e-16
 -8.326672684688674e-17 
  1.6653345369377348e-16
bb[:,end]-b_leg
8-element Array{Float64,1}:
  8.881784197001252e-16 
 -4.440892098500626e-16 
  1.1102230246251565e-16
 -4.440892098500626e-16 
  2.7755575615628914e-16
 -1.6653345369377348e-16
  4.440892098500626e-16 
 -6.661338147750939e-16 

Let's plot these values to get a better feeling.

using Plots
gr()
plot(Γ,bb',yaxis=:log10, w=3, legend=false)
zs, os = zeros(N), ones(N)
scatter!(zs,b_cheb,marker=:x)
scatter!(os,b_leg,marker=:circle)

xlabel!("gamma")
ylabel!("beta")
title!("N=$N Recurrence Coefficients Interpolating from Chebyshev to Legendre")
0.000.250.500.751.0010-0.50 10-0.25 100.00 100.25 100.50 N=8 Recurrence Coefficients Interpolating from Chebyshev to Legendregammabeta

The crosses denote the values of the β recursion coefficients for Chebyshev polynomials; the circles the β recursion coefficients for Legendre polynomials. The interpolating line in between stands for the β recursion coefficients of $w(t;\gamma)$.