# Numerical Integration

The goal of this tutorial is to solve an integral using Gauss quadrature,

$$$I := \int_{0}^{1} f(t) \mathrm{d} t \approx \sum_{k=1}^n w_k f(t_k),$$$

where we choose $f(t) = \sin(t)$, and $n = 5$.

Make sure to check out this tutorial too.

### Variant 0

julia> using PolyChaos

julia> n = 5;

julia> f(t) = sin(t);

julia> variant0 = integrate(f, op)
0.4596976941320484

julia> print("Numerical error: $(abs(1 - cos(1) - variant0))") Numerical error: 1.8868240303504535e-13 with negligible numerical errors. ### Variant 1 Let us now solve the same problem, while elaborating what is going on under the hood. At first, we load the package by calling julia> using PolyChaos Now we define a measure, specifically the uniform measure$\mathrm{d}\lambda(t) = w(t) \mathrm{d} t$with the weight$w$defined as $$$w: \mathcal{W} = [0,1] \rightarrow \mathbb{R}, \quad w(t) = 1.$$$ This measure can be defined using the composite type Uniform01Measure: julia> measure = Uniform01Measure() Uniform01Measure(PolyChaos.w_uniform01, (0.0, 1.0), true) Next, we need to compute the quadrature rule relative to the uniform measure. To do this we use the composite type Quad. julia> quadRule1 = Quad(n-1, measure) ┌ Warning: For measures of type Uniform01Measure the quadrature rule should be based on the recurrence coefficients. └ @ PolyChaos ~/Documents/Code/JuliaDev/PolyChaos/src/typesQuad.jl:58 Quad{Float64,Array{Float64,1}}("quadgp", 4, [1.0, 0.8535533905932737, 0.5, 0.14644660940672627, 0.0], [0.033333333333333354, 0.26666666666666666, 0.4, 0.26666666666666666, 0.033333333333333354]) julia> nw(quadRule1) 5×2 Array{Float64,2}: 1.0 0.0333333 0.853553 0.266667 0.5 0.4 0.146447 0.266667 0.0 0.0333333 This creates a quadrature rule quadRule_1 relative to the measure measure. The function nw() prints the nodes and weights. To solve the integral we call integrate() julia> variant1 = integrate(f, quadRule1) 0.4596977927043755 julia> print("Numerical error:$(abs(1 - cos(1) - variant1))")
Numerical error: 9.857251526135258e-8

### Revisiting Variant 0

Why is the error from variant 0 so much smaller? It's because the quadrature rule for variant 0 is based on the recurrence coefficients of the polynomials that are orthogonal relative to the measure measure. Let's take a closer look First, we compute the orthogonal polynomials using the composite type OrthoPoly, and we set the keyword addQuadrature to false.

julia> op = Uniform01OrthoPoly(n, addQuadrature=false)
Uniform01OrthoPoly{Array{Float64,1},Uniform01Measure,EmptyQuad{Float64}}(5, [0.5, 0.5, 0.5, 0.5, 0.5, 0.5], [1.0, 0.08333333333333333, 0.06666666666666667, 0.06428571428571428, 0.06349206349206349, 0.06313131313131314], Uniform01Measure(PolyChaos.w_uniform01, (0.0, 1.0), true), EmptyQuad{Float64}())

Note how op has a field EmptyQuad, i.e. we computed no quadrature rule. The resulting system of orthogonal polynomials is characterized by its recursion coefficients $(\alpha, \beta)$, which can be extracted with the function coeffs().

julia> coeffs(op)
6×2 Array{Float64,2}:
0.5  1.0
0.5  0.0833333
0.5  0.0666667
0.5  0.0642857
0.5  0.0634921
0.5  0.0631313

Now, the quadrature rule can be constructed based on op, and the integral be solved.

julia> quadRule2 = Quad(n, op)
Quad{Float64,Array{Float64,1}}("golubwelsch", 5, [0.046910077030667935, 0.23076534494715842, 0.49999999999999994, 0.7692346550528418, 0.9530899229693321], [0.11846344252809445, 0.23931433524968332, 0.28444444444444444, 0.23931433524968337, 0.1184634425280949])

5×2 Array{Float64,2}:
0.0469101  0.118463
0.230765   0.239314
0.5        0.284444
0.769235   0.239314
0.95309    0.118463

0.4596976941320484

julia> print("Numerical error: \$(abs(1 - cos(1) - variant0_revisited))")
Numerical error: 1.8818280267396403e-13

### Comparison

We see that the different variants provide slightly different results:

julia> 1-cos(1) .- [variant0 variant1 variant0_revisited]
1×3 Array{Float64,2}:
-1.88183e-13  -9.85725e-8  -1.88183e-13

with variant0 and variant0_revisited being the same and more accurate than variant1. The increased accuracy is based on the fact that for variant0 and variant0_revisited the quadrature rules are based on the recursion coefficients of the underlying orthogonal polynomials. The quadrature for variant1 is based on an general-purpose method that can be significantly less accurate, see also the next tutorial.