The MCSimMod package can be used to define ordinary
differential equation (ODE) models and solve initial value problems
(IVPs) for such models. To define an ODE model, one first creates a
“model specification” as a text string or a text file. (More details
about the structure of the model specification text are provided in a
separate vignette: “Model Specification”.) The model specification text
can then be used to create a model object (i.e., an instance of the
MCSimMod Model
class) via the
MCSimMod function createModel()
.
Consider, for example, a simple ODE model given by \[\begin{align} \frac{\textrm{d}}{\textrm{d}t}y(t) &= m, \\ y(0) &= y_0. \end{align}\] We can create a model object corresponding to this ODE model using a text string that both provides the ODE for the state variable, \(y\), and sets the (default) values of the model parameters to \(y_0 = 2\) and \(m = 0.5\).
First, we load the MCSimMod package as follows.
Then, we define the model string as follows.
mod_string <- "
States = {y};
y0 = 2;
m = 0.5;
Initialize {
y = y0;
}
Dynamics {
dt(y) = m;
}
End.
"
Next, we create a Model
object called model
using the following command.
Once the model object is created, we can “load” the model (so that it’s ready for use in a given R session) as follows.
Then, we can perform a simulation that provides results for the desired output times (\(t = 0, 0.1, 0.2, \ldots, 20.0\)) using the following commands.
The final command shown above performs the model simulation and
stores the simulation results in a “matrix” data structure called
out
. The first five rows of this data structure are shown
below. Note that the independent variable, which is \(t\) in the case of the linear model we’ve
created here, is always labeled “time” in the output data structure.
time | y |
---|---|
0.0 | 2.00 |
0.1 | 2.05 |
0.2 | 2.10 |
0.3 | 2.15 |
0.4 | 2.20 |
We can examine the parameter values and initial conditions that were used for this simulation with the following commands.
Finally, we can create a visual representation of the simulation results. For example, we can plot the value of \(y\) vs. time (\(t\)) using the following command.
# Plot simulation results.
plot(out[, "time"], out[, "y"],
type = "l", lty = 1, lwd = 2, xlab = "Time",
ylab = "y"
)