# Problem 11: # polyroot( c(15, 3/10, 1) ) # Problem 12: # polyroot( c(5e5, 1500, 1) ) A = matrix(data=c(1, 1, -500, -1000), nrow=2, ncol=2, byrow=TRUE) solve(A, c(1e-6, 0)) # Problem 28: # t = seq(0, 10, length.out=100) ut = sqrt(2)*sin(sqrt(2)*t) upt = 2*cos(sqrt(2)*t) #postscript("../../WriteUp/Graphics/Chapter3/chap_3_sect_7_prob_28_plot.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 2)) ymin = min( c(ut, upt) ) ymax = max( c(ut, upt) ) plot(t, ut, type='l', col='blue', ylim=c(ymin, ymax), ylab='u and u prime', main='time series') lines(t, upt, type='l', col='red') legend('bottomright', c('u', 'u prime'), lty=c(1, 1), col=c('blue', 'red')) grid() plot(ut, upt, type='l', col='black', xlab='u', ylab='u prime', main='phase plot') grid() par(mfrow=c(1, 1)) #dev.off() # Problem 29: # polyroot( c(2, 0.25, 1) ) t = seq(0, 10, length.out=100) A = sqrt(127)/8 ut = (2/A) * exp(-t/8) * sin(A*t) upt = -(1/(4*A)) * exp(-t/8) * sin(A*t) + 2 * exp(-t/8) * cos(A*t) #postscript("../../WriteUp/Graphics/Chapter3/chap_3_sect_7_prob_29_plot.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 2)) ymin = min( c(ut, upt) ) ymax = max( c(ut, upt) ) plot(t, ut, type='l', col='blue', ylim=c(ymin, ymax), ylab='u and u prime', main='time series') lines(t, upt, type='l', col='red') legend('bottomright', c('u', 'u prime'), lty=c(1, 1), col=c('blue', 'red')) grid() plot(ut, upt, type='l', col='black', xlab='u', ylab='u prime', main='phase plot') grid() par(mfrow=c(1, 1)) #dev.off() # Problem 32: # library(deSolve) dydx = function(t, state, parameters){ with(as.list(c(state, parameters)),{ # rate of change dy = v dv = - y - eps * y^3 list(c(dy, dv)) }) } t = seq( 0, 28, length.out=1000 ) state = c( y=0, v=1 ) res = ode(y=state, times=t, fun=dydx, parms=c( eps=0.0 )) ts = res[,'time'] ys_0 = res[,'y'] res = ode(y=state, times=t, fun=dydx, parms=c( eps=0.1 )) ys_p1 = res[,'y'] res = ode(y=state, times=t, fun=dydx, parms=c( eps=0.2 )) ys_p2 = res[,'y'] res = ode(y=state, times=t, fun=dydx, parms=c( eps=0.3 )) ys_p3 = res[,'y'] #postscript("../../WriteUp/Graphics/Chapter3/chap_3_sect_7_prob_32_positive_eps_plots.eps", onefile=FALSE, horizontal=FALSE) plot( ts, ys_0, col='black', type='l', xlab='t', ylab='y(t)' ) lines( ts, ys_p1, col='blue', type='l') lines( ts, ys_p2, col='green', type='l') lines( ts, ys_p3, col='red', type='l') grid() #dev.off() state = c( y=0, v=1 ) res = ode(y=state, times=t, fun=dydx, parms=c( eps=0.0 )) ts_2 = res[,'time'] ys_0 = res[,'y'] res = ode(y=state, times=t, fun=dydx, parms=c( eps=-0.1 )) ys_p1 = res[,'y'] res = ode(y=state, times=t, fun=dydx, parms=c( eps=-0.2 )) ys_p2 = res[,'y'] res = ode(y=state, times=t, fun=dydx, parms=c( eps=-0.3 )) ys_p3 = res[,'y'] #postscript("../../WriteUp/Graphics/Chapter3/chap_3_sect_7_prob_32_negative_eps_plots.eps", onefile=FALSE, horizontal=FALSE) ymin = min( c(ys_0, ys_p1, ys_p2, ys_p3) ) ymax = max( c(ys_0, ys_p1, ys_p2, ys_p3) ) plot( ts, ys_0, col='black', type='l', xlab='t', ylab='y(t)', ylim=c(ymin, ymax) ) lines( ts, ys_p1, col='blue', type='l') lines( ts, ys_p2, col='green', type='l') lines( ts, ys_p3, col='red', type='l') grid() #dev.off()