This is an old version, view current version.

16.4 Coding the DAE system function

The first step in coding an DAE system in Stan is defining the DAE residual function. The system functions require a specific signature so that the solvers know how to use them properly.

The first argument to the residual function is time, passed as a real; the second argument to the residual function is the system state \(y\), passed as a vector, the third argument to the residual function is the state derivative \(y'\), also passed as a vector. The residual function’s return value is a vector of the same size as state and stae derivatives. Additional arguments can be included in the residual function to pass other information into the solve (these will be passed through the function that starts the DAE solution). These argument can be parameters (in our example, the reaction rate coefficient \(\alpha\), \(\beta\), and \(\gamma\)), data, or any quantities that are needed to define the DAE.

The above reaction be coded using the following function in Stan (see the user-defined functions chapter for more information on coding user-defined functions).

 vector chem(real t, vector yy, vector yp,
                 real alpha, real beta, real gamma) {
    vector[3] res;
    res[1] = yp[1] + alpha * yy[1] - beta * yy[2] * yy[3];
    res[2] = yp[2] - alpha * yy[1] + beta * yy[2] * yy[3] + gamma * yy[2] * yy[2];
    res[3] = yy[1] + yy[2] + yy[3] - 1.0;
    return res;
  }
}

The function takes in a time t (a real), the system state yy (a vector), state derivative yp (a vector), as well as parameter alpha (a real), beta (a real), and gamma (a real). The function returns a vector of the residuals at time t. The DAE coded here does not explicitly depend on t, however one still needs to specify t as an argument.

Strict signature

The types in the DAE residual function are strict. The first argument is the time passed as a real, the second argument is the state passed as a vector, the third argument is the state derivative passed as a vector, and the return type is a vector. A model that does not have this signature will fail to compile. The fourth argument onwards can be any type, granted all the argument types match the types of the respective arguments in the solver call.

All of these are possible DAE signatures:

vector my_dae1(real t, vector y, vector yp, real a0);
vector my_dae2(real t, vector y, vector yp, array[] int a0, vector a1);
vector my_dae3(real t, vector y, vector yp, matrix a0, array[] real a1, row_vector a2);

but these are not allowed:

vector my_dae1(real t, array[] real y, vector yp); 
// Second argument is not a vector
array[] real my_dae2(real t, vector y, vector yp); 
// Return type is not a vector
vector my_dae3(real t, vector y); 
// First argument is not a real and missing the third argument