Different ODE solutions by integrate_ode_bdf and integrate_ode_rk45

Modeling
1

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
}

Rplot
Rplot657Ă—534 54 KB

3 Likes
2

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

3

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
4
  • 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
5

Hi,

It is not uncommon to observe different results using different solvers. One, sometimes both, stray away from the true answer.

Rescaling and adjusting tolerances address the same problem. When is a state variable zero (atol), and when is a change in that state variable not changing between steps (rtol).

I would recommend always rescaling, and adjust tolerances thereafter. You pay the price of very stringent tolerances with integration time. In my own work I always use a stiff solver (bdf). The bdf solver’s accuracy outweighs it’s computational costs.

Good luck

1 Like