# System of ODE with variable parameter

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.

[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)

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)

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.