-
Notifications
You must be signed in to change notification settings - Fork 3
/
trajectory.R
263 lines (263 loc) · 11.6 KB
/
trajectory.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
#' Phase plane trajectory plotting
#'
#' Performs numerical integration of the chosen ODE system, for a user specified
#' set of initial conditions. Plots the resulting solution(s) in the phase
#' plane.
#'
#' @param deriv A function computing the derivative at a point for the ODE
#' system to be analysed. Discussion of the required structure of these
#' functions can be found in the package vignette, or in the help file for the
#' function \code{\link[deSolve]{ode}}.
#' @param y0 The initial condition(s). In the case of a one-dimensional system,
#' this can either be a \code{\link[base]{numeric}} \code{\link[base]{vector}}
#' of \code{\link[base]{length}} one, indicating the location of the dependent
#' variable initially, or a \code{\link[base]{numeric}}
#' \code{\link[base]{vector}} indicating multiple initial locations of the
#' independent variable. In the case of a two-dimensional system, this can
#' either be a \code{\link[base]{numeric}} \code{\link[base]{vector}} of
#' \code{\link[base]{length}} two, reflecting the location of the two
#' dependent variables initially, or it can be \code{\link[base]{numeric}}
#' \code{\link[base]{matrix}} where each row reflects a different initial
#' condition. Alternatively this can be specified as \code{\link[base]{NULL}},
#' and then \code{\link[graphics]{locator}} can be used to specify initial
#' condition(s) on a plot. In this case, for one-dimensional systems, all
#' initial conditions are taken at tlim[1], even if not selected so on the
#' graph. Defaults to \code{\link[base]{NULL}}.
#' @param n If \code{y0} is left \code{NULL}, such initial conditions can be
#' specified using \code{\link[graphics]{locator}}, \code{n} sets the number of
#' initial conditions to be chosen. Defaults to \code{NULL}.
#' @param tlim Sets the limits of the independent variable for which the
#' solution should be plotted. Should be a \code{\link[base]{numeric}}
#' \code{\link[base]{vector}} of \code{\link[base]{length}} two. If
#' \code{tlim[2] > tlim[1]}, then \code{tstep} should be negative to indicate a
#' backwards trajectory.
#' @param tstep The step length of the independent variable, used in numerical
#' integration. Decreasing the absolute magnitude of \code{tstep} theoretically
#' makes the numerical integration more accurate, but increases computation
#' time. Defaults to \code{0.01}.
#' @param parameters Parameters of the ODE system, to be passed to \code{deriv}.
#' Supplied as a \code{\link[base]{numeric}} \code{\link[base]{vector}}; the
#' order of the parameters can be found from the \code{deriv} file. Defaults to
#' \code{NULL}.
#' @param system Set to either \code{"one.dim"} or \code{"two.dim"} to indicate
#' the type of system being analysed. Defaults to \code{"two.dim"}.
#' @param col The colour(s) to plot the trajectories in. Should be a
#' \code{\link[base]{character}} \code{\link[base]{vector}}. Will be reset
#' accordingly if it is not of the \code{\link[base]{length}} of the number of
#' initial conditions. Defaults to \code{"black"}.
#' @param add Logical. If \code{TRUE}, the trajectories added to an existing
#' plot. If \code{FALSE}, a new plot is created. Defaults to \code{TRUE}.
#' @param \dots Additional arguments to be passed to plot.
#' @inheritParams .paramDummy
#' @return Returns a list with the following components (the exact make up is
#' dependent on the value of \code{system}):
#' \item{add}{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{n}{As per input.}
#' \item{parameters}{As per input.}
#' \item{system}{As per input.}
#' \item{tlim}{As per input.}
#' \item{tstep}{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.}
#' \item{x}{In the two-dimensional system case, a \code{\link[base]{numeric}}
#' \code{\link[base]{matrix}} whose columns are the numerically computed values
#' of the first dependent variable for each initial condition.}
#' \item{y}{In the two-dimensional system case, a \code{\link[base]{numeric}}
#' \code{\link[base]{matrix}} whose columns are the numerically computed values
#' of the second dependent variable for each initial condition.
#' In the one-dimensional system case, a \code{\link[base]{numeric}}
#' \code{\link[base]{matrix}} whose columns are the numerically computed values
#' of the dependent variable for each initial condition.}
#' \item{y0}{As per input, but converted to a \code{\link[base]{numeric}}
#' \code{\link[base]{matrix}} if supplied as a vector initially.}
#' @author Michael J Grayling
#' @seealso \code{\link[deSolve]{ode}}, \code{\link[graphics]{plot}}
#' @export
#' @examples
#' # Plot the flow field, nullclines and several trajectories for the
#' # one-dimensional autonomous ODE system logistic
#' logistic_flowField <- flowField(logistic,
#' xlim = c(0, 5),
#' ylim = c(-1, 3),
#' parameters = c(1, 2),
#' points = 21,
#' system = "one.dim",
#' add = FALSE)
#' logistic_nullclines <- nullclines(logistic,
#' xlim = c(0, 5),
#' ylim = c(-1, 3),
#' parameters = c(1, 2),
#' system = "one.dim")
#' logistic_trajectory <- trajectory(logistic,
#' y0 = c(-0.5, 0.5, 1.5, 2.5),
#' tlim = c(0, 5),
#' parameters = c(1, 2),
#' system = "one.dim")
#'
#' # Plot the velocity field, nullclines and several trajectories for the
#' # two-dimensional autonomous ODE system simplePendulum
#' simplePendulum_flowField <- flowField(simplePendulum,
#' xlim = c(-7, 7),
#' ylim = c(-7, 7),
#' parameters = 5,
#' points = 19,
#' add = FALSE)
#' y0 <- matrix(c(0, 1, 0, 4, -6, 1, 5, 0.5, 0, -3),
#' 5, 2, byrow = TRUE)
#' \donttest{
#' simplePendulum_nullclines <- nullclines(simplePendulum,
#' xlim = c(-7, 7),
#' ylim = c(-7, 7),
#' parameters = 5,
#' points = 500)
#' }
#' simplePendulum_trajectory <- trajectory(simplePendulum,
#' y0 = y0,
#' tlim = c(0, 10),
#' parameters = 5)
trajectory <- function(deriv, y0 = NULL, n = NULL, tlim, tstep = 0.01,
parameters = NULL, system = "two.dim", col = "black",
add = TRUE,
state.names =
if (system == "two.dim") c("x", "y") else "y", ...) {
if (tstep == 0) {
stop("tstep is equal to 0")
}
if (tlim[1] == tlim[2]) {
stop("tlim[1] is equal to tlim[2]")
}
if (all(tlim[1] > tlim[2], tstep > 0)) {
stop("tstep must be negative if tlim[1] > tlim[2]")
}
if (all(tlim[1] < tlim[2], tstep < 0)) {
stop("tstep must be positive if tlim[1] < tlim[2]")
}
if (!(system %in% c("one.dim", "two.dim"))) {
stop("system must be set to either \"one.dim\" or \"two.dim\"")
}
if (!is.vector(col)) {
stop("col is not a vector as required")
}
if (!is.logical(add)) {
stop("add must be logical")
}
if (all(is.null(y0), is.null(n))) {
stop("Both y0 and n cannot be NULL")
}
if (all(!is.null(y0), !is.null(n))) {
warning("n is non-NULL whilst y0 has also been specified")
}
if (all(is.null(y0), !add)) {
stop("y0 cannot be NULL when add is set to FALSE")
}
if (is.null(y0)) {
y0 <- locator(n = n)
if (system == "two.dim") {
re.set <- matrix(0, ncol = 2, nrow = n)
for (i in 1:n) {
re.set[i, ] <- c(y0$x[i], y0$y[i])
}
y0 <- re.set
}
if (system == "one.dim") {
re.set <- numeric(n)
for (i in 1:n) {
re.set[i] <- y0$y[i]
}
y0 <- re.set
}
}
if (all(!is.vector(y0), !is.matrix(y0))) {
stop("y0 is neither a number, vector or matrix as required")
}
if (is.vector(y0)) {
y0 <- as.matrix(y0)
}
if (all(system == "one.dim", dim(y0) > 1)) {
stop("For system equal to \"one.dim\", y0 must contain either a vector or ",
"a matrix where either nrow(y0) or ncol(y0) is 1")
}
if (all(system == "two.dim", !any(dim(y0) == 2))) {
stop("For system equal to \"two.dim\", y0 must contain either a vector of ",
"length two or a matrix where either nrow(y0) or ncol(y0) is 2")
}
if (system == "one.dim") {
if (ncol(y0) > nrow(y0)) {
y0 <- t(y0)
}
} else {
if (all(nrow(y0) == 2, ncol(y0) != 2)) {
y0 <- t(y0)
}
}
if (nrow(y0) > length(col)) {
col <- rep(col, nrow(y0))
message("Note: col has been reset as required")
} else if (nrow(y0) < length(col)) {
col <- col[1:nrow(y0)]
message("Note: col has been reset as required")
}
t <- seq(tlim[1], tlim[2], tstep)
x <- matrix(0, length(t), nrow(y0))
if (system == "two.dim") {
y <- x
}
method <- ifelse(tstep > 0, "ode45", "lsoda")
for (i in 1:nrow(y0)) {
phase.trajectory <- deSolve::ode(times = t,
y = stats::setNames(c(y0[i, ]),
state.names),
func = deriv,
parms = parameters,
method = method)
x[, i] <- phase.trajectory[, 2]
if (system == "two.dim") {
y[, i] <- phase.trajectory[, 3]
}
if (all(!add, i == 1)) {
if (system == "one.dim") {
graphics::plot(t, x[, i], col = col[i], type = "l", ...)
} else {
graphics::plot(x[, i], y[, i], col = col[i], type = "l", ...)
}
} else {
if (system == "one.dim") {
lines(t, x[, i], col = col[i], type = "l", ...)
} else {
lines(x[, i], y[, i], col = col[i], type = "l", ...)
}
}
}
if (system == "one.dim") {
graphics::points(rep(tlim[1], nrow(y0)), y0, col = col, ...)
return(list(add = add,
col = col,
deriv = deriv,
n = n,
parameters = parameters,
system = system,
t = t,
tlim = tlim,
tstep = tstep,
y = x,
y0 = y0))
} else {
graphics::points(y0[, 1], y0[, 2], col = col, ...)
return(list(add = add,
col = col,
deriv = deriv,
n = n,
parameters = parameters,
system = system,
t = t,
tlim = tlim,
tstep = tstep,
x = x,
y = y,
y0 = y0))
}
}