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) - parameters$gamma * yt list(dy) } diff_eq_params = list(gamma=1.0) #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_3_prob_25_plot.eps', onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 3)) ts = seq(0, 20, length.out=500) res_1 = ode(y=c(0, 2), times=ts, fun=my_yprime, parms=list(gamma=1/4) ) xs_1 = res_1[, 2] # the numerical function x(t) res_2 = ode(y=c(0, 5), times=ts, fun=my_yprime, parms=list(gamma=1/4)) 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) mean_level = c() for( v in vs ){ res = ode(y=c(0, v), times=ts, fun=my_yprime, parms=list(gamma=1/4)) xs = res[, 2] # the numerical function x(t) mean_level = c(mean_level, mean(xs) ) } plot(vs, mean_level, type='p', pch=19, xlab='v', ylab='maximum amplitude') grid() # Part (c): # estimate_vc = function(gamma=1/4) { # Estimate the value of v_c using logic like in the plots above # vs = seq(2, 5, length.out=100) mean_level = c() for( v in vs ){ res = ode(y=c(0,v), times=ts, fun=my_yprime, parms=list(gamma=gamma ) ) xs = res[, 2] # the numerical function x(t) mean_level = c(mean_level, mean(xs)) } largest_jump_idx = which(diff(mean_level) == max(diff(mean_level))) return(vs[largest_jump_idx]) } print(sprintf('for gamma= %f v_c= %f', 1/4, estimate_vc(1/4))) # Use the above function to estimate v_c for several values of gamma: # gamma_list = seq(0.25, 0.8, length.out=50) vcs = c() for( g in gamma_list ){ vc_est = estimate_vc(g) #print(sprintf('for gamma= %f v_c= %f', g, vc_est)) vcs = c(vcs, vc_est) } plot(gamma_list, vcs, type='p', pch=19, xlab='gamma', ylab='v_c') grid() #dev.off() par(mfrow=c(1, 1))