if( !require('phaseR') ){ install.packages('phaseR') } library(phaseR) my_yprime = function(t, y, parameters) { xt = y[1] yt = y[2] dy = rep(NA, length(y)) dy[1] = yt dy[2] = -4*sin(xt) list(dy) } diff_eq_params = list() #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_3_prob_24_plot.eps', onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 2)) ts = seq(0, 5, length.out=500) res_1 = ode(y=c(0, 2), times=ts, fun=my_yprime, parms=diff_eq_params) xs_1 = res_1[, 2] # the numerical function x(t) res_2 = ode(y=c(0, 5), times=ts, fun=my_yprime, parms=diff_eq_params) xs_2 = res_2[, 2] # the numerical function x(t) Xs = cbind(xs_1, xs_2) matplot(ts, Xs, type='l', lty=1, lwd=2, xlab='t', ylab='x(t)') grid() # Part (b): Estimate the value of v_c: # vs = seq(2, 5, length.out=100) max_amplitude = c() for( v in vs ){ res = ode(y=c(0, v), times=ts, fun=my_yprime, parms=diff_eq_params) xs = res[, 2] # the numerical function x(t) max_amplitude = c(max_amplitude, max(xs) ) } plot(vs, max_amplitude, type='p', pch=19, xlab='v', ylab='maximum amplitude') grid() #dev.off() par(mfrow=c(1, 1))