[R-sig-dyn-mod] stiff ode - accuracy?

[Tim Keitt]

On Fri, Apr 29, 2016 at 9:50 AM, Thomas Petzoldt <
thomas.petzoldt at tu-dresden.de> wrote:

> so I
> wonder if and how other "different software" handles this.
>

I get the same result with odeintr:

library(odeintr)

sys = '
dxdt[0] = k * x[0] * x[1] / (K + x[1]);
dxdt[1] = -0.75 * k * x[0] * x[1] / (K + x[1]);
'

jac = '
J(0, 0) = k * x[1]/(K + x[1]);
J(0, 1) = k * x[0]/(K + x[1]) - k * x[0] * x[1]/((K + x[1]) * (K + x[1]));
J(1, 0) = -0.75 * k * x[1]/(K + x[1]);
J(1, 1) = -0.75 * k * x[0]/(K + x[1]) - -0.75 * k * x[0] * x[1]/((K + x[1])
*(K + x[1]));
dfdt[0] = 0.0;
dfdt[1] = 0.0;
'

compile_implicit("stiff", sys, jacobian = jac,
                 pars = c(k = 0.3, K = 1e-6),
                 atol = 1e-6, rtol = 1e-6,
                 const = TRUE, rebuild = TRUE)

times = seq(from = 0, to = 20, by = 0.01)

ini = c(0.05, 5)

res = stiff_at(ini, times)

par(mfrow = c(2, 2))

plot(X1 ~ Time, res, type = "l", lwd = 2)
plot(X2 ~ Time, res, type = "l", lwd = 2)

compile_implicit("stiff", sys, jacobian = jac,
                 pars = c(k = 0.3, K = 1e-6),
                 atol = 1e-7, rtol = 1e-7,
                 const = TRUE, rebuild = TRUE)

res = stiff_at(ini, times)

plot(X1 ~ Time, res, type = "l", lwd = 2)
plot(X2 ~ Time, res, type = "l", lwd = 2)


http://www.keittlab.org/

	[[alternative HTML version deleted]]