Ctrl+K

Computing definite integrals

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
Contents