if( !require('phaseR') ){ install.packages('phaseR') } library(phaseR) if( !require('nleqslv') ){ install.packages('nleqslv') } library(nleqslv) my_jacobian = function(x, y) { # returns the Jacobian at the point (x, y) # J = matrix(c(-2 - 3*x^2 - y^2, -1 - 2*x*y, 1 + 2*x*y, -1 + x^2 + 3*y^2), nrow=2, ncol=2, byrow=TRUE) } my_yprime = function(t, y, parameters) { xt = y[1] yt = y[2] dy = rep(NA, length(y)) dy[1] = -2*xt - yt - xt*(xt^2 + yt^2) dy[2] = xt - yt + yt*(xt^2 + yt^2) list(dy) } diff_eq_params = list() print('Solving f(x) = 0 numerically to find the roots:') rts_1 = nleqslv(c(1, 1), function(x){my_yprime(NA, x, diff_eq_params)[[1]]}, jac=function(x){my_jacobian(x[1], x[2])}) print(rts_1\$x) rts_2 = nleqslv(c(-1, -1), function (x){my_yprime(NA, x, diff_eq_params)[[1]]}, jac=function(x){my_jacobian(x[1], x[2])}) print(rts_2\$x) # Classify the critical points of this system: CP = data.frame (x=c(rts_1\$x[1], 0, rts_2\$x[1]), y=c(rts_1\$x[2], 0, rts_2\$x[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) } #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_3_prob_15_plot.eps', onefile=FALSE, horizontal=FALSE) t.end = 5.0 L = 4 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') # Plot some trajectories with initial condition near the critical points: # trajectory(my_yprime, y0 = c(1, 1), t.end = t.end, parameters = diff_eq_params, col='black', pch=19) trajectory(my_yprime, y0 = c(-1, 1), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) trajectory(my_yprime, y0 = c(1, -1), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) trajectory(my_yprime, y0 = c(-1, -1), t.end = t.end, parameters = diff_eq_params, col='green', pch=19) #dev.off()