if( !require('phaseR') ){ install.packages('phaseR') } library(phaseR) my_jacobian = function(x, y) { # returns the Jacobian at the point (x, y) # J = matrix(c(-1 + 2*y, 2*x, -2*x, 1-2*y), nrow=2, ncol=2, byrow=TRUE) } # The critical points of this system: CP = data.frame(x=c(-1/2, 0, 0, 1/2), y=c(1/2, 0, 1, 1/2)) As = mapply(my_jacobian, CP\$x, CP\$y, SIMPLIFY=FALSE) for (ii in 1:length(As)){ es = eigen(As[[ii]]) print(sprintf('CP: x= %f; y= %f', CP\$x[ii], CP\$y[ii])) print('Jacobian=') print(As[[ii]]) print('eigenvalues=') print(es\$values) } my_yprime = function(t, y, parameters) { xt = y[1] yt = y[2] dy = rep(NA, length(y)) dy[1] = -xt + 2*xt*yt dy[2] = yt - xt^2 - yt^2 list(dy) } diff_eq_params = list() #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_2_prob_11_plot.eps', onefile=FALSE, horizontal=FALSE) t.end = 0.5 L = 2 flowField(my_yprime, x.lim = c(-L, +L), y.lim = c(-L, +L), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x_1', ylab='x_2') trajectory(my_yprime, y0 = c(1.0, 1.0), t.end = t.end, parameters = diff_eq_params, col='black', pch=19) trajectory(my_yprime, y0 = c(1.0, -1.0), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) trajectory(my_yprime, y0 = c(-1.0, 1.0), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory(my_yprime, y0 = c(-1.0, -1.0), t.end = t.end, parameters = diff_eq_params, col='green', pch=19) #dev.off()