8.1.2. 1d integration

Let us now integrate the following function:

(1)\[f(x) = a\sin(x) + b\cos(kx)\]

in 30 intervals in a range \((-\pi,\pi)\). The analytic integral we will use as a cross check reads as follows:

(2)\[\int f(x)\,dx = -a\cos(x) + \frac{b}{k}\sin(kx) + C.\]

The example script is below. We will use the Gauss-Legendre quadrature rule with 3 points for each bin.

01_integral1d.py

import load
from gna.env import env
import gna.constructors as C
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec
from gna.bindings import common
import ROOT as R

# Make a variable for global namespace
ns = env.globalns

# Create a parameter in the global namespace
pa = ns.defparameter('a', central=1.0, free=True, label='weight 1')
pb = ns.defparameter('b', central=0.2, free=True, label='weight 2')
pk = ns.defparameter('k', central=8.0, free=True, label='argument scale')

# Print the list of parameters
ns.printparameters(labels=True)
print()

# Define binning and integration orders
nbins = 30
x_edges = np.linspace(-1.0*np.pi, 1.0*np.pi, nbins+1, dtype='d')
orders = 3

# Initialize integrator
integrator = R.IntegratorGL(nbins, orders, x_edges)
int_points = integrator.points.x

# Create integrable: a*sin(x) + b*cos(k*x)
cos_arg = C.WeightedSum( ['k'], [int_points] )
sin_t = R.Sin(int_points)
cos_t = R.Cos(cos_arg.sum.sum)
fcn   = C.WeightedSum(['a', 'b'], [sin_t.sin.result, cos_t.cos.result])

integrator.hist.f(fcn.sum.sum)

# Print objects
integrator.print()
print()

fcn.print()
print()

cos_t.print()
print()

# Label transformations
integrator.points.setLabel('Sampler\n(Gauss-Legendre)')
integrator.hist.setLabel('Integrator\n(convolution)')
cos_arg.sum.setLabel('kx')
sin_t.sin.setLabel('sin(x)')
cos_t.cos.setLabel('cos(kx)')
fcn.sum.setLabel('a sin(x) + b cos(kx)')

# Compute the integral analytically as a cross check
# int| a*sin(x) + b*cos(k*x) = a*cos(x) - b/k*sin(kx) + C
a, b, k = (p.value() for p in (pa, pb, pk))
F = - a*np.cos(x_edges) + (b/k)*np.sin(k*x_edges)
integral_a = F[1:] - F[:-1]
bar_width = (x_edges[1:]-x_edges[:-1])*0.5
bar_left = x_edges[:-1]+bar_width

integral_n = integrator.hist.hist.data()

# Do some plotting
fig = plt.figure()
gs = gridspec.GridSpec(nrows=4, ncols=1, hspace=0.0)
ax = plt.subplot(gs[:3,0])
ax.minorticks_on()

ax.set_ylabel( 'f(x)' )
ax.set_title(r'$a\,\sin(x)+b\,\sin(kx)$')

# Draw the function and integrals
fcn.sum.sum.plot_vs(int_points, 'o-', label='function at integration points', alpha=0.5, markerfacecolor='none', markersize=4.0)
integrator.hist.hist.plot_bar(label='integral: numerical', divide=2, shift=0)
ax.bar(bar_left, integral_a, bar_width, label='integral: analytical', align='edge')

# Freeze axis limits and draw bin edges
ax.autoscale(enable=True, axis='y')
ymin, ymax = ax.get_ylim()
ax.vlines(integrator.points.xedges.data(), ymin, ymax, linestyle='--', alpha=0.4, linewidth=0.5)

ax.legend(loc='lower right')

# Add difference
ax = plt.subplot(gs[3,0], sharex=ax)
ax.set_xlabel('x')
diff_factor=1.e-7
ax.set_ylabel(r'diff., $\times10^{-7}$')

diff = integral_n-integral_a
ax.bar(x_edges[:-1], diff/diff_factor, bar_width*2.0, align='edge')

# Freeze axis limits and draw bin edges
ax.autoscale(enable=True, axis='y')
ymin, ymax = ax.get_ylim()
ax.vlines(integrator.points.xedges.data(), ymin, ymax, linestyle='--', alpha=0.4, linewidth=0.5)

plt.show()

Before going into details let us first look at the graph that follows the procedure described in Introduction.

Computatinal graph used to integrate function \(f(x)=a\sin(x)+b\cos(kx)\).

