if( !require('phaseR') ){ install.packages('phaseR') } library(deSolve) library(phaseR) A = matrix(c(-4/5, 2, -1, 6/5), nrow=2, ncol=2, byrow=T) print(eigen(A)) # Part (b) Plot (x1, x2): # my_yprime = function(t, y, parameters){ dy = parameters\$A %*% y list(dy) } diff_eq_params = list(A=A) #postscript("../../WriteUp/Graphics/Chapter7/chap_7_sect_6_prob_12_x1_x2_plot.eps", onefile=FALSE, horizontal=FALSE) flowField(my_yprime, x.lim=c(-10, 10), y.lim=c(-10, 10), parameters = diff_eq_params, points=21, add=FALSE, xlab='x_1', ylab='x_2', main='Problem 12') trajectory(my_yprime, y0 = c(3, 8), t.end = 3.0, parameters = diff_eq_params, col='black') trajectory(my_yprime, y0 = c(-2.5, -5), t.end = 2.0, parameters = diff_eq_params, col='black') trajectory(my_yprime, y0 = c(5, -5), t.end = 1.0, parameters = diff_eq_params, col='black') trajectory(my_yprime, y0 = c(-2, 2), t.end = 1.0, parameters = diff_eq_params, col='black') grid() #dev.off() # Part (b) Plot (t, x1(t)) and (t, x2(t)): # t = seq( 0, 15, length.out=1000 ) state = c( y=c(3, 8) ) parms = list(A=A) res = ode(y=state, times=t, fun=my_yprime, parms=parms) ts = res[,'time'] #postscript("../../WriteUp/Graphics/Chapter7/chap_7_sect_6_prob_12_t_x1_N_t_x2.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 2)) plot(ts, res[,'y1'], type='l', lty=1, xlab='t', ylab='y_1(t)') lines(c(0), c(3), type='p', pch=19, cex=1.5) grid() plot(ts, res[,'y2'], type='l', lty=1, xlab='t', ylab='y_2(t)') lines(c(0), c(8), type='p', pch=19, cex=1.5) grid() par(mfrow=c(1, 1)) #dev.off()