`deSolve`

pitfallBeware!

As suggested in the `ode()`

examples, I typically use a `with()`

statement in my gradient function for convenience, so I don't have to refer to the order of the parameters in the parameter vector or the order of the variables in the state vector. I thought this meant that I didn't have to worry about ordering at all, but that's **not true**. Here's an example/explanation.

Set up a simple linear model (\( x \) exponentially increasing, \( y \) exponentially decreasing):

```
library("deSolve")
gradfun <- function(t, statevec, params) {
g <- with(as.list(c(statevec, params)), c(x = a * x, y = b * y))
list(g)
}
params0 <- c(a = log(2), b = -log(2))
tvec <- 0:5
```

Starting at `x=0`

, `y=1`

, the right way:

```
ode(c(x = 0, y = 1), tvec, gradfun, params0)
```

```
## time x y
## 1 0 0 1.00000
## 2 1 0 0.50000
## 3 2 0 0.25000
## 4 3 0 0.12500
## 5 4 0 0.06250
## 6 5 0 0.03125
```

As should happen, \( x \) remains at zero and \( y \) decreases exponentially back toward zero.

What happens if we get the order of the initial conditions wrong (even though the names are right?)

```
ode(c(y = 1, x = 0), tvec, gradfun, params0)
```

```
## time y x
## 1 0 1.0000 0.0000
## 2 1 0.7692 -0.6390
## 3 2 0.1835 -0.9830
## 4 3 -0.4870 -0.8734
## 5 4 -0.9327 -0.3607
## 6 5 -0.9479 0.3185
```

The system is now *cross-wired*, equivalent to

\[ \begin{split} \dot x & = -y \\ \dot y & = x \end{split} \]

```
cbind(cos(log(2) * tvec), -sin(log(2) * tvec))
```

```
## [,1] [,2]
## [1,] 1.0000 0.0000
## [2,] 0.7692 -0.6390
## [3,] 0.1835 -0.9830
## [4,] -0.4870 -0.8734
## [5,] -0.9327 -0.3607
## [6,] -0.9479 0.3185
```

What is happening here?

- The gradient gets computed according to the
*names*of the vector; because of the way we set the initial values vector, \( y \) is the first element and \( x \) is the second element. - all other internal stuff deals with the vectors by position (i.e. ignoring names), so the first computed gradient, which is supposed to be applied to the \( x \) component, gets applied to the first element (which is the \( y \) component) instead, so we end up with the cross-wired solution.

This is a workable solution, I think: it depends on the user having assigned appropriate names to the gradient vector (this is good practice in any case).

```
myCheckFunc <- deparse(deSolve:::checkFunc)
myCheckFunc <- c(myCheckFunc[1:3], " if (!all(names(tmp[[1]])==names(y))) {",
"warning('name mismatch between y and results of grad()')", " }", myCheckFunc[4:length(myCheckFunc)])
myCheckFunc <- eval(parse(text = myCheckFunc), envir = environment(deSolve:::checkFunc))
assignInNamespace("checkFunc", myCheckFunc, ns = "deSolve")
```

```
ode(c(y = 1, x = 0), tvec, gradfun, params0)
```

```
## Warning: name mismatch between y and results of grad()
```

```
## time y x
## 1 0 1.0000 0.0000
## 2 1 0.7692 -0.6390
## 3 2 0.1835 -0.9830
## 4 3 -0.4870 -0.8734
## 5 4 -0.9327 -0.3607
## 6 5 -0.9479 0.3185
```

No warning if the gradient vector is *not* named (it would be possible to add this check, i.e. warn if `y`

has names but the result of `grad()`

doesn't).

```
gradfun2 <- function(t, statevec, params) {
g <- with(as.list(c(statevec, params)), c(a * x, b * y))
list(g)
}
ode(c(y = 1, x = 0), tvec, gradfun2, params0)
```

```
## time y x
## 1 0 1.0000 0.0000
## 2 1 0.7692 -0.6390
## 3 2 0.1835 -0.9830
## 4 3 -0.4870 -0.8734
## 5 4 -0.9327 -0.3607
## 6 5 -0.9479 0.3185
```