I’ve been looking into using DifferentialEquations.jl to solve some IVPs (chemical kinetics problems more precisely, something I have extensive experience working with using the scientific python stack).
I was delighted to see there’s a DSL for chemical kinetics: DiffEqBiological
So I went ahead trying to model a problem of mine which has been demanding for Boost’s odeint, GNU Scientific Library, as well as Sundials’ CVode (unless tweaking and tuning tolerances & constraints extensively).
I find that the Rosenbrock23 algorithm as implemented in DifferentialEquations fared best. Here’s my problem definition:
using DifferentialEquations, DiffEqBiological, Plots
rn = @reaction_network begin k1*A^2, 2B + 2A ⇒ C + 2D k2, 2E --> C k3, 2F --> G k4*E*A, E + B + A ⇒ C + D k5, A + F --> D k6, E + F --> B k7, A + G --> F + D k8, A + H --> I k9*I*A, 2B + I + A ⇒ G + 2D k10, A + J --> K k11, E + G --> F + B k12, E + H --> J k13, E + J --> 2F k14, E + I --> K k15, F + G --> J + B k16, F + I --> H + D k17, F + J --> H + B k18, 2J --> G + H k19*I*J, B + J + I ⇒ G + H + D k20*G, 2G ⇒ H + 2B k21, F + K --> B + I k22, M + G --> D + J k23, M + K --> D + I k24, M + C --> D + E k25, B --> L + D k26, L + D --> B k27, G --> L + K k28, L + K --> G k29, G + D --> K + B k30, K + B --> G + D k31, F --> L + M k32, L + M --> F k33, F + D --> M + B k34, M + B --> F + D k35, J --> L + I k36, L + I --> J k37, J + D --> I + B k38, I + B --> J + D k39, E --> L + A k40, L + A --> E k41, E + D --> A + B k42, A + B --> E + D k43, E + B --> C + F k44, C + F --> E + B k45, M + H --> N k46, N --> M + H k47, S + D --> B + R k48, B + R --> S + D k49, R + F --> B + Q k50, R + M --> Q + D k51, E + R --> S + A k52, S + A --> E + R k53, G + Q --> J + R k54, K + Q --> R + I k55, J + Q --> R + H k56, E + Q --> R k57*H*T^2, 2B + 4T + H ⇒ 4L + 4V k58*E*T, E + 2T ⇒ P + F k59, U + F --> L + W k60, U + I --> V + H k61*U^2, B + 2U ⇒ 2L + V + W k62, V + F --> U + D k63, W + A --> X k64*X, B + X ⇒ U + 2D k65, S + X --> R + U + D k66, L + X --> U + D k67, X + H --> W + I k68, P + A --> O + M k69*A*V, B + V + A ⇒ T + 2D k70, Q + H --> B + T k71, U + M --> W k72, T + I --> W k73*T*U, B + T + U ⇒ 2L + 2V k74, 0 ⇒ A k75, 0 ⇒ G k76, 0 ⇒ C k77, 0 ⇒ F k78, 0 ⇒ E k79, 0 ⇒ L k80, B ⇒ 0 k81, 0 ⇒ H end k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21 k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36 k37 k38 k39 k40 k41 k42 k43 k44 k45 k46 k47 k48 k49 k50 k51 k52 k53 k54 k55 k56 k57 k58 k59 k60 k61 k62 k63 k64 k65 k66 k67 k68 k69 k70 k71 k72 k73 k74 k75 k76 k77 k78 k79 k80 k81 p = Dict([(:k1, 7.226e+09), (:k2, 5.126e+09), (:k3, 4.806e+09), (:k4, 2.754e+10), (:k5, 3.546e+10), (:k6, 1.092e+10), (:k7, 1.357e+10), (:k8, 2.285e+10), (:k9, 1.295e+10), (:k10, 1.295e+10), (:k11, 3.635e+07), (:k12, 1.304e+10), (:k13, 1.135e+10), (:k14, 1.135e+10), (:k15, 2.911e+07), (:k16, 1.095e+10), (:k17, 8.828e+09), (:k18, 8.365e+05), (:k19, 1.01e+08), (:k20, 1.281e-07), (:k21, 4.057e+09), (:k22, 4.057e+09), (:k23, 7.825e+08), (:k24, 1.276e+08), (:k25, 2.086e-05), (:k26, 1.174e+11), (:k27, 0.09327), (:k28, 5.007e+10), (:k29, 1.326e+10), (:k30, 1.264e+06), (:k31, 0.09327), (:k32, 5.007e+10), (:k33, 1.326e+10), (:k34, 1.264e+06), (:k35, 7.692e+05), (:k36, 5.007e+10), (:k37, 1.326e+10), (:k38, 0.1533), (:k39, 5.752), (:k40, 2.092e+10), (:k41, 2.421e+07), (:k42, 15.64), (:k43, 4.481e-05), (:k44, 3.937e+07), (:k45, 3.739e+09), (:k46, 2593), (:k47, 1.326e+10), (:k48, 2.407e+05), (:k49, 7.982e+07), (:k50, 7.982e+06), (:k51, 9.935e+05), (:k52, 6.659e+05), (:k53, 8.974e+07), (:k54, 8.974e+07), (:k55, 8.974e+07), (:k56, 1.994e+10), (:k57, 5.983e+06), (:k58, 1.097e+10), (:k59, 4.487e+09), (:k60, 9.972e+07), (:k61, 9.972e+07), (:k62, 9.672e+08), (:k63, 9.6e+09), (:k64, 9.972e+04), (:k65, 1.994e+08), (:k66, 1.994e+10), (:k67, 1.994e+08), (:k68, 9.6e+09), (:k69, 3.576e+09), (:k70, 1.097e+09), (:k71, 3.49e+09), (:k72, 6.681e+09), (:k73, 7.24e+06), (:k74, 4.266e-08), (:k75, 1.104e-08), (:k76, 6.793e-09), (:k77, 4.353e-08), (:k78, 9.373e-09), (:k79, 4.266e-08), (:k80, 6.561e-08), (:k81, 4.066e-13)]) ics = Dict([(:B, 55.4), (:R, 0.1), (:L, 1e-11), (:D, 0.001), (:A, 1e-24), (:C, 1e-24), (:E, 1e-24), (:F, 1e-24), (:G, 1e-24), (:H, 1e-24), (:I, 1e-24), (:J, 1e-24), (:K, 1e-24), (:M, 1e-24), (:N, 1e-24), (:O, 1e-24), (:P, 1e-24), (:Q, 1e-24), (:S, 1e-24), (:T, 1e-24), (:U, 1e-24), (:V, 1e-24), (:W, 1e-24), (:X, 1e-24)])
Solving this succeeds (albeit at somewhat high computational cost, and I’d actually like tspan to go up to 1e13 – something which currently fails):
tspan = (0., 3.6e7) u0 = Array([get(ics, k, 0.0) for k=keys(speciesmap(rn))]) parr = Array([p[k] for k=keys(paramsmap(rn))]) oprob = ODEProblem(rn, u0, tspan, parr) sol = solve(oprob, reltol=1e-9, abstol=1e-10, Rosenbrock23() ) sol.retcode
sol.destats I see ~100k rhs evaluations and ~30k jacobian creations (to be compared with ~30k rhs evals & ~3k jacobian evaluations when solving using my hand-crafted CVode solver).
Plotting the result gives a hint why the computational cost is high (high frequency oscillatory behavior in a region which should be smooth):
gr() lgt = log10.(clamp.(sol.t[2:end], 0, Inf)) lgu = log10.(clamp.(hcat(sol.u[2:end]...)', 0, Inf)) plot(lgt, lgu, fmt = :png, legend=:topleft, legendfontsize=8)
lowering the tolerances seem to make the solution deteriorate.
So my question is: what other options do DifferentialEquations/DiffEqBiological offer for tuning the solution behavior (e.g. are positivity constraints implied from using DiffEqBiological, or does one need to pass that infomation on when solving)?