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')
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?