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
then the definite integral (1) is the value of $x$ at $t=B$.
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
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