System of ODE with variable parameter

[efus]

Hello everyone, new R studio user here.

I need a little bit of help to understand how to do what I would like to do and I'm quite sure R can do.

To make things simple, I have a system of ODE with one parameter that is calculated from the result of the same ODE with an algebraic equation.

Now, If I put this therm explicit in the ODE system, no problem, all work.
What I would like to do is to put the equation in the ODE as a variable, calculating it in a separate place at every iteration.

Probably code can explain things better than how I'm able to do.
I also add a screen from Smath that is quite clear regarding this.

Any advice is welcome!

[image]

In R this work

> 
> ## LSODA
> 
> SPCmod <- function(t, y, parms) {
>   with(as.list(c(parms, y)), {
>     CO2.i<-y[1]
>     O2.i<-y[2]       
>     dCO2.i <- (((A.p * p.CO2 * P.atm / L.f) * (0.01 * (CO2.0 - CO2.i)) + (W.s * ((V.mCO2*O2.i) / (K.mCO2 + (1 + CO2.i / K.iCO2) * O2.i)))) * 100 / V.f)  #ODE atmosphere evolution O2           
>     dO2.i <- (((A.p * p.O2 * P.atm / L.f) * (0.01 * (O2.0 - O2.i)) - (W.s * ((V.mO2*O2.i) / (K.mO2 + (1 + CO2.i / K.iO2) * O2.i)))) * 100 / V.f )        #ODE atmosphere evolution O2 
>      res <- c(dCO2.i, dO2.i)
>     
>     list(res)
>   })
> }
> 
> ## Parameters 
> 
> parms  <- c(A.p = 0.069, CO2.0 = 0.03, O2.0 = 21, V.f = 492, V.mO2 = 22.71, V.mCO2 = 17.64, 
>             W.s = 0.200,
>             K.iO2 = 14.42, K.iCO2 = 11.99, K.mO2 = 7.63, K.mCO2 = 5.08, P.atm= 1, 
>             p.O2= 5.6*10^(-3) ,
>             p.CO2= 20.9*10^(-3)  , L.f = 45*10^(-6))
> 
> ## vector of timesteps
> 
> times  <- seq(1, 5000, length=50)
> 
> 
> ## Start values for steady state
> y <- xstart <- c(dCO2.i = 0.03, dO2.i = 21)
> 
> ## Solving
> out <-  lsoda(xstart, times, SPCmod, parms) 
> 
> ## Plotting
> mf <- par("mfrow")
> plot(out, main = c("CO2", "O2"))
> par(mfrow = mf)
> tail(out)

This don't

> SPCmod <- function(t, y, parms) {
>   with(as.list(c(parms, y)), {
>     CO2.i<-y[1]
>     O2.i<-y[2]       
>     
>     
>     dCO2.i <- (((A.p * p.CO2 * P.atm / L.f) * (0.01 * (CO2.0 - CO2.i)) + (W.s * tr.CO2)) * 100 / V.f)  #ODE atmosphere evolution O2           
>     dO2.i <- (((A.p * p.O2 * P.atm / L.f) * (0.01 * (O2.0 - O2.i)) - (W.s * tr.O2)) * 100 / V.f )       #ODE atmosphere evolution O2 
>     
>     res <- c(dCO2.i, dO2.i)
>     
>     list(res)
>   })
> }
> 
> ## Parameters 
> 
> parms  <- c(A.p = 0.069, CO2.0 = 0.03, O2.0 = 21, V.f = 492, V.mO2 = 22.71, V.mCO2 = 17.64, 
>             W.s = 0.200,
>             K.iO2 = 14.42, K.iCO2 = 11.99, K.mO2 = 7.63, K.mCO2 = 5.08, P.atm= 1, 
>             p.O2= 5.6*10^(-3) ,
>             p.CO2= 20.9*10^(-3)  , L.f = 45*10^(-6), tr.CO2 = (V.mCO2*O2.i) / ((K.mCO2 + (1 + CO2.i / K.iCO2) * O2.i)) 
>             ,tr.O2 = (V.mO2*O2.i) / ((K.mO2 + (1 + CO2.i / K.iO2) * O2.i))
>  )
> 
> 
> ## vector of timesteps
> 
> times  <- seq(1, 5000, length=50)
> 
> 
> ## Start values for steady state
> y <- xstart <- c(dCO2.i = 0.03, dO2.i = 21)
> 
> ## Solving
> out <-  lsoda(xstart, times, SPCmod, parms) 
> 
> ## Plotting
> mf <- par("mfrow")
> plot(out, main = c("CO2", "O2"))
> par(mfrow = mf)
> tail(out)

[efus]

Ok...playing a little I have found a way to make it work.

Really don't know why, so, again, any advice is still welcome.

## LSODA

SPCmod <- function(t, y, parms) {
  with(as.list(c(parms, y)), {
    CO2.i<-y[1]
    O2.i<-y[2]       
    
    tr.CO2 = (V.mCO2*O2.i) / ((K.mCO2 + (1 + CO2.i / K.iCO2) * O2.i)) 
    
    tr.O2 = (V.mO2*O2.i) / ((K.mO2 + (1 + CO2.i / K.iO2) * O2.i))
    
    dCO2.i <- (((A.p * p.CO2 * P.atm / L.f) * (0.01 * (CO2.0 - CO2.i)) + (W.s * tr.CO2)) * 100 / V.f)  #ODE atmosphere evolution O2           
    dO2.i <- (((A.p * p.O2 * P.atm / L.f) * (0.01 * (O2.0 - O2.i)) - (W.s * tr.O2)) * 100 / V.f )       #ODE atmosphere evolution O2 
    
    res <- c(dCO2.i, dO2.i)
    
    list(res)
  })
}

## Parameters 

parms  <- c(A.p = 0.069, CO2.0 = 0.03, O2.0 = 21, V.f = 492, V.mO2 = 22.71, V.mCO2 = 17.64, 
            W.s = 0.200,
            K.iO2 = 14.42, K.iCO2 = 11.99, K.mO2 = 7.63, K.mCO2 = 5.08, P.atm= 1, 
            p.O2= 5.6*10^(-3) ,
            p.CO2= 20.9*10^(-3)  , L.f = 45*10^(-6)
 )


## vector of timesteps

times  <- seq(1, 5000, length=50)


## Start values for steady state
y <- xstart <- c(dCO2.i = 0.03, dO2.i = 21)

## Solving
out <-  lsoda(xstart, times, SPCmod, parms) 

## Plotting
mf <- par("mfrow")
plot(out, main = c("CO2", "O2"))
par(mfrow = mf)
tail(out)