if( !require('phaseR') ){ install.packages('phaseR') } library(phaseR) library(deSolve) source('../Chapter3/utils.R') my_yprime = function(t, y, parameters) { xt = y[1] yt = y[2] dy = rep(NA, length(y)) dy[1] = yt dy[2] = -xt + parameters$mu * ( 1 - xt^2 ) * yt list(dy) } # Parts (b-c): # #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_7_prob_17_part_b_N_c_plot.eps', onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 3)) t.end = 20.0 L = 4 diff_eq_params = list(mu=-1) flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('mu= %.2f', diff_eq_params$mu) ) trajectory(my_yprime, y0 = c(0.5, 1.75), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory (my_yprime, y0 = c(-0.5, -1.5), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(0, 0, pch=19) diff_eq_params = list(mu=-2) flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('mu= %.2f', diff_eq_params$mu) ) trajectory(my_yprime, y0 = c(0.5, 1.75), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory(my_yprime, y0 = c(-0.5, -1.5), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(0, 0, pch=19) diff_eq_params = list(mu=-0.5) flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('mu= %.2f', diff_eq_params$mu) ) trajectory(my_yprime, y0 = c(0.5, 1.75), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory(my_yprime, y0 = c(-0.5, -1.5), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(0, 0, pch=19) par(mfrow=c(1, 1)) #dev.off() # Part (d): # #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_7_prob_17_part_d_plot.eps', onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 3)) t.end = 10.0 L = 4 diff_eq_params = list(mu=1.0) flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('mu= %.2f', diff_eq_params$mu) ) trajectory(my_yprime, y0 = c(0.5, 1.75), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory(my_yprime, y0 = c(-0.5, -1.5), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(0, 0, pch=19) diff_eq_params = list(mu=0.1) flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('mu= %.2f', diff_eq_params$mu) ) trajectory(my_yprime, y0 = c(0.5, 1.75), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory(my_yprime, y0 = c(-0.5, -1.5), t.end = t.end, , parameters = diff_eq_params, col='blue', pch=19) points(0, 0, pch=19) diff_eq_params = list(mu=0.01) flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('mu= %.2f', diff_eq_params$mu) ) trajectory(my_yprime, y0 = c(0.5, 1.75), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory(my_yprime, y0 = c(-0.5, -1.5), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(0, 0, pch=19) par(mfrow=c(1, 1)) #dev.off()