Added method as an option where relevant

mjg211 · Nov 26, 2019 · 4a9b72e · 4a9b72e
1 parent 4f58357
commit 4a9b72e
Show file tree

Hide file tree

Showing 10 changed files with 55 additions and 23 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: phaseR
 Type: Package
 Title: Phase Plane Analysis of One- And Two-Dimensional Autonomous ODE Systems
-Version: 2.1.3
+Version: 2.2
 Authors@R: c(person(given   = "Michael J",
                     family  = "Grayling",
                     email   = "michael.grayling@newcastle.ac.uk",

diff --git a/R/drawManifolds.R b/R/drawManifolds.R
@@ -32,13 +32,16 @@
 #' @param add.legend Logical. If \code{TRUE}, a legend is added to the plots.
 #' Defaults to \code{TRUE}.
 #' @inheritParams .paramDummy
+#' @param method Passed to \code{\link[deSolve]{ode}}. See there for further
+#' details. Defaults to \code{"lsoda"}.
 #' @param ... Additional arguments to be passed to plot.
 #' @return Returns a \code{\link[base]{list}} with the following components:
 #' \item{add.legend}{As per input.}
 #' \item{col}{As per input, but with possible editing if a
 #' \code{\link[base]{character}} \code{\link[base]{vector}} of the wrong
 #' \code{\link[base]{length}} was supplied.}
 #' \item{deriv}{As per input.}
+#' \item{method}{As per input.}
 #' \item{parameters}{As per input.}
 #' \item{stable.1}{A \code{\link[base]{numeric}} \code{\link[base]{matrix}}
 #' whose columns are the numerically computed values of the dependent variables
@@ -59,7 +62,8 @@
 #' @export
 drawManifolds <- function(deriv, y0 = NULL, parameters = NULL, tstep = 0.1,
                           tend = 100, col = c("green", "red"),
-                          add.legend = TRUE, state.names = c("x", "y"), ...) {
+                          add.legend = TRUE, state.names = c("x", "y"),
+                          method = "lsoda", ...) {
   if (is.null(y0)) {
     y0         <- locator(n = 1)
     y0         <- c(y0$x, y0$y)
@@ -109,29 +113,33 @@ drawManifolds <- function(deriv, y0 = NULL, parameters = NULL, tstep = 0.1,
     eps        <- 1e-2
     ymax       <- 0.5 + max(abs(ystar))
     maxtime.1  <- max(tend, log(2500*ymax)/abs(eigenvalues[i1]))
-    out.1      <- deSolve::ode(times = seq(0, maxtime.1, tstep),
-                               y     = setNames(ystar + eps*v1, state.names),
-                               func  = deriv,
-                               parms = parameters)
+    out.1      <- deSolve::ode(times  = seq(0, maxtime.1, tstep),
+                               y      = setNames(ystar + eps*v1, state.names),
+                               func   = deriv,
+                               parms  = parameters,
+                               method = method)
     graphics::lines(out.1[, 2], out.1[, 3], type = "l", col = col[2], ...)
     unstable.1 <- out.1[, 1:3]
-    out.2      <- deSolve::ode(times = seq(0, maxtime.1, tstep),
-                               y     = setNames(ystar - eps*v1, state.names),
-                               func  = deriv,
-                               parms = parameters)
+    out.2      <- deSolve::ode(times  = seq(0, maxtime.1, tstep),
+                               y      = setNames(ystar - eps*v1, state.names),
+                               func   = deriv,
+                               parms  = parameters,
+                               method = method)
     graphics::lines(out.2[, 2], out.2[, 3], type = "l", col = col[2], ...)
     unstable.2 <- out.2[, 1:3]
     maxtime.2  <- max(tend, log(2500*ymax)/abs(eigenvalues[i2]))
-    out.3      <- deSolve::ode(times = -seq(0, maxtime.2, tstep),
-                               y     = setNames(ystar + eps*v2, state.names),
-                               func  = deriv,
-                               parms = parameters)
+    out.3      <- deSolve::ode(times  = -seq(0, maxtime.2, tstep),
+                               y      = setNames(ystar + eps*v2, state.names),
+                               func   = deriv,
+                               parms  = parameters,
+                               method = method)
     graphics::lines(out.3[, 2], out.3[, 3], type = "l", col = col[1], ...)
     stable.1   <- out.3[, 1:3]
-    out.4      <- deSolve::ode(times = -seq(0, maxtime.2, tstep),
-                               y     = setNames(ystar - eps*v2, state.names),
-                               func  = deriv,
-                               parms = parameters)
+    out.4      <- deSolve::ode(times  = -seq(0, maxtime.2, tstep),
+                               y      = setNames(ystar - eps*v2, state.names),
+                               func   = deriv,
+                               parms  = parameters,
+                               method = method)
     graphics::lines(out.4[, 2], out.4[, 3], type = "l", col = col[1], ...)
     stable.2   <- out.4[, 1:3]
     if (add.legend) {
@@ -141,6 +149,7 @@ drawManifolds <- function(deriv, y0 = NULL, parameters = NULL, tstep = 0.1,
     return(list(add.legend = add.legend,
                 col        = col,
                 deriv      = deriv,
+                method     = method,
                 parameters = parameters,
                 stable.1   = stable.1,
                 stable.2   = stable.2,

diff --git a/R/numericalSolution.R b/R/numericalSolution.R
@@ -40,6 +40,8 @@
 #' @inheritParams .paramDummy
 #' @param xlab Label for the x-axis of the resulting plot.
 #' @param ylab Label for the y-axis of the resulting plot.
+#' @param method Passed to \code{\link[deSolve]{ode}}. See there for further
+#' details. Defaults to \code{"ode45"}.
 #' @param \dots Additional arguments to be passed to
 #' \code{\link[graphics]{plot}}.
 #' @return Returns a \code{\link[base]{list}} with the following components:
@@ -49,6 +51,7 @@
 #' \code{\link[base]{character}} \code{\link[base]{vector}} of the wrong
 #' \code{\link[base]{length}} was supplied.}
 #' \item{deriv}{As per input.}
+#' \item{method}{As per input.}
 #' \item{parameters}{As per input.}
 #' \item{t}{A \code{\link[base]{numeric}} \code{\link[base]{vector}} containing
 #' the values of the independent variable at each integration step.}
@@ -75,7 +78,8 @@ numericalSolution <- function(deriv, y0 = NULL, tlim, tstep = 0.01,
                               parameters = NULL, type = "one",
                               col = c("red", "blue"), add.grid = TRUE,
                               add.legend = TRUE, state.names = c("x", "y"),
-                              xlab = "t", ylab = state.names, ...) {
+                              xlab = "t", ylab = state.names, method = "ode45",
+                              ...) {
   if (any(tlim < 0)) {
     stop("tlim contains negative values")
   }
@@ -133,7 +137,7 @@ numericalSolution <- function(deriv, y0 = NULL, tlim, tstep = 0.01,
                                    y      = stats::setNames(y0, state.names),
                                    func   = deriv,
                                    parms  = parameters,
-                                   method = "ode45")
+                                   method = method)
   x                <- phase.trajectory[, 2]
   y                <- phase.trajectory[, 3]
   if (type == "one") {
@@ -167,6 +171,7 @@ numericalSolution <- function(deriv, y0 = NULL, tlim, tstep = 0.01,
               add.legend = add.legend,
               col        = col,
               deriv      = deriv,
+              method     = method,
               parameters = parameters,
               t          = t,
               tlim       = tlim,

diff --git a/man/drawManifolds.Rd b/man/drawManifolds.Rd
diff --git a/man/flowField.Rd b/man/flowField.Rd
diff --git a/man/nullclines.Rd b/man/nullclines.Rd
diff --git a/man/numericalSolution.Rd b/man/numericalSolution.Rd
diff --git a/man/simplePendulum.Rd b/man/simplePendulum.Rd
diff --git a/man/trajectory.Rd b/man/trajectory.Rd
diff --git a/man/vanDerPol.Rd b/man/vanDerPol.Rd