Tutorial 3. Estimating the multivariate integral

<< Tutorial 2. Estimating the volume of a body

Back to the list of tutorials

Tutorial 4. Handler functions >>

In this tutorial we are going to demonstrate yet another application of \((t, m, s)\)-nets which is numerical calculation of definite integrals over cuboid regions of \(s\)-dimensional space. Say, we are given a \(5\)-dimensional integral over a unit cube \(I^5\) where \(I = [0,1)\):

\[J(\overrightarrow{x}) = \int_{I^5} \frac{ x^{(1)} \sin{(x^{(2)}\pi)} }{ (x^{(3)} - 5) \cos{(x^{(4)}x^{(5)})} } d\overrightarrow{x}\]

The function \(f(\overrightarrow{x})\) inside the integral cannot be plotted since it depends on \(5\) variables, however, an important note here is that this function is Lipschitz-continuous over the whole \(I^5\). Knowing this, let us estimate the value of \(J(\overrightarrow{x})\) with the help of the formula:

\[J(\overrightarrow{x}) \approx \frac{1}{2^m} \sum_{i=1}^{2^m} f(N_i)\]

where \(N_i\) is the \(i\)-th point of a \((t, m, s)\)-net.

1. Basic workflow

We yet again assume the same file structure as in tutorial 1. This time, for the sake of variety we are going to utilise Sobol nets that are represented with the tms::Sobol class. These nets have the same constructor as Niederreiter nets, however, due to some inner changes in construction algorithms may potentially expose additional useful features.

Let us define a \((t, 30, 5)\)-net, that is, a five-dimensional net with \(2^{30} = 1,073,741,824\) points.

#include "libraries/tms-nets/tms-nets.hpp"
#include <cstdint> // needed for uint*_t types
 
int main()
{
    uint32_t const                      tms_net_point_count = 1UL << 30;
    tms::Sobol                          tms_net(30, 5);
 
    return 0;
}

2. Count values of function at each point of a net

We now implement a function f that will calculate the value of function \(f(\overrightarrow{x})\) at the given point of five-dimensional space.

#include "libraries/tms-nets/tms-nets.hpp"
#include <cmath>
#include <cstdint> // needed for uint*_t types
#define pi 3.14159265358979
 
long double f(tms::Point x)
{
    return x[0] * sin(x[1] * pi) / ( (x[2] - 5) * cos(x[3] * x[4]) );
}
 
int main()
{
    uint32_t const                      tms_net_point_count = 1UL << 30;
    tms::Sobol                          tms_net(30, 5);
 
    return 0;
}

Let us then apply our formula within the main function.

#include "libraries/tms-nets/tms-nets.hpp"
#include <cmath>
#include <cstdint> // needed for uint*_t types
#include <iostream>
#define pi 3.14159265358979
 
long double f(tms::Point x)
{
    return x[0] * sin(x[1] * pi) / ( (x[2] - 5) * cos(x[3] * x[4]) );
}
 
int main()
{
    uint32_t const                      tms_net_point_count = 1UL << 30;
    tms::Sobol                          tms_net(30, 5);
 
    long double                         answer = 0.0;
 
    for (uint32_t point_i = 0; point_i < tms_net_point_count; ++point_i)
        answer += f(tms_net.generate_point(point_i));
    answer /= tms_net_point_count;
 
    std::cout << "J(x) = " << answer << '\n';
 
    return 0;
}

The expected output is:

J(x) = -0.0757313

The analytical solution for \(J(\overrightarrow{x})\) cannot be expressed via elementary functions and is expressed as follows:

\[J(\overrightarrow{x}) = \frac{\log 4 - \log 5}{\pi} \int_0^1 \frac{1}{t} \log{\left( \frac{2}{\cot{\frac{t}{2}} - 1} + 1 \right)} dt \approx -0.07573128645316662\]

As we see, our net provides us with a great accuracy. You may try comparing its accuracy with accuracy of regular Monte Carlo method and see the difference on your own.

<< Tutorial 2. Estimating the volume of a body

Back to the list of tutorials

Tutorial 4. Handler functions >>