Differential Eq using deSolve in R

Question

My apologies, for being unclear earlier. I now understand the function a bit more, but could use some assistance on a few aspects.

I would like to get back a relationship of conversion ( X ) versus volume ( V ), or the other way around would be fine as well. It would seem to me that the traditional "times" term is what I want to replace with an X sequence from 0 - 1, X is conversion remember so bounded by 0 and 1.0

Below, rw is the reaction rate, and is a function of the partial pressures at any given moment, which are described as P.w, P.x, P.y, and P.z which themselves are functions of the initial conditions (P.w0, v.0) and the conversion, again X.

Thank you in advance

rm(list = ls())

weight <- function( Vols, State, Pars ) {


  with(as.list(c(State, Pars)), {

    y = 1
    delta = 2
    ya.0 = 0.4
    eps = ya.0 * delta

    temp = 800
    R = 8.314

    k.2 = exp( (35000 / ( R*temp )) - 7.912 )

    K.3 = exp( 4.084  /     temp    - 4.33  )

    P.w <- P.w0 * ( 1 -   X ) * y / ( 1 + eps * X )
    P.x <- P.w0 * ( 1 - 2*X ) * y / ( 1 + eps * X )
    P.y <- P.w0 * ( 1 +   X ) * y / ( 1 + eps * X )
    P.z <- P.w0 * ( 1 + 4*X ) * y / ( 1 + eps * X )

    r.w <- k.2 * ( K.3 * P.w * P.x ^ 2 - P.y * P.z^4 )

    F.w0 <- P.w0 * v.0 / ( R * temp )

    dX.dq <- r.w / F.w0
    res <- dX.dq

    return(list(res))

  })
}

pars <- c( y = 1, 
           P.w0 = 23, 
           v.0 = 120 )

yini <- c( X = 0 )

vols <- seq( 0 , 100 , by = 1 )

out <- ode( yini , vols , weight , pars )

I don't know wher eyour problem comes from, but please don't use `T` as a variable name. `T` is a boolean in R. — Heroka, Feb 28 '16 at 20:01
Your input of `V` in `ode` is `0`, hence it is divided by 0 in the `vol.func` function? — Sixiang.Hu, Feb 28 '16 at 20:59

score 0 · Answer 1 · answered Feb 28 '16 at 21:00

Just running

vol.func(0,0,params)

i.e., evaluating the gradient at the initial conditions, gives NaN. The proper way to diagnose this is to divide your complex gradient expressions up into separate terms and see which one is causing trouble. I'm not going to go through this in detail, but as @Sixiang.Hu points out in comments above, you're dividing by V in your gradient function, which will cause infinite values if the numerator is finite or NaN values if the numerator is zero ...

More generally, it's not clear whether you understand that the first argument to the gradient function (your vol.func) is supposed to be the current time, not a value of the state variable. Perhaps V is supposed to be your state variable, and X should be a parameter ...?

Differential Eq using deSolve in R

1 Answers1