Computing definite integrals

In this short tutorial, we will show how it is possible to use heyoka.py to compute definite integrals.

Consider the integral

\[ \int_A^B f\left( t \right)\,dt, \tag{1} \]

where \(f\left( t \right)\) is a differentiable function of \(t\) and \(A, B\in\mathbb{R}\). This integral can be seen as as the solution of the time-dependent differential equation

\[ \frac{dx}{dt} = f\left( t \right), \tag{2} \]

where \(x\) is a dummy state variable. Indeed, the integration of (2) by quadrature between \(t=A\) and \(t=B\) yields:

\[ x\left(B\right) - x\left(A\right) = \int_A^B f\left(t\right) \, dt. \]

Note that we are always free to choose \(x\left( A \right) = 0\), because the dynamics of \(x\) in (2) does not depend on the value of \(x\) itself. Thus, provided that we set up a numerical integration of (2) in which

• we set \(t=A\) as initial time coordinate,

• we set \(x\left( A \right) = 0\) as initial condition,

then the definite integral (1) is the value of \(x\) at \(t=B\).

Examples

Let us start easy:

import heyoka as hy, numpy as np
x = hy.make_vars("x")

# Integrate sin(t) between 0 and pi.
ta = hy.taylor_adaptive([(x, hy.sin(hy.time))], [0.])
ta.propagate_until(np.pi)

# Print the result.
ta.state[0]

2.0

As expected, \(\int_0^\pi \sin t\, dt = 2\). Let’s try to change the integration range:

# Reset the state.
ta.state[0] = 0

# New integration limits: from 1 to 2.
ta.time = 1
ta.propagate_until(2.)

ta.state[0]

0.9564491424152821

Indeed, \(\int_1^2 \sin t\, dt = \cos 1 - \cos 2 = 0.9564491424152821\ldots\).

Let us try with a more complicated function:

\[ \int_\sqrt{2}^\sqrt{3} \frac{\sin \left( \cos t \right)\cdot \operatorname{erf}{t}}{\log\left(\sqrt{t}\right)}\,dt. \]

ta = hy.taylor_adaptive([(x, hy.sin(hy.cos(hy.time))*hy.erf(hy.time)/hy.log(hy.sqrt(hy.time)))],
                        [0.], time = np.sqrt(2))
ta.propagate_until(np.sqrt(3))
ta.state[0]

0.012382281847117892

The result matches the value produced by Wolfram Alpha.

Note that, since heyoka.py supports integration backwards in time, flipping around the integration limits also works as expected:

ta.state[0] = 0
ta.time = np.sqrt(3)
ta.propagate_until(np.sqrt(2))
ta.state[0]

-0.012382281847117878

Let us also perform the integration in extended precision:

ta = hy.taylor_adaptive([(x, hy.sin(hy.cos(hy.time))*hy.erf(hy.time)/hy.log(hy.sqrt(hy.time)))],
                        [np.longdouble(0)], time = np.sqrt(np.longdouble(3)), fp_type=np.longdouble)
ta.propagate_until(np.sqrt(np.longdouble(2)))
ta.state[0]

-0.012382281847117883605

Limitations and caveats

This method for the computation of definite integrals inherits all the peculiarities and caveats of heyoka.py. For instance, the computation will fail if the derivative of \(f\left( t \right)\) becomes infinite within the integration interval:

# Cannot compute the area of a semi-circle!
ta = hy.taylor_adaptive([(x, hy.sqrt(1 - hy.time**2))],
                        [0.], time = -1.)
ta.propagate_until(1.)
ta.state[0]

nan