Integrals

Pylians provide routines to carry out numerical integrals in a more efficient way than scipy integrate. The philosophy is that to compute the integral \(\int_a^b f(x)dx\) the user passes the integrator two arrays, one with some values of \(x\) between a and b and another with the values of \(f(x)\) at those positions. The integrator will interpolate internally the input data to evaluate the function at an arbitrary position \(x\). Pylians implements in c the fortran odeint function and wrap in python through cython.

For instance, to compute \(\int_0^5 (3x^3+2x+5) dx\) one would do

import numpy as np
import integration_library as IL

# integral value, its limits and precision parameters
yinit = np.zeros(1, dtype=np.float64)
x1    = 0.0
x2    = 5.0
eps   = 1e-8
h1    = 1e-10
hmin  = 0.0

# integral method and integrand function
function = 'linear'
bins     = 1000
x        = np.linspace(x1, x2, bins)
y        = 3*x**2 + 2*x + 5

I = IL.odeint(yinit, x1, x2, eps, h1, hmin, x, y, function, verbose=True)[0]

The value of integral is stored in I. The odeint routine needs the following ingredients:

  • yinit. The value of the integral is stored in this variable. Should be a 1D double numpy array with one single element equal to 0. If several integrals are being computed sequentially this variable need to be declared for each integral.

  • x1. Lower limit of the integral.

  • x2. Upper limit of the integral.

  • eps. Maximum local relative error tolerated. Typically set its value to be 1e8-1e10 lower than the value of the integral. Verify the convergence of the results by studying the dependence of the integral on this number.

  • h1. Initial guess for the first time step (cannot be 0).

  • hmin. Minimum allowerd time step (can be 0).

  • verbose. Set it to True to print some information on the integral computation.

There are two main methods to carry out the integral, depending on how the interpolation is performed.

  • function = 'linear'. The function is evaluated by interpolating linearly the input values of \(x\) and \(y=f(x)\).

    • x. 1D double numpy array containing the input, equally spaced, values of \(x\).

    • y. 1D double numpy array containing the values of \(y=f(x)\) at the x array positions.

  • function = 'log'. The function is evaluated by interpolating logaritmically the input values of \(x\) and \(y=f(x)\).

    • x. 1D double numpy array containing the input, equally spaced in log, values of \(\log_{10}(x)\).

    • y. 1D double numpy array containing the values of \(y=f(x)\) at the x array positions.

An example of using the log-interpolation to compute the integral \(\int_1^2 e^x dx\) is this

import numpy as np
import integration_library as IL

# integral value, its limits and precision parameters
yinit = np.zeros(1, dtype=np.float64)
x1    = 1.0
x2    = 2.0
eps   = 1e-10
h1    = 1e-12
hmin  = 0.0

# integral method and integrand function
function = 'log'
bins     = 1000
x        = np.logspace(np.log10(x1), np.log10(x2), bins)
y        = np.exp(x)

I = IL.odeint(yinit, x1, x2, eps, h1, hmin, np.log10(x), y,
              function, verbose=True)[0]

Warning

Be careful when using the log-interpolation, since the code will crash if a zero or negative value is encounter.

The user can create its own function to avoid evaluating the integrand via interpolations. This function has to be placed in the file library/integration_library/integration_library.pyx (see linear and sigma functions as examples). After that, a new function call has to be created in the function odeint of that file (see linear and log as examples).