tests/testthat/test-example-3-1.R

In nlmixrdevelopment/RxODE: Facilities for Simulating from ODE-Based Models

## Example 3.1 from
## "Solving Differential Equations in R" by Soetaert et al (2012)
## https://cran.r-project.org/web/packages/diffEq/vignettes/ODEinR.pdf Example #1
rxodeTest(
  {
    df.test <- structure(c(0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.2, 6.4, 6.6, 6.8, 7, 7.2, 7.4, 7.6, 7.8, 8, 8.2, 8.4, 8.6, 8.8, 9, 9.2, 9.4, 9.6, 9.8, 10, 10.2, 10.4, 10.6, 10.8, 11, 11.2, 11.4, 11.6, 11.8, 12, 12.2, 12.4, 12.6, 12.8, 13, 13.2, 13.4, 13.6, 13.8, 14, 14.2, 14.4, 14.6, 14.8, 15, 15.2, 15.4, 15.6, 15.8, 16, 16.2, 16.4, 16.6, 16.8, 17, 17.2, 17.4, 17.6, 17.8, 18, 18.2, 18.4, 18.6, 18.8, 19, 19.2, 19.4, 19.6, 19.8, 20, 2, 2.339221, 2.716438, 3.129638, 3.574846, 4.046088, 4.535599, 5.034256, 5.532232, 6.019753, 6.487848, 6.928975, 7.337435, 7.709526, 8.043484, 8.339239, 8.59808, 8.822277, 9.014729, 9.178661, 9.317383, 9.43412, 9.531896, 9.613469, 9.681302, 9.737552, 9.784099, 9.822541, 9.85424, 9.880346, 9.901823, 9.919476, 9.933976, 9.945878, 9.955645, 9.963656, 9.970225, 9.975609, 9.980021, 9.983637, 9.986599, 9.989026, 9.991013, 9.992641, 9.993974, 9.995066, 9.99596, 9.996692, 9.997292, 9.997782, 9.998184, 9.998514, 9.998784, 9.999004, 9.999185, 9.999332, 9.999453, 9.999552, 9.999633, 9.9997, 9.999754, 9.999798, 9.999834, 9.999864, 9.999888, 9.999909, 9.999926, 9.999941, 9.999951, 9.99996, 9.999967, 9.999973, 9.999978, 9.999982, 9.999986, 9.999988, 9.99999, 9.999992, 9.999993, 9.999994, 9.999995, 9.999996, 9.999997, 9.999997, 9.999998, 9.999998, 9.999999, 9.999999, 9.999999, 9.999999, 9.999999, 9.999999, 10, 10, 10, 10, 10, 10, 10, 10, 10, 0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.2, 6.4, 6.6, 6.8, 7, 7.2, 7.4, 7.6, 7.8, 8, 8.2, 8.4, 8.6, 8.8, 9, 9.2, 9.4, 9.6, 9.8, 10, 10.2, 10.4, 10.6, 10.8, 11, 11.2, 11.4, 11.6, 11.8, 12, 12.2, 12.4, 12.6, 12.8, 13, 13.2, 13.4, 13.6, 13.8, 14, 14.2, 14.4, 14.6, 14.8, 15, 15.2, 15.4, 15.6, 15.8, 16, 16.2, 16.4, 16.6, 16.8, 17, 17.2, 17.4, 17.6, 17.8, 18, 18.2, 18.4, 18.6, 18.8, 19, 19.2, 19.4, 19.6, 19.8, 20, 12, 11.580181, 11.257713, 11.006774, 10.809502, 10.653184, 10.528525, 10.428613, 10.348214, 10.283304, 10.230763, 10.188146, 10.153517, 10.125341, 10.102387, 10.083672, 10.068401, 10.055933, 10.045748, 10.037425, 10.03062, 10.025056, 10.020505, 10.016782, 10.013736, 10.011243, 10.009203, 10.007533, 10.006167, 10.005049, 10.004133, 10.003383, 10.00277, 10.002268, 10.001857, 10.00152, 10.001244, 10.001018, 10.000833, 10.000682, 10.000559, 10.000458, 10.000375, 10.000307, 10.000252, 10.000206, 10.000169, 10.000139, 10.000114, 10.000093, 10.000076, 10.000062, 10.00005, 10.000041, 10.000034, 10.000027, 10.000022, 10.000018, 10.000015, 10.000012, 10.00001, 10.000008, 10.000007, 10.000006, 10.000005, 10.000004, 10.000003, 10.000003, 10.000002, 10.000002, 10.000001, 10.000001, 10.000001, 10.000001, 10.000001, 10.000001, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.2, 6.4, 6.6, 6.8, 7, 7.2, 7.4, 7.6, 7.8, 8, 8.2, 8.4, 8.6, 8.8, 9, 9.2, 9.4, 9.6, 9.8, 10, 10.2, 10.4, 10.6, 10.8, 11, 11.2, 11.4, 11.6, 11.8, 12, 12.2, 12.4, 12.6, 12.8, 13, 13.2, 13.4, 13.6, 13.8, 14, 14.2, 14.4, 14.6, 14.8, 15, 15.2, 15.4, 15.6, 15.8, 16, 16.2, 16.4, 16.6, 16.8, 17, 17.2, 17.4, 17.6, 17.8, 18, 18.2, 18.4, 18.6, 18.8, 19, 19.2, 19.4, 19.6, 19.8, 20, 2, 2.339223, 2.716446, 3.129649, 3.574856, 4.046097, 4.535606, 5.034264, 5.53224, 6.019761, 6.487856, 6.928986, 7.337447, 7.70954, 8.043498, 8.339252, 8.598091, 8.822286, 9.014735, 9.178665, 9.317385, 9.43412, 9.531895, 9.613468, 9.681301, 9.737555, 9.784101, 9.822543, 9.854241, 9.880347, 9.901823, 9.919476, 9.933977, 9.94588, 9.955647, 9.963657, 9.970225, 9.97561, 9.980022, 9.983637, 9.986599, 9.989026, 9.991013, 9.992641, 9.993974, 9.995066, 9.99596, 9.996692, 9.997292, 9.997782, 9.998184, 9.998513, 9.998783, 9.999003, 9.999184, 9.999332, 9.999453, 9.999552, 9.999633, 9.9997, 9.999754, 9.999799, 9.999835, 9.999865, 9.99989, 9.99991, 9.999926, 9.999939, 9.99995, 9.999959, 9.999967, 9.999973, 9.999978, 9.999982, 9.999985, 9.999988, 9.99999, 9.999992, 9.999993, 9.999995, 9.999995, 9.999996, 9.999997, 9.999998, 9.999998, 9.999998, 9.999999, 9.999999, 9.999999, 9.999999, 9.999999, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.2, 6.4, 6.6, 6.8, 7, 7.2, 7.4, 7.6, 7.8, 8, 8.2, 8.4, 8.6, 8.8, 9, 9.2, 9.4, 9.6, 9.8, 10, 10.2, 10.4, 10.6, 10.8, 11, 11.2, 11.4, 11.6, 11.8, 12, 12.2, 12.4, 12.6, 12.8, 13, 13.2, 13.4, 13.6, 13.8, 14, 14.2, 14.4, 14.6, 14.8, 15, 15.2, 15.4, 15.6, 15.8, 16, 16.2, 16.4, 16.6, 16.8, 17, 17.2, 17.4, 17.6, 17.8, 18, 18.2, 18.4, 18.6, 18.8, 19, 19.2, 19.4, 19.6, 19.8, 20, 12, 11.580174, 11.257712, 11.006774, 10.809504, 10.653181, 10.528522, 10.428611, 10.348211, 10.283303, 10.230764, 10.188146, 10.153518, 10.125341, 10.102388, 10.083673, 10.068402, 10.055933, 10.045748, 10.037424, 10.03062, 10.025055, 10.020504, 10.016781, 10.013735, 10.011243, 10.009203, 10.007533, 10.006167, 10.005048, 10.004133, 10.003384, 10.00277, 10.002268, 10.001857, 10.00152, 10.001244, 10.001019, 10.000834, 10.000683, 10.000559, 10.000458, 10.000375, 10.000307, 10.000251, 10.000206, 10.000168, 10.000138, 10.000113, 10.000092, 10.000076, 10.000062, 10.000051, 10.000042, 10.000034, 10.000028, 10.000023, 10.000019, 10.000015, 10.000013, 10.00001, 10.000008, 10.000007, 10.000006, 10.000005, 10.000004, 10.000003, 10.000003, 10.000002, 10.000002, 10.000001, 10.000001, 10.000001, 10.000001, 10.000001, 10.000001, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10), .Dim = c(101L, 8L), .Dimnames = list(NULL, c("time", "y", "time", "y", "time", "y", "time", "y")))

    ms <- c("liblsoda", "lsoda", "dop853")
    ode <- RxODE(model = "d/dt(y) = r * y * (1.0 - y/K);")
    .rxWithOptions(list(digits = 6), {
      for (meth in ms) {
        context(sprintf("Example 3-1 (%s)", meth))

        ## create event table with times at which we observe the system
        et <- eventTable(time.units = NA)
        et$add.sampling(seq(from = 0, to = 20, by = 0.2))

        ## same model, different initial values (2 and 12)
        out1 <- ode$solve(params = c(r = 1, K = 10), events = et, inits = c(y = 2))
        out2 <- ode$solve(params = c(r = 1, K = 10), events = et, inits = c(y = 12))

        ## matplot(x = out1[,1], y = cbind(out1[,2], out2[,2]), type = "l",
        ##    main = "logistic growth", xlab="time", ylab="")

        ## Now use a non-stiff solver
        out1.ns <- NULL
        test_that("stiff option adds warning", {
          expect_warning({
            out1.ns <<- ode$solve(params = c(r = 1, K = 10), events = et, inits = c(y = 2), stiff = FALSE)
          })
        })
        out2.ns <-
          ode$solve(params = c(r = 1, K = 10), events = et, inits = c(y = 12), method = "dop853")

        df <- round(cbind(out1, out2, out1.ns, out2.ns), 3)

        test_that("Runs example 3.1 correctly", {
          expect_equal(df, round(df.test, 3))
        })

        head(cbind(out1, out2, out1.ns, out2.ns), n = 15)
      }
    })
  },
  silent = TRUE,
  test = "lvl2"
)