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] = 3 * (xt + yt - xt^3/3 - parameters$k) dy[2] = (-1/3) * (xt + 0.8 * yt - 0.7) list(dy) } my_jacobian = function(x, y) { # returns the Jacobian at the point (x, y) # J = matrix(c(3* (1-x^2), 3, -1/3, -4/15), nrow=2, ncol=2, byrow=TRUE) } # Part (b): # #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_7_prob_21_part_b_plot.eps', onefile=FALSE, horizontal=FALSE) # k=0: # diff_eq_params = list(k=0) pr = polyroot (c(diff_eq_params$k - 7/8, 1/4, 0, 1/3)) cp_x = Re(pr[1]) cp_y = 7/8 - (5/4)*cp_x CP = data.frame(x=cp_x, y=cp_y) As = mapply(my_jacobian, CP$x, CP$y, SIMPLIFY=FALSE) for (ii in 1:length(As)){ es = eigen(As[[ii]]) print(sprintf('k= %f; CP: x= %f; y= %f', diff_eq_params$k, CP$x[ii], CP$y[ii])) print('Jacobian=') print(As[[ii]]) print('eigenvalues=') print(es$values) } par(mfrow=c(1, 2)) t.end = 20.0 flowField(my_yprime, x.lim = c(-5, 5), y.lim = c(-2, 7), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('k= %.2f', diff_eq_params$k) ) trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(CP$x, CP$y, pch=19) # k=0.5: # diff_eq_params = list(k=0.5) pr = polyroot (c(diff_eq_params$k - 7/8, 1/4, 0, 1/3)) cp_x = Re(pr[1]) cp_y = 7/8 - (5/4)*cp_x CP = data.frame(x=cp_x, y=cp_y) As = mapply(my_jacobian, CP$x, CP$y, SIMPLIFY=FALSE) for (ii in 1:length(As)){ es = eigen(As[[ii]]) print(sprintf('k= %f; CP: x= %f; y= %f', diff_eq_params$k, CP$x[ii], CP$y[ii])) print('Jacobian=') print(As[[ii]]) print('eigenvalues=') print(es$values) } flowField(my_yprime, x.lim = c(-5, 5), y.lim = c(-2, 7), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('k= %.2f', diff_eq_params$k) ) trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(CP$x, CP$y, pch=19) par(mfrow=c(1, 1)) #dev.off() # Part (c): # ks = seq( 0.25, 0.45, length.out=20 ) xs = c() ys = c() eigen_1 = c() eigen_2 = c() for( k in ks ){ pr = polyroot(c(k - 7/8, 1/4, 0, 1/3)) cp_x = Re(pr[1]) cp_y = 7/8 - (5/4) *cp_x J = my_jacobian(cp_x, cp_y) es = eigen(J) xs = c(xs, cp_x) ys = c(ys, cp_y) eigen_1 = c(eigen_1, es$values[1]) eigen_2 = c(eigen_2, es$values[2]) } R = data.frame(k=ks, CP_X=xs, CP_Y=ys, eigen_1=eigen_1, eigen_2=eigen_2) print(round(R, 3)) k0 = 0.346 #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_7_prob_21_part_c_plot.eps', onefile=FALSE, horizontal=FALSE) diff_eq_params = list(k=k0) pr = polyroot (c(diff_eq_params$k - 7/8, 1/4, 0, 1/3)) cp_x = Re(pr[1]) cp_y = 7/8 - (5/4)*cp_x CP = data.frame(x=cp_x, y=cp_y) As = mapply(my_jacobian, CP$x, CP$y, SIMPLIFY=FALSE) for (ii in 1:length(As)){ es = eigen(As[[ii]]) print(sprintf('k= %f; CP: x= %f; y= %f', diff_eq_params$k, CP$x[ii], CP$y[ii])) print('Jacobian=') print(As[[ii]]) print('eigenvalues=') print(es$values) } t.end = 20.0 flowField(my_yprime, x.lim = c(-5, 5), y.lim = c(-2, 7), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('k= %.2f', diff_eq_params$k) ) trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(CP$x, CP$y, pch=19) #dev.off() # Part (d): # ##postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_7_prob_21_part_d_phase_plots.eps', onefile=FALSE, horizontal=FALSE) t.end = 30.0 ks = c(0.4, 0.5, 0.6) par(mfrow=c(1, length(ks))) for( k in ks ){ diff_eq_params = list(k=k) pr = polyroot(c(diff_eq_params$k - 7/8, 1/4, 0, 1/3)) cp_x = Re(pr[1]) cp_y = 7/8 - (5/4)*cp_x CP = data.frame(x=cp_x, y=cp_y) flowField(my_yprime, x.lim = c(-3.0, 6), y.lim = c(-0.5, 7), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x', ylab='y', main=sprintf('k= %.2f', diff_eq_params$k) ) trajectory(my_yprime, y0 = c(3, 6), t.end = t.end, parameters = diff_eq_params, col='blue', pch=19) points(CP$x, CP$y, pch=19) } par(mfrow=c(1, 1)) ##dev.off() # Plot x(t) as a function of t: # t.end = 50.0 ts = seq(0, t.end, length.out=500) ks = c(0.4, 0.5, 0.6) Xs = c() # the numerical solutions: Ts = c() # the estimated periods for( k in ks ){ diff_eq_params = list(k=k) res = ode(y=c(3, 6), times=ts, fun=my_yprime, parms=diff_eq_params) xs = res[, 2] # the numerical function x(t) Xs = cbind(Xs, xs) Ts = c(Ts, my_period(ts, Xs) ) } colnames(Xs) = ks ##postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_7_prob_21_part_d_x_vs_t_plots.eps', onefile=FALSE, horizontal=FALSE) matplot(ts, Xs, type='l', lty=1, lwd=2, xlab='t', ylab='x(t)') grid() legend('topright', sprintf('k= %.2f', ks), lty=c(1, 1, 1), lwd=2, col=c('black','red','green') ) ##dev.off() # Display the k and period values: # R = data.frame(k=ks, period=Ts) print(R) # Part (e): # ks = seq( 3.3, 4.0, length.out=20 ) xs = c() ys = c() eigen_1 = c() eigen_2 = c() for( k in ks ){ pr = polyroot(c(k - 7/8, 1/4, 0, 1/3)) cp_x = Re(pr[1]) cp_y = 7/8 - (5/4) *cp_x J = my_jacobian(cp_x, cp_y) es = eigen(J) xs = c(xs, cp_x) ys = c(ys, cp_y) eigen_1 = c(eigen_1, es$values[1]) eigen_2 = c(eigen_2, es$values[2]) } R = data.frame(k=ks, CP_X=xs, CP_Y=ys, eigen_1=eigen_1, eigen_2=eigen_2) print(round(R, 3))