Different ODE solutions by integrate_ode_bdf and integrate_ode_rk45

I am working with an ODE system, where quantities live on the nanoscale. Rescaling of the states and parameter has helped to resolve the issue, but I would like to double-check my understanding: while working on the original scale, I was obtaining extremely different solutions produced by integrate_ode_rk45 and integrate_ode_bdf. I assume that it is due to the fact that in the rdk45 case the numerical error is larger than the modelled quantities. Is that true? The relevant code and resulting plots are below. Any insights and comments are welcome.


R code:

n0 <- rep(0, 7)
n0[1] <- 10^3 *10^(-9)
n0[2] <- 10^4 *10^(-9)
n0[3] <- 250  *10^(-9)
n0[4] <- 0
n0[5] <- 0
n0[6] <- 0
n0[7] <- 0

# parameters
kp1 <- 500  *10^(-9)
kp2 <- 5000 *10^(-9)
km1 <- 10^(-1)
km2 <- 10^(-1)
kp3 <- 10^(2)
kdeg <- 10^(-1)


#times
t	<- seq(0, 60, length.out = 1000)
t0 <- t[1] 
t <- t[-1]
N_times <- length(t);

data <-  list (N_times  = N_times, n0 = n0, t0 = t0, t = t,
               kp1=kp1, kp2=kp2, km1=km1, km2=km2, kp3=kp3, kdeg=kdeg)

fit <- stan(file = 'ode_generate_bdf.stan', 
            data = data, 
            algorithm="Fixed_param",
            chains = 1,
            iter =1)

s <- rstan::extract(fit)
str(s$n)
n <- data.frame(s$n[1,1:N_times,1:7])


par(mfrow=c(2,2))

plot(1, xlim=c(min(t), max(t)), ylim=c(min(n[,c(1)]),max(n[,c(1)])), type='n', main='n1, integrate_ode_bdf' )
lines(t, n[,1], col="red")

plot(1, xlim=c(min(t), max(t)), ylim=c(min(n[,c(7)]),max(n[,c(7)])), type='n', main='n7, integrate_ode_bdf' )
lines(t, n[,7], col="blue")


fit <- stan(file = 'ode_generate_rdk45.stan', 
            data = data, 
            algorithm="Fixed_param",
            chains = 1,
            iter =1)

s <- rstan::extract(fit)
str(s$n)
n <- data.frame(s$n[1,1:N_times,1:7])
head(n)

plot(1, xlim=c(min(t), max(t)), ylim=c(min(n[,c(1)]),max(n[,c(1)])), type='n', main='n1, integrate_ode_rdk45' )
lines(t, n[,1], col="red")

plot(1, xlim=c(min(t), max(t)), ylim=c(min(n[,c(7)]),max(n[,c(7)])), type='n', main='n7, integrate_ode_rdk45' )
lines(t, n[,7], col="blue")

And Stan code:

functions {
  real[] ode_protac(real t, real[] n, real[] theta, 
             real[] x_r, int[] x_i) {
      
      real kp1 = theta[1];
      real km1 = theta[2];
      real kp2 = theta[3];
      real km2 = theta[4];
      real kp3 = theta[5];
      real kdeg = theta[6];
      
      real dndt[7];
            
      dndt[1] =  - kp1 * n[1] * n[3] + km1 *n[4] - kp1 * n[1] * n[5] + km1 * n[6];
      dndt[2] =  - kp2 * n[2] * n[3] + km2 *n[5] - kp2 * n[2] * n[4] + km2 * n[6];
      dndt[3] =  - kp1 * n[1] * n[3] + km1 *n[4] - kp2 * n[2] * n[3] + km2 * n[5];
      dndt[4] =    kp1 * n[1] * n[3] - km1 *n[4] - kp2 * n[2] * n[4] + km2 * n[6];
      dndt[5] =    kp2 * n[2] * n[3] - km2 *n[5] - kp1 * n[1] * n[5] + km1 * n[6] + kp3 * n[6];
      dndt[6] =    kp2 * n[2] * n[4] - km2 *n[6] + kp1 * n[1] * n[5] - km1 * n[6] - kp3 * n[6];
      dndt[7] =    kp3 * n[6] - kdeg*n[7];  
      return dndt;
  }
}
data {
  int<lower=1> N_times;
  real n0[7];
  real t0;
  real t[N_times];
  real kp1;
  real kp2;
  real km1;
  real km2;
  real kp3;
  real kdeg;
}
transformed data {
}
parameters {
}
transformed parameters{
}
model {
}
generated quantities {
  real n[N_times, 7];
  real theta[6];
  real x_r[0];
  int x_i[0];
  
  theta[1] = kp1;
  theta[2] = km1;
  theta[3]=kp2;
  theta[4]=km2;
  theta[5]=kp3;
  theta[6]=kdeg;
    
  n = integrate_ode_bdf(ode_protac, n0, t0, t, theta, x_r, x_i); //![Rplot|615x500](upload://1wYcg1HSPf1J58bNfjMorvGF0DJ.png) or integrate_ode_rk45
}

3 Likes

Tagging @wds15 and @bbbales2 , two of our ODE aficionados.

The ode integrators have by default only a absolute accuracy of about 1e-8. Rescaling is the right thing to do and I would anyway recommend to specify the tolerances manually.

3 Likes
  • Why the bdf solution spits garbage(nonsmoothness):
    By the 1st equation your unknown is at scale 10^20, set the BDF solver’s tolerance & max_step at
integrate_ode_bdf(..., 1.e-19, 1.e-20, 1000000);

should do the job.

  • Why the rk45 solution spits garbage(oscillations):
    The issue in n[7] is actually caused by n[6] because the 7th equation is stable in terms of n[7]. The 6th equation gives n[6] two distinct reaction scales: that by kp1 & kp2 at scale 10^-6, and that by km1 & km2 at scale 10^-1. Rk45 solver is unfit for this stiff problem. But since we are doing exploration, my guess if you decrease tolerance to something like
integrate_ode_rk45(..., 1.e-30, 1.e-35, 10000000);

you should be able to recover BDF solver’s result, albeit inefficiently.

5 Likes