
rxode2 Jacobian specification and Stiff Systems
2023-08-22
Source:vignettes/rxode2-stiff.Rmd
rxode2-stiff.Rmd
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 rxode2 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 rxode2 these can be specified by the following:
## rxode2 2.0.13.9000 using 1 threads (see ?getRxThreads)
## no cache: create with `rxCreateCache()`
Vtpol2 <- rxode2({
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
})
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
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 rxode2 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
## # i 194 more rows
While this is not stiff at mu=1, mu=1000 is a stiff system
## -- Solved rxode2 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
## # i 194 more rows
While this is easy enough to do, it is a bit tedious. If you have rxode2 setup appropriately, you can use the computer algebra system sympy to calculate the Jacobian automatically.
This is done by the rxode2 option calcJac
option:
Vtpol <- rxode2({
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);