The sampler provide an array of \(x\), all the points needed to compute the integral. The points are used as input to the chain of 4 transformations implementing the integrand function (1). The output of the function is then passed to the integrator, that does the convolution and produces the histogram. The unbound sampler output contains the bins’ edges.

We now skip the definition of the parameters as it was done in previous examples and go to the initialization of the integrator.

nbins = 30
x_edges = np.linspace(-1.0*np.pi, 1.0*np.pi, nbins+1, dtype='d')
orders = 3
integrator = R.IntegratorGL(nbins, orders, x_edges)
int_points = integrator.points.x

Here we define the 31 edges in a range \((-\pi,\pi)\) for 30 bins and choose the integration orders (3). The orders may be specified for each bin in separate, in this case the orders variable should contain a numpy array of size 30 with integers. Since we are using raw [1] constructor we explicitly have specified the data type ‘d’ for bin edges.

In the last line we get the output with integration points.

The function (1) is defined as follows:

cos_arg = C.WeightedSum( ['k'], [int_points] )
sin_t = R.Sin(int_points)
cos_t = R.Cos(cos_arg.sum.sum)
fcn   = C.WeightedSum(['a', 'b'], [sin_t.sin.result, cos_t.cos.result])

We create multiplication by \(k\) as WeightedSum with one item: integration points (\(x\)) weighted by \(k\). Then we create \(\sin\) and \(\cos\) and sum them in another instance of WeightedSum. The last step is to pass the output of the function to the integrator:

integrator.hist.f(fcn.sum.sum)

That is it. The integration chain is ready. This chain will react to the change of the variables. Moreover, if, for example, only the value of \(a\) is changed, the branch with \(cos(kx)\) will not be executed since it has the cached result for all the values of \(x\). In this case only \(sin(x)\), the following WeightedSum and the integrator will be triggered.

The status of the integrator object with two transformations is the following:

[obj] IntegratorGL: 2 transformation(s), 0 variables
 0 [trans] points: 0 input(s), 2 output(s)
     0 [out] x: array 1d, shape 90, size  90
     1 [out] xedges: array 1d, shape 31, size  31
 1 [trans] hist: 1 input(s), 1 output(s)
     0 [in]  f <- [out] sum: array 1d, shape 90, size  90
     0 [out] hist: hist,  30 bins, edges -3.14159265359->3.14159265359, width 0.209439510239

31 bin edges, 30 bins, 3 points per bin: 90 integration points in total.

In order to cross check the integrator let us also implement the analytical solution for antiderivative \(F\) from (2).

a, b, k = (p.value() for p in (pa, pb, pk))
F = - a*np.cos(x_edges) + (b/k)*np.sin(k*x_edges)
integral_a = F[1:] - F[:-1]

Here we obtain the values of the parameters, compute the antiderivative for each bin edge and then compute the integrals as differences. The result our the computations is shown on the figure below.

The Gauss-Legendre quadrature application to the function (1).

The line-marker plot represents the value of the function computed at the integration points. Dashed lines show the bin edges. One may see that since Gauss-Legendre quadrature is used the integration points are located unevenly within bins. Two integrals are plotted as bar charts stacked left-to-right: the output of the integrator and the analytical solution. The lower plot shows the difference between the integrator output and the solution.

8.1.2.1. Different integrators

For 1d case there are 3 integrators implemented: Gauss-Legendre, rectangular and trapezoidal. They may be used in a similar fashion:

integrator = R.IntegratorGL(nbins, orders, x_edges)
integrator = R.IntegratorRect(nbins, orders, x_edges)
integrator = R.IntegratorTrap(nbins, orders, x_edges)

all the other procedures are similar. The rectangular integration has an extra option that specifies the rectangles’ alignment within the interval. The default is center.

integrator = R.IntegratorRect(nbins, orders, x_edges, 'left')
integrator = R.IntegratorRect(nbins, orders, x_edges, 'center')
integrator = R.IntegratorRect(nbins, orders, x_edges, 'right')

A final note on the integration order. For each integrator order means the number of points needed to integrate a single bin. Therefore for \(N\) bins each integrated with order of \(O_i\) the total number of integration points will be:

\[N_p = \sum\limits_i^N O_i.\]

In case of the trapezoidal rule \(N-1\) points are on the adjacent bin edges. Thus the total number of points is:

\[N_p = \sum\limits_i^N O_i - (N-1).\]

[1]

The constructors defined in C++ require the C++ datatypes to be passed. Numpy arrays should have the appropriate datatype specified, ‘d’ for doubles and ‘i’ for integers. On the other hand, the constructors defined in the gna.constructors library do implicit conversion.