Runtime parameters#
Added in version 0.2.0.
In all the examples we have seen so far, numerical constants have always been hard-coded to fixed values when constructing the expressions of the ODEs.
While this approach leads to efficient code generation by embedding the values of the constants directly in the source code and by opening up further optimisation opportunities for the compiler, on the other hand it requires the construction of a new integrator if the values of the numerical constants change. In some applications (e.g., parametric studies), this overhead can become a performance bottleneck, especially for large ODE systems.
In order to avoid having to re-create a new integrator if the value of a constant
in an expression changes, heyoka’s expression system provides a node type,
called param (or runtime parameter), which
represents mathematical constants whose value is not known at the time of
construction of the expression.
In this section, we will illustrate the usage of runtime parameters in heyoka via a pendulum ODE system in which the values of the gravitational constant \(g\) and of the pendulum’s length \(l\) are left undetermined at the time of the definition of the system:
Let’s start with the construction of the integrator:
// Create the symbolic variables x and v.
auto [x, v] = make_vars("x", "v");
// Create the integrator object
// in double precision.
auto ta = taylor_adaptive<double>{// Definition of the ODE system:
// x' = v
// v' = -g/l * sin(x)
{prime(x) = v, prime(v) = -par[0] / par[1] * sin(x)},
// Initial conditions
// for x and v.
{0.05, 0.},
// Values of the runtime parameters
// g and l. If not provided,
// all parameters will be inited
// to zero.
kw::pars = {9.8, 1.}};
With respect to the previous examples, where \(g/l\) had been
hard-coded to 9.8, now \(g/l\) is represented as par[0] / par[1].
The syntax par[i] indicates a runtime parameter that is stored
at the index i in an array of parameter values. The array of
parameter values is optionally passed to the constructor as the keyword argument
kw::pars (which, in this case, contains the values 9.8 for \(g\) and 1.
for \(l\)).
If an array of parameter values is not passed to the constructor,
heyoka will infer that the ODE system contains 2 parameters and will then initialise
the array of parameter values with zeroes.
If we try to print the integrator to screen, the output will confirm the presence of runtime parameters:
Tolerance : 2.2204460492503131e-16
Taylor order : 20
Dimension : 2
Time : 0.0000000000000000
State : [0.050000000000000003, 0.0000000000000000]
Parameters : [9.8000000000000007, 1.0000000000000000]
Note that the array of parameter values, like the state vector and the time coordinate, is stored as a data member in the integrator object.
The period of a pendulum of length \(1\,\mathrm{m}\) on Earth is \(\sim 2\,\mathrm{s}\) in the small-angle approximation, as it can be confirmed via numerical integration:
// Integrate for ~1 period.
ta.propagate_until(2.0074035758801299);
std::cout << ta << '\n';
Tolerance : 2.2204460492503131e-16
Taylor order : 20
Dimension : 2
Time : 2.0074035758801299
State : [0.050000000000000003, 7.5784060331002885e-17]
Parameters : [9.8000000000000007, 1.0000000000000000]
As you can see, after 1 period the state of the system went back to the initial conditions.
We are now going to move to Mars, where the gravitational acceleration on the surface is \(\sim 3.72\,\mathrm{m}/\mathrm{s}^2\) (instead of Earth’s \(\sim 9.8\,\mathrm{m}/\mathrm{s}^2\)). First we reset the time coordinate:
// Reset the time coordinate.
ta.set_time(0.);
Then we change the value of the gravitational constant \(g\), which, as explained above, is stored at index 0 in the array of parameter values:
// Change the value of g.
ta.get_pars_data()[0] = 3.72;
Note that, like the for the state data, the get_pars_data()
function returns a naked pointer
that can be used to modify the parameter values. Another function
of the integrator object, get_pars(),
returns a const reference to the std::vector
holding the parameter values.
Because gravity is weaker on Mars, the period of a \(1\,\mathrm{m}\) pendulum increases to \(\sim 3.26\,\mathrm{s}\). We can confirm this via numerical integration:
// Integrate for ~1 period on Mars.
ta.propagate_until(3.2581889116828258);
std::cout << ta << '\n';
Tolerance : 2.2204460492503131e-16
Taylor order : 20
Dimension : 2
Time : 3.2581889116828258
State : [0.050000000000000003, 2.1864533707994132e-16]
Parameters : [3.7200000000000002, 1.0000000000000000]
Full code listing#
#include <iostream>
#include <heyoka/heyoka.hpp>
using namespace heyoka;
int main()
{
// Create the symbolic variables x and v.
auto [x, v] = make_vars("x", "v");
// Create the integrator object
// in double precision.
auto ta = taylor_adaptive<double>{// Definition of the ODE system:
// x' = v
// v' = -g/l * sin(x)
{prime(x) = v, prime(v) = -par[0] / par[1] * sin(x)},
// Initial conditions
// for x and v.
{0.05, 0.},
// Values of the runtime parameters
// g and l. If not provided,
// all parameters will be inited
// to zero.
kw::pars = {9.8, 1.}};
// Let's print to screen the integrator object.
std::cout << ta << '\n';
// Integrate for ~1 period.
ta.propagate_until(2.0074035758801299);
std::cout << ta << '\n';
// Reset the time coordinate.
ta.set_time(0.);
// Change the value of g.
ta.get_pars_data()[0] = 3.72;
// Integrate for ~1 period on Mars.
ta.propagate_until(3.2581889116828258);
std::cout << ta << '\n';
}