Computing Gauss-Turán qudrature rules

In this section we look at how to use the algorithm in GaussTuranQuadrature.jl to compute new quadrature rules. For the mathematical background of this algorithm see the theory section. We use the example presented in this article, with the functions

\[\begin{equation} \varphi_j(x) = x^{p_j}, \quad j = 1, \ldots, 2(s + 1)n \end{equation}\]

where

\[\begin{equation} p_j = \begin{cases} \frac{j}{2} - 1 \text{ if $j$ is even} \\ \frac{j - 1}{2} - \frac{1}{3} \text{ if $j$ is odd} \end{cases}, \end{equation}\]

and the trivial weight function $w(x) = 1$.

Computing the rule

Below is some setup of the inputs for computing the rule.

using TaylorDiff: derivatives

p(n, s) = [j % 2 == 0 ? j / 2 - 1 : (j - 1) / 2 - 1 / 3 for j in 1:(2 * (s + 1) * n)]

n = 5
s = 1

# Powers of the input functions
P = p(n, s)

# Integrals over (0, 1) of the input functions
rhs = @. 1 / (P + 1)

# The functions themselves and their derivatives
ϕ = (x, j) -> derivatives(x -> x^P[j], x, 1.0, Val(2s + 1)).value;

#3 (generic function with 1 method)

Now we can compute the rule.

using Optim
using TaylorDiff
using PreallocationTools
using GaussTuranQuadrature

I, res = GaussTuranComputeRule(ϕ, n, s, rhs)

I

Gauss-Turán quadrature rule on the interval (0.0, 1.0), with 5 nodes and derivatives up to order 2.
* Nodes:
5-element Vector{Float64}:
 0.0018183674053714562
 0.045953252552785824
 0.22342395164593046
 0.5808639413677223
 0.921794918354219
* Weights:
5×3 Matrix{Float64}:
 0.0121793   2.9136e-5    2.39045e-8
 0.0955793   0.000898476  2.00418e-5
 0.276105    0.00439912   0.00054227
 0.392025   -0.000255524  0.00162487
 0.224111   -0.0054717    0.000305731

And evaluate the rule.

# Create a function which computes exp(x) once and fills a tuple of length 2s+1 with that value
f = x -> (val = Base.exp(x); ntuple(_ -> val, 2s + 1))

evaluation = I(f)
error = abs(evaluation - (Base.exp(1) - 1))
error

4.440892098500626e-16

We can also have a look at the optimization statistics.

res

 * Status: success

 * Candidate solution
    Final objective value:     1.639804e-12

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 2.23e-06 ≰ 1.0e-08

 * Work counters
    Seconds run:   1  (vs limit Inf)
    Iterations:    139
    f(x) calls:    329
    ∇f(x) calls:   329

Computing Gauss-Turán qudrature rules

Computing the rule

Obtaining higher accuracy