RxODE Jacobian specification and Stiff Systems

2019-10-16

Stiff ODEs with Jacobian Specification

Occasionally, you may come across a stiff differential equation, that is a differential equation that is numerically unstable and small variations in parameters cause different solutions to the ODEs. One way to tackle this is to choose a stiff-solver, or hybrid stiff solver (like the default LSODA). Typically this is enough. However exact Jacobian solutions may increase the stability of the ODE. (Note the Jacobian is the derivative of the ODE specification with respect to each variable). In RxODE you can specify the Jacobian with the df(state)/dy(variable)= statement. A classic ODE that has stiff properties under various conditions is the Van der Pol differential equations.

In RxODE these can be specified by the following:

library(RxODE)

Vtpol2 <- RxODE({
    d/dt(y)  = dy
    d/dt(dy) = mu*(1-y^2)*dy - y
    ## Jacobian
    df(y)/dy(dy)  = 1
    df(dy)/dy(y)  = -2*dy*mu*y - 1
    df(dy)/dy(dy) = mu*(1-y^2)
    ## Initial conditions
    y(0) = 2
    dy(0) = 0
    ## mu
    mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
})

et <- eventTable();
et$add.sampling(seq(0, 10, length.out=200));
et$add.dosing(20, start.time=0);

s1 <- Vtpol2 %>%  solve(et, method="lsoda")
print(s1)

## ___________________________ Solved RxODE object ___________________________
## -- Parameters ($params): --------------------------------------------------
##  
## mu 
##  1 
## -- Initial Conditions ($inits): -------------------------------------------
##  y dy 
##  2  0 
## -- First part of data (object): -------------------------------------------
## # A tibble: 200 x 3
##     time     y      dy
##    <dbl> <dbl>   <dbl>
## 1 0       22    0     
## 2 0.0503  22.0 -0.0456
## 3 0.101   22.0 -0.0456
## 4 0.151   22.0 -0.0456
## 5 0.201   22.0 -0.0456
## 6 0.251   22.0 -0.0456
## # ... with 194 more rows
## ___________________________________________________________________________

While this is not stiff at mu=1, mu=1000 is a stiff system

s2 <- Vtpol2 %>%  solve(c(mu=1000), et)
print(s2)

## ___________________________ Solved RxODE object ___________________________
## -- Parameters ($params): --------------------------------------------------
##  
##   mu 
## 1000 
## -- Initial Conditions ($inits): -------------------------------------------
##  y dy 
##  2  0 
## -- First part of data (object): -------------------------------------------
## # A tibble: 200 x 3
##     time     y         dy
##    <dbl> <dbl>      <dbl>
## 1 0       22    0        
## 2 0.0503  22.0 -0.0000455
## 3 0.101   22.0 -0.0000455
## 4 0.151   22.0 -0.0000455
## 5 0.201   22.0 -0.0000455
## 6 0.251   22.0 -0.0000455
## # ... with 194 more rows
## ___________________________________________________________________________

While this is easy enough to do, it is a bit tedious. If you have RxODE setup appropriately, that is you have:

• Python installed in your system
• sympy installed in your system
• SnakeCharmR installed in R

You can use the computer algebra system sympy to calculate the Jacobian automatically.

This is done by the RxODE option calcJac option:

Vtpol <- RxODE({
    d/dt(y)  = dy
    d/dt(dy) = mu*(1-y^2)*dy - y
    ## Initial conditions
    y(0) = 2
    dy(0) = 0
    ## mu
    mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
}, calcJac=TRUE)

To see the generated model, you can use rxCat():

> rxCat(Vtpol)
d/dt(y)=dy;
d/dt(dy)=mu*(1-y^2)*dy-y;
y(0)=2;
dy(0)=0;
mu=1;
df(y)/dy(y)=0;
df(dy)/dy(y)=-2*dy*mu*y-1;
df(y)/dy(dy)=1;
df(dy)/dy(dy)=mu*(-Rx_pow_di(y,2)+1);