Adaptive solvers not working properly?

[justinmcgrath]

When solving this model, the solver doesn't seem to choose appropriate time steps. When I change the error tolerances and max steps options, it doesn't do any better. Changing output_step_size to a smaller value gives more values in the results, but overall it's the same. Interpolating the drivers to many more points forces it to use a smaller steps size and that works. It doesn't seem like it's supposed to work this way though. Is something wrong with this, or am I missing something?

tibrary(BioCro)
library(nitrogenAssim)
library(lattice)

model = list(
    direct_modules = list(),
    differential_modules = list('nitrogenAssim:n_uptake_and_distribution'),
    ode_solver = list(type='boost_rkck54', output_step_size = 0.02, adaptive_rel_error_tol = 1e-10, adaptive_ab_error_tol=1e-10, adaptive_max_steps=20000),
    initial_values = list(
        #Root_Biomass = 0.6,
        Root_Biomass = 0.3,
        #N = 200,
        N = 290,
        #Rn = 0.0002,
        Rn = 7e-14,
        #Shoot_Biomass = 1.2,
        Shoot_Biomass = 0.27145,
        #Sn = 0.00014,
        Sn = 0.00316,
        R1 = .001,
        Seed_N_Conc = 0.001
    ),
    parameters = list(
         timestep=1,
         Cpm = 20,
         #RGr = 0.004,
         RGr = 0.00275,
         #k = 0.00003,
         k = 9e-5,
         #SGr = 0.0045,
         SGr = 0.00484,
         #Cpms = 55,
         Cpms = 65,
         #Vmax = 0.003,
         Vmax = 5.18e-5,
         #sV = 0.013,
         sV = 1.55e-5,
         #F1 = 0.0001,
         F1 = 1.3e-5,
         #F6 = 39,
         F6 = 0.85,
         #F7 = 20
         F7 = 835
    )
)

read_text_table = function(text, ...) {utils::read.table(textConnection(text), ...)}

volume = read_text_table(header=TRUE, text=
"time    Volume
24      0.004
96      0.00386667
191.76  0.00381111
192     0.004
264     0.00391667
359.76  0.00371667
360     0.004
432     0.0037
527.76  0.00266667
528     0.004
600     0.00245
695.76  0.00151667
696     0.004
768     0.0026
863.76  0.00082222
864     0.004
936     0.00133333
1031.76 0.0002
1032    0.004
1104    0.0023
1199.76 0.0011
1200    0.004
1272    0.0022
10000   0.0022")  # Assume previous value to simulate a longer period

interp_volume = as.data.frame(setNames(with(volume, approx(time, Volume, seq(24, max(volume$time), length=1000))), c('time', 'Volume')))  # Change the length here to 10000 to get the correct result.

N_in_moles = read_text_table(header=TRUE, text=
"time   N_in_moles
24      0.04
192     0.04
264     0.03876737
360     0.04
432     0.03540856
528     0.04
600     0.03039316
696     0.04
768     0.01227705
864     0.04
936     0.04232669
1032    0.04
1104    0.02666792
1200    0.04
1272    0.02600718"
)

N_fun = with(N_in_moles, approxfun(time, N_in_moles, rule=2))
interp_volume$N_in_moles = N_fun(interp_volume$time)

results = with(model, {run_biocro(
    initial_values,
    parameters,
    interp_volume,
    direct_modules,
    differential_modules,
    ode_solver,
    verbose=TRUE
)})

x11(); xyplot(N ~ doy, results, type='b')
x11(); xyplot(Shoot_Biomass ~ time, results, type='b')

[justinmcgrath]

It looks like increasing the number of interpolating points isn't making it more correct. It looks like it has an effect similar to increasing the range of the prediction. With 1000, it has the first part of the curve, with 10000 it has more of the curve, and with 30000 it's adding more of the asymptote. I'm not sure how that would happen.

[justinmcgrath]

This happens because it gets the time step from the "timestep" parameter. You need to calculate the time step from the interpolated data and set that as the "timestep". I think we've talked about this before. It would be nice to change how this works.