mjg211/phaser

Development version of phaseR, an R package for phase plane analysis of one- and two-dimensional autonomous ODE systems

phaseR

Phase plane analysis of one- and two-dimensional autonomous ODE systems

Description

phaseR provides functions to perform a qualitative analysis of one- and two-dimensional autonomous ordinary differential equation (ODE) systems, using phase plane methods. Programs are available to identify and classify equilibrium points, plot the direction field, and plot trajectories for multiple initial conditions. In the one-dimensional case, a program is also available to plot the phase portrait. Whilst in the two-dimensional case, programs are additionally available to plot nullclines and stable/unstable manifolds of saddle points. Many example systems are provided for the user.

Getting started

You can install the released version of phaseR from CRAN with:

install.packages("phaseR")

Alternatively, the latest development version available from GitHub can be installed with:

devtools::install_github("mjg211/phaseR")

An introductory example of how to make use of the package’s core functionality can be found below. More detailed support is available in the package vignette, which can be accessed with vignette("introduction", package = "phaseR"). For further help, please contact Michael Grayling at [email protected].

Example

As a basic example, we consider analysing the non-linear two-dimensional system of ODEs provided in phaseR via example12(). By hand, we typically first locate the nullclines and then identify the equilibrium points. Following this, we produce a plot from which trajectories can be sketched. This can all be seamlessly carried out in phaseR with:

example12_flowField   <- flowField(example12,
                                   xlim = c(-4, 4),
                                   ylim = c(-4, 4),
                                   add  = FALSE)
example12_nullclines  <- nullclines(example12,
                                    xlim   = c(-4, 4), 
                                    ylim   = c(-4, 4),
                                    points = 500)
y0                    <- matrix(c( 2,  2,
                                  -3,  0,
                                   0,  2,
                                   0, -3), 
                                nrow  = 4,
                                ncol  = 2,
                                byrow = TRUE)
example12_trajectory  <- trajectory(example12,
                                    y0   = y0,
                                    tlim = c(0, 10))
#> Note: col has been reset as required

It appears that both of the equilibria are unstable. We could verify this by hand, but we can also perform this analysis in phaseR using stability():

example12_stability_1 <- stability(example12,
                                   ystar = c(1, 1))
#> tr = 3, Delta = 4, discriminant = -7, classification = Unstable focus
example12_stability_2 <- stability(example12,
                                   ystar = c(-1, -1),
                                   h     = 1e-8)
#> tr = -1, Delta = -4, discriminant = 17, classification = Saddle

References

Grayling MJ (2014) phaseR: An R package for phase plane analysis of autonomous ODE systems. The R Journal 6(2):43-51. DOI: 10.32614/RJ-2014-023.

phaser's Issues

Package is removed from CRAN?

Running into the problem that students can't use phaseR anymore.

https://cran.r-project.org/web/packages/phaseR/index.html

Archived on 2019-09-08 as check problems were not corrected in time.

Time to get that .travis.yml to work ;)

vignette("phaseR") does not work

Just a (self?) reminder that vignette("phaseR"), as suggested in ?phaseR, does not work.
The call that work is vignette("my-vignette", package = "phaseR")

Missing value where TRUE/FALSE needed

I am trying to draw the phase diagram of the Solow economic growth model. Using y^a with a that is a non-natural number, I get:

