Initial conditions for high-dimensional systems of ODEs that depend on a small number of parameters

Modeling
1

I was wondering how the Stan ODE solver deals with initial conditions that depend on parameters. Suppose we have a system of ODEs of dimension n, where n can be quite large (50, say). The initial state depends on the parameters y(t_0) = y_0(\theta) for parameter vector \theta \in \mathbb{R}^k where k \ll n.

When I integrate my system of ODEs in Stan with the sensitivity method:

vector[n] y0 = initial_state_func(theta);
array[m] vector[n] y = ode_rk45(my_ode_sys, y0, t0, ts, theta);

Is this internally converted to a system of n \times (k+n) ODEs, as there are k parameters and n initial conditions? Or does Stan “know” that y_0 depends only on k parameters?

In some cases, it is possible to transform the system of ODEs such that the initial state does not depend on the parameters.

transformed data {
    vector[n] z0 = rep_vector(1.0, n);
    // ...
}
model {
    array[m] vector[n] z = ode_rk45(my_transformed_ode_sys, z0, t0, ts, theta);
    array[m] vector[n] y;
    for ( i in 1:m ) {
        vector y = my_inverse_transform_func(z[i], theta);
    } 
}

Which might considerably speed things up. I can do some experiments, but was wondering if anyone had a simple answer.

Thanks!!

2

I don’t know a simple answer but for high-dimensional systems you should probably use the adjoint ODE solver.