Different ODE solutions by integrate_ode_bdf and integrate_ode_rk45

Modeling
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