Error in if (any(dx[i, j] != 0, dy[i, j] != 0)) { : 
  missing value where TRUE/FALSE needed

The code I use is:

solow <- function(t, y, parameters) {
  list(y^0.5)
}
example2_flowField     <- flowField(solow,
                                    xlim   = c(0, 4),
                                    ylim   = c(0, 4),
                                    system = "one.dim",
                                    add    = FALSE,
                                    xlab   = "t")

Can somebody help me?

Trajectory 1D system starts at t = 0, instead of (t, y) clicked

Came across this when using the phasePlaneAnalysis on a 1D function. It's not really a problem, but it might not be what you expect to happen.

Error when writing equations in matrix form

Thanks for developing and maintaining the package!

I run into the following error when trying to plot a phase diagram on a 2-dim system. The function is:

## Pollution model:
pollution <- function(t, y, params){
  with(as.list(c(y, params)), {
    x = numeric(length = n)
    x = y
    pollutant <- u - s*x + v * (x^alpha/(z^alpha + x^alpha))
    outflow <-  (A_ij * delta_ij) %*% x
    inflow <-  t(A_ij * delta_ij) %*% x

    dx <- pollutant + (inflow - outflow)
    
    return(list(c(dx)))
  })
}

With the following parameters:

n <- 1 # number of systems
delta_ij <- matrix(runif(n^2, min = 0.02, max = 0.05), ncol = n)
A_ij <- matrix(rbinom(n^2, 1, prob = 0.3), n, n)

## Parameters: 
params <- list(
    n = n,                             # number of systems (e.g. lakes)
    u = rep(0.5, n),             # pollution load from humans
    s = rep(2.5, n),             # internal loss rate (sedimentation)
    v = rep(10, n),              # max level of internal nutrient release
    z = rep(2.2, n),             # threshold
    alpha = 4 ,                    # sharpness of the shift
    delta_ij = delta_ij,        # matrix of difussion terms
    A_ij = A_ij                     # adjacency matrix
)

phaseR works on the one-dimensional case, eg:

fld <- flowField(pollution, xlim = c(0,5), ylim = c(-2,7), parameters = params,
          system = "one.dim", add = FALSE)
nll <- nullclines(pollution, xlim = c(0,5), ylim = c(-2,7), parameters = params,
          system = "one.dim", add = TRUE)
trj <- trajectory(pollution, y0 = c(-0.5, 1, 3, 5), tlim = c(0,5), parameters = params,
          system = "one.dim", add = TRUE)

But as soon as I change n = 2 I get the following error:

fld <- flowField(pollution, xlim = c(0,5), ylim = c(-2,7), system = "two.dim", add = FALSE, parameters = params)

Error in (A_ij * delta_ij) %*% x : non-conformable arguments

I don't quite understand the origin of the error. The system works with n = 1, 2, or 100s when doing numerical simulations with deSolve, but for some reason fails on the two-dimensional case in phaseR.

If I rewrite the function in a non-matrix form such as:

pollution2 <- function(t,y,params){
    with(as.list(c(y, params)), {
    x = numeric(length = n)
    x = y
    dx <- u[1] - s[1]*x[1] + v[1] * (x[1]^alpha/(z[1]^alpha + x[1]^alpha)) + (A_ij[2,1] * delta_ij[2,1]) %*% x[1] - (A_ij[1,2] * delta_ij[1,2]) %*% x[2]
    
    dy <- u[2] - s[2]*x[2] + v[2] * (x[2]^alpha/(z[2]^alpha + x[2]^alpha)) + (A_ij[1,2] * delta_ij[1,2]) %*% x[2] - (A_ij[1,2] * delta_ij[1,2]) %*% x[1]
    
    return(list(c(dx, dy)))
  })
}

Trying to get the flow field throughs a different error:

fld <- flowField(pollution2, xlim = c(0,5), ylim = c(-2,7), system = "two.dim", add = FALSE, parameters = params)

Error in if (any(dx[i, j] != 0, dy[i, j] != 0)) { : missing value where TRUE/FALSE needed

Which seems similar to issue #13
Any help is greatly appreciated!

Below my session info if useful:

session_info()
─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.2.2 (2022-10-31)
 os       macOS Ventura 13.5.2
 system   x86_64, darwin17.0
 ui       RStudio
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Stockholm
 date     2023-09-23
 rstudio  2023.06.0+421 Mountain Hydrangea (desktop)
 pandoc   2.14.2 @ /usr/local/bin/pandoc

─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────
 package        * version date (UTC) lib source
 assertthat       0.2.1   2019-03-21 [1] CRAN (R 4.2.0)
 backports        1.4.1   2021-12-13 [1] CRAN (R 4.2.0)
 broom            1.0.2   2022-12-15 [1] CRAN (R 4.2.2)
 cachem           1.0.7   2023-02-24 [1] CRAN (R 4.2.0)
 callr            3.7.3   2022-11-02 [1] CRAN (R 4.2.0)
 cellranger       1.1.0   2016-07-27 [1] CRAN (R 4.2.0)
 cli              3.6.1   2023-03-23 [1] CRAN (R 4.2.0)
 coda             0.19-4  2020-09-30 [1] CRAN (R 4.2.0)
 colorspace       2.1-0   2023-01-23 [1] CRAN (R 4.2.0)
 crayon           1.5.2   2022-09-29 [1] CRAN (R 4.2.0)
 DBI              1.1.3   2022-06-18 [1] CRAN (R 4.2.0)
 dbplyr           2.2.1   2022-06-27 [1] CRAN (R 4.2.0)
 deSolve        * 1.34    2022-10-22 [1] CRAN (R 4.2.0)
 devtools       * 2.4.5   2022-10-11 [1] CRAN (R 4.2.0)
 digest           0.6.33  2023-07-07 [1] CRAN (R 4.2.0)
 dplyr          * 1.1.1   2023-03-22 [1] CRAN (R 4.2.0)
 ellipsis         0.3.2   2021-04-29 [1] CRAN (R 4.2.0)
 fansi            1.0.4   2023-01-22 [1] CRAN (R 4.2.0)
 fastmap          1.1.1   2023-02-24 [1] CRAN (R 4.2.0)
 forcats        * 0.5.2   2022-08-19 [1] CRAN (R 4.2.0)
 fs               1.6.3   2023-07-20 [1] CRAN (R 4.2.0)
 gargle           1.2.1   2022-09-08 [1] CRAN (R 4.2.0)
 generics         0.1.3   2022-07-05 [1] CRAN (R 4.2.0)
 ggplot2        * 3.4.2   2023-04-03 [1] CRAN (R 4.2.0)
 glue             1.6.2   2022-02-24 [1] CRAN (R 4.2.0)
 googledrive      2.0.0   2021-07-08 [1] CRAN (R 4.2.0)
 googlesheets4    1.0.1   2022-08-13 [1] CRAN (R 4.2.0)
 gtable           0.3.3   2023-03-21 [1] CRAN (R 4.2.0)
 haven            2.5.1   2022-08-22 [1] CRAN (R 4.2.0)
 hms              1.1.2   2022-08-19 [1] CRAN (R 4.2.0)
 htmltools        0.5.5   2023-03-23 [1] CRAN (R 4.2.0)
 htmlwidgets      1.6.0   2022-12-15 [1] CRAN (R 4.2.2)
 httpuv           1.6.9   2023-02-14 [1] CRAN (R 4.2.0)
 httr             1.4.5   2023-02-24 [1] CRAN (R 4.2.0)
 igraph           1.3.5   2022-09-22 [1] CRAN (R 4.2.0)
 jsonlite         1.8.7   2023-06-29 [1] CRAN (R 4.2.0)
 knitr            1.43    2023-05-25 [1] CRAN (R 4.2.0)
 later            1.3.0   2021-08-18 [1] CRAN (R 4.2.0)
 lattice          0.20-45 2021-09-22 [1] CRAN (R 4.2.2)
 lifecycle        1.0.3   2022-10-07 [1] CRAN (R 4.2.0)
 lubridate        1.9.2   2023-02-10 [1] CRAN (R 4.2.0)
 magrittr         2.0.3   2022-03-30 [1] CRAN (R 4.2.0)
 memoise          2.0.1   2021-11-26 [1] CRAN (R 4.2.0)
 mime             0.12    2021-09-28 [1] CRAN (R 4.2.0)
 miniUI           0.1.1.1 2018-05-18 [1] CRAN (R 4.2.0)
 modelr           0.1.10  2022-11-11 [1] CRAN (R 4.2.0)
 munsell          0.5.0   2018-06-12 [1] CRAN (R 4.2.0)
 network          1.18.0  2022-10-06 [1] CRAN (R 4.2.0)
 phaseR         * 2.2.1   2022-09-02 [1] CRAN (R 4.2.0)
 pillar           1.9.0   2023-03-22 [1] CRAN (R 4.2.0)
 pkgbuild         1.4.0   2022-11-27 [1] CRAN (R 4.2.0)
 pkgconfig        2.0.3   2019-09-22 [1] CRAN (R 4.2.0)
 pkgload          1.3.2.1 2023-07-08 [1] CRAN (R 4.2.0)
 prettyunits      1.1.1   2020-01-24 [1] CRAN (R 4.2.0)
 processx         3.8.2   2023-06-30 [1] CRAN (R 4.2.0)
 profvis          0.3.7   2020-11-02 [1] CRAN (R 4.2.0)
 promises         1.2.0.1 2021-02-11 [1] CRAN (R 4.2.0)
 ps               1.7.5   2023-04-18 [1] CRAN (R 4.2.0)
 purrr          * 1.0.2   2023-08-10 [1] CRAN (R 4.2.0)
 R6               2.5.1   2021-08-19 [1] CRAN (R 4.2.0)
 Rcpp             1.0.10  2023-01-22 [1] CRAN (R 4.2.0)
 readr          * 2.1.3   2022-10-01 [1] CRAN (R 4.2.0)
 readxl           1.4.1   2022-08-17 [1] CRAN (R 4.2.0)
 remotes          2.4.2.1 2023-07-18 [1] CRAN (R 4.2.0)
 reprex           2.0.2   2022-08-17 [1] CRAN (R 4.2.0)
 rlang            1.1.1   2023-04-28 [1] CRAN (R 4.2.0)
 rstudioapi       0.14    2022-08-22 [1] CRAN (R 4.2.0)
 rvest            1.0.3   2022-08-19 [1] CRAN (R 4.2.0)
 scales           1.2.1   2022-08-20 [1] CRAN (R 4.2.0)
 sessioninfo      1.2.2   2021-12-06 [1] CRAN (R 4.2.0)
 shiny            1.7.4   2022-12-15 [1] CRAN (R 4.2.2)
 statnet.common   4.7.0   2022-09-08 [1] CRAN (R 4.2.0)
 stringi          1.7.12  2023-01-11 [1] CRAN (R 4.2.0)
 stringr        * 1.5.0   2022-12-02 [1] CRAN (R 4.2.0)
 tibble         * 3.2.1   2023-03-20 [1] CRAN (R 4.2.0)
 tidyr          * 1.3.0   2023-01-24 [1] CRAN (R 4.2.0)
 tidyselect       1.2.0   2022-10-10 [1] CRAN (R 4.2.0)
 tidyverse      * 1.3.2   2022-07-18 [1] CRAN (R 4.2.0)
 timechange       0.2.0   2023-01-11 [1] CRAN (R 4.2.0)
 tzdb             0.3.0   2022-03-28 [1] CRAN (R 4.2.0)
 urlchecker       1.0.1   2021-11-30 [1] CRAN (R 4.2.0)
 usethis        * 2.1.6   2022-05-25 [1] CRAN (R 4.2.0)
 utf8             1.2.3   2023-01-31 [1] CRAN (R 4.2.0)
 vctrs            0.6.3   2023-06-14 [1] CRAN (R 4.2.0)
 waydown        * 1.1.0   2021-03-17 [1] CRAN (R 4.2.0)
 withr            2.5.0   2022-03-03 [1] CRAN (R 4.2.0)
 xfun             0.40    2023-08-09 [1] CRAN (R 4.2.0)
 xml2             1.3.5   2023-07-06 [1] CRAN (R 4.2.0)
 xtable           1.8-4   2019-04-21 [1] CRAN (R 4.2.0)

 [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library

Method "ode45" not adapted for some models trajectories

phaseR/R/trajectory.R

Line 209 in 4f58357

method <- ifelse(tstep > 0, "ode45", "lsoda")

First of thank you for your time-saving package. I give a reproducible exemple below: for some initial conditions on this model, method "ode45" computing time for the trajectory is exponentially longer than method "euler". I am not sure whether the choice of "ode45" was guided by some particular reason? Maybe you could add the possibility to provide a method in the call to trajectory?

library(tictoc)
library(deSolve)

params <- c(s = 0.35, gamma = 2.5, n = 1.5, alpha = 3, beta = 2)

model <- function(t, y, parameters){
  with(as.list(parameters),{
    g <- y[1] 
    gw <- y[2]
    dg <- gamma*s*(1/gw-1/g)
    dgw <- -gw^2/s*(alpha*s*(1/g-1/gw) - beta*(g-n))
    
    list(c(dg,dgw))
  })
}
yini <- c(g = 0.27, gw = 1)
times <- seq(0, 5, by = 0.1)
tic()
out <- ode(yini, times, model, params, method = "ode45")
toc()
24.501 sec elapsed

tic()
out <- ode(yini, times, model, params, method = "euler")
toc()
0.003 sec elapsed

deriv function has inflexible `parameters` argument

Sorry for the vague title, it is a mistake in my own code:

library(phaseR)
#> -------------------------------------------------------------------------------
#> phaseR: Phase plane analysis of one- and two-dimensional autonomous ODE systems
#> -------------------------------------------------------------------------------
#> 
#> v.2.1: For an overview of the package's functionality enter: ?phaseR
#> 
#> For news on the latest updates enter: news(package = "phaseR")
ex2.alt <- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
      dA = A*(1 - A)*(2 - A)
      list(c(dA))
    })
}
example2_stability_1 <- stability(ex2.alt, ystar = 0, system = "one.dim", 
                                  state.names = "A")
#> discriminant = 2, classification = Unstable
                                  
ex2.alt.alt <- function(t, state, pars) {
    with(as.list(c(state, pars)), {
      dA = A*(1 - A)*(2 - A)
      list(c(dA))
    })
}
example2_stability_1 <- stability(ex2.alt.alt, ystar = 0, system = "one.dim", 
                                  state.names = "A")
#> Error in deriv(0, stats::setNames(ystar + h, state.names), parameters = parameters): unused argument (parameters = parameters)

^{Created on 2019-09-08 by the reprex package (v0.3.0)}

Apparently, it is not possible to change the name of the parameters argument in your deriv
The offending lines are https://github.com/mjg211/phaseR/blob/615b31c7ac/R/stability.R#L100-L101
Solution is easy: just change parameters = parameters into parameters (will make a pull request when I have time, but if you're earlier maybe you can take care of this when resolving #6 )

Possible typo in RJournal manuscript

Hi! Thanks for this fantastic package.

Not sure if this is the place for this, but I think there might be a typo in the associated R Journal manuscript (found here)

I think this is the other way around, right? (e.g. dx/dt = -x is stable at x=0, and k<0; dx/dt = x is unstable at x=0, and k>0)

Is it possible to plot the output of phaseR in ggplot?

Thank you a lot for all the great effort you have put in making phaseR and for the updates that you offer. I want to ask you if possible to plot the outcomes of the phaseR analysis into ggplot2?

thank you,
Loukas