145 lines
3.6 KiB
Plaintext
145 lines
3.6 KiB
Plaintext
---
|
|
title: "BCR 06 Sept"
|
|
author: "Jonathan Herrewijnen"
|
|
date: "September 6, 2018"
|
|
output: html_document
|
|
---
|
|
|
|
```{r setup, include=FALSE}
|
|
knitr::opts_chunk$set(echo = TRUE)
|
|
```
|
|
|
|
## R Markdown
|
|
|
|
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.
|
|
|
|
When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
|
|
|
|
```{r cars}
|
|
summary(cars)
|
|
```
|
|
|
|
## Including Plots
|
|
|
|
You can also embed plots, for example:
|
|
|
|
```{r pressure, echo=FALSE}
|
|
plot(pressure)
|
|
```
|
|
|
|
Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
|
|
|
|
//CBR Two dimensional phase portraits, equilibria and stability
|
|
|
|
```{r}
|
|
#install.packages("phaseR")
|
|
library("phaseR")
|
|
|
|
# define 2D ODE model
|
|
model = function(timepoint, state, parameters) {
|
|
with(as.list(c(state, parameters)), {
|
|
dx = -2 * x - 3 * y
|
|
dy = -x - 2 * y
|
|
list(c(dx, dy))
|
|
})
|
|
}
|
|
# for setting the axis limits
|
|
xlim = c(-10, 10)
|
|
ylim = c(-10, 10)
|
|
# no parameters in the ODE, so we set pars as an empty vector
|
|
pars = vector()
|
|
|
|
# draw vector field
|
|
model.flowField = flowField(xlab = "x", ylab = "y", model, xlim = xlim, ylim = ylim,
|
|
parameters = pars, points = 19, add = FALSE, state.names = c("x", "y"))
|
|
|
|
# draw nullclines
|
|
model.nullclines = nullclines(model, xlim = xlim, ylim = ylim, parameters = pars,
|
|
points = 500, state.names = c("x", "y"))
|
|
|
|
# draw trajectory
|
|
model.trajectory = trajectory(model, y0 = c(-5, -5), tlim = c(0, 10), parameters = pars,
|
|
state.names = c("x", "y"))
|
|
|
|
```
|
|
|
|
//Numerical Solution
|
|
|
|
```{r}
|
|
# perform numerical simulation and plot end result
|
|
model.numericalSolution = numericalSolution(model, y0 = c(-5, -5), tlim = c(0,
|
|
20), type = "one", parameters = pars, col = c("green", "orange"), ylim = ylim,
|
|
add.legend = F)
|
|
|
|
# make legend for plot
|
|
legend("topright", col = c("green", "orange"), legend = c("x", "y"), lty = 1)
|
|
|
|
# after the next command, click mouse to find values of a point within the
|
|
# phase portrait (e.g., to approximate the value of an equilibrium)
|
|
|
|
eq = locator(n = 1)
|
|
eq
|
|
```
|
|
|
|
//QUESTION 1
|
|
|
|
The 0, 0 location was determined from the first phase portrait!
|
|
Not from the second click portrait!
|
|
|
|
```{r}
|
|
#?stability
|
|
|
|
stability(model, ystar = c(0, 0))
|
|
#stability(deriv, ystar = NULL, parameters = NULL, system = "two.dim", h = 1e-07, summary = TRUE, state.names = c("x", "y"))
|
|
|
|
|
|
```
|
|
|
|
//Question 2
|
|
Meaning of T and Delta.
|
|
|
|
T=(I THINK!) the value of the equilibrium at dx/dt=0 -> So for the equilibrium, in our case the equilibrium is stable.
|
|
Delta="In the two dimensional system case, value of the Jacobians determinant at y.star. ->So the Jacobian matrix
|
|
|
|
//Question 3: PhasePlaneAnalysis
|
|
|
|
```{r}
|
|
# perform phase plane analysis
|
|
phasePlaneAnalysis(model, xlim = xlim, ylim = ylim, parameters = pars)
|
|
```
|
|
|
|
//QUESTION 4-8
|
|
```{r}
|
|
library("phaseR")
|
|
# define 2D ODE model
|
|
model = function(timepoint, state, parameters) {
|
|
with(as.list(c(state, parameters)), {
|
|
|
|
#Question 4
|
|
dx = -x + (2 * y)
|
|
dy = (-5 * x) - y
|
|
|
|
#Question 5
|
|
#dx = -3*x -2*y
|
|
#dy = -2*x -y
|
|
|
|
list(c(dx, dy))
|
|
})
|
|
}
|
|
|
|
# for setting the axis limits
|
|
xlim = c(-10, 10)
|
|
ylim = c(-10, 10)
|
|
|
|
# no parameters in the ODE, so we set pars as an empty vector
|
|
pars = vector()
|
|
|
|
#Determine Stability
|
|
stability(model, ystar = c(0, 0))
|
|
|
|
# perform phase plane analysis
|
|
phasePlaneAnalysis(model, xlim = xlim, ylim = ylim, parameters = pars)
|
|
```
|
|
|
|
|