(1) Introduction to MCSimMod

Dustin Kapraun

2025-04-09

The MCSimMod package can be used to define ordinary differential equation (ODE) models and solve initial value problems for such models. Models can be specified using text that that is supplied either as a string or in a file. The text is used to create a model object (i.e., an instance of the MCSimMod Model class). The R commands that follow create a model object corresponding to a simple ODE model given by \[\begin{align} \frac{\textrm{d}}{\textrm{d}t}y(t) &= m, \\ y(0) &= y_0, \end{align}\] using a 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\). (More details about the structure of the model specification text are provided in a separate tutorial.)

First, we load the MCSimMod package as follows.

library(MCSimMod)

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.

model <- createModel(mString = mod_string)

Once the model object is created, we can “load” the model (so that it’s ready for use in a given R session) and perform a simulation that provides results for the desired output times (\(t = 0, 0.1, 0.2, \ldots, 20.0\)) using the following commands.

model$loadModel()
#> C compilation complete. Full details are available in the file C:\Users\dkapraun\AppData\Local\Temp\Rtmp4YuNTn\compiler_output.txt.
#> Hash created and saved in the file C:\Users\dkapraun\AppData\Local\Temp\Rtmp4YuNTn\mcsimmod_369c7a4573eb_model.md5.
times <- seq(from = 0, to = 20, by = 0.1)
out <- model$runModel(times)

The final command shown above performs a 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.

We can examine the parameter values and initial conditions that were used for this simulation with the following commands.

model$parms
#>  y0   m 
#> 2.0 0.5
model$Y0
#> y 
#> 2

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"
)