Demonstration that Sundials and DifferentialEquations.jl generated data are within 9 digits of accuracy.

[ChrisRackauckas]

using ParameterEstimation
using ModelingToolkit, DifferentialEquations
solver = Tsit5()

@parameters mu_N mu_EE mu_LE mu_LL mu_M mu_P mu_PE mu_PL delta_NE delta_EL delta_LM rho_E rho_P
@variables t N(t) E(t) S(t) M(t) P(t) y1(t) y2(t) y3(t) y4(t)
D = Differential(t)
states = [N, E, S, M, P]
parameters = [
    mu_N,
    mu_EE,
    mu_LE,
    mu_LL,
    mu_M,
    mu_P,
    mu_PE,
    mu_PL,
    delta_NE,
    delta_EL,
    delta_LM,
    rho_E,
    rho_P,
]
@named model = ODESystem([
                             D(N) ~ -N * mu_N - N * P * delta_NE,
                             D(E) ~ N * P * delta_NE - E^2 * mu_EE -
                                    E * delta_EL + E * P * rho_E,
                             D(S) ~ S * delta_EL - S * delta_LM - S^2 * mu_LL -
                                    E * S * mu_LE,
                             D(M) ~ S * delta_LM - mu_M * M,
                             D(P) ~ P^2 * rho_P - P * mu_P - E * P * mu_PE -
                                    S * P * mu_PL,
                         ], t, states, parameters)
measured_quantities = [y1 ~ N, y2 ~ E, y3 ~ S + M, y4 ~ P]

ic = [1.0, 1.0, 1.0, 1.0, 1.0]
time_interval = [0.0, 1.0]
datasize = 10
p_true = [1, 1.3, 1.1, 1.2, 1.1, 1, 0.5, 1.0, 1.0, 1.0, 1.0, 0.9, 1.2] # True Parameters
data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval,
                                              p_true, ic, datasize; solver = solver)

using Sundials
data_sample2 = ParameterEstimation.sample_data(model, measured_quantities, time_interval,
                                              p_true, ic, datasize; solver = CVODE_BDF())

@show data_sample[N] - data_sample2[N]
@show data_sample[E] - data_sample2[E]
@show data_sample[M + S] - data_sample2[M + S]
@show data_sample[P] - data_sample2[P]
@show data_sample["t"] - data_sample2["t"]

data_sample[N] - data_sample2[N] = [0.0, 5.041704831398874e-10, 2.8465840795632857e-10, 5.82355275113855e-11, -2.0006651890724925e-10, -2.724666603448611e-10, -3.192306419208535e-10, -3.31299987443856e-10, -5.182496654043689e-10, -5.951490689160011e-10]

data_sample[E] - data_sample2[E] = [0.0, -8.092885250832182e-10, -3.5877745219181634e-10, 1.302885577203483e-10, 6.291624998766565e-10, 7.014300251739769e-10, 6.842985067478935e-10, 6.195496893290908e-10, 8.133770323937028e-10, 7.955597847164597e-10]

data_sample[M + S] - data_sample2[M + S] = [0.0, 1.227966395234148e-9, 7.136780055816416e-10, 1.9830870279236024e-10, -4.093141381389387e-10, -6.53945786410759e-10, -8.094991343909896e-10, -9.063720884938675e-10, -1.2748408995122418e-9, -1.4841798900988579e-9]

data_sample[P] - data_sample2[P] = [0.0, 2.7097513122242844e-10, 1.0363487845665986e-10, -4.4600545479056564e-11, -1.900888335626405e-10, -2.2204182936746975e-10, -2.3077273425542444e-10, -2.208757066135547e-10, -2.9855290462066364e-10, -3.1491209639966655e-10]

data_sample["t"] - data_sample2["t"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]