mu_given_gamma = function(gamma){ sqrt(4 - gamma^2)/2 } u_fn = function(t, gamma){ mu = mu_given_gamma(gamma) exp(-(gamma*t)/2) * ( 2 * cos( mu * t ) + (gamma/mu) * sin( mu * t ) ) } compute_tau = function(gamma, t_min=0, t_max=75, y_thresh=0.01){ # # Computes the time tau where |u(t)| < y_thresh for all t > tau # mu = mu_given_gamma(gamma) t = seq( t_min, t_max, length.out=50000 ) y = exp(-(gamma*t)/2) * ( 2 * cos( mu * t ) + (gamma/mu) * sin( mu * t ) ) # Walk backwards until we get the first FALSE: # below_t = abs(y) < y_thresh t_index = length(y) while( below_t[t_index] ){ #print(c(gamma, t_index, below_t[t_index], y[t_index], y_thresh)) t_index = t_index-1 } tau = t[t_index+1] list(t=t, y=y, tau=tau) } approx_tau_fn = function(gamma){ # # An approximation to tau see Part (e) # (2/gamma)*log(400/sqrt(4-gamma^2)) } # Part (a): # y_thresh = 0.01 gamma = 0.25 ct_res = compute_tau(gamma) tau = ct_res\$tau print(sprintf('tau= %.3f', tau)) tau = approx_tau_fn(gamma) print(sprintf('gamma= %.2f; approximate tau= %.2f', gamma, tau)) #postscript("../../WriteUp/Graphics/Chapter3/sect_7_prob_25_plots.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 3)) ts = seq(0, 50, length.out=500) ys = u_fn(ts, gamma) plot( ts, ys, type='l', main=sprintf('Problem 25: gamma: %.2f; tau= %.2f', gamma, tau), xlab='t', ylab='y(t)') abline(v=tau, col='black') abline(h=c(-y_thresh, y_thresh), col='blue') grid() # Part (b): # gammas = seq(0.5, 1.5, length.out=100) taus = sapply(gammas, function(x){ compute_tau(x, gamma)\$tau }) plot(gammas, taus, type='b', xlab='gamma', ylab='tau', main=sprintf('%.2f < gamma < %.2f', min(gammas), max(gammas))) lines(gammas, approx_tau_fn(gammas), type='l', col='red') grid() # Part (c): # gammas = seq(1.5, 1.99, length.out=100) taus = sapply(gammas, function(x){ compute_tau(x, gamma)\$tau }) min_tau_indx = which.min(taus) tau0 = taus[min_tau_indx] gamma0 = gammas[min_tau_indx] ts = sprintf('gamma0= %.3f; tau0= %.3f', gamma0, tau0) print(ts) plot(gammas, taus, type='b', xlab='gamma', ylab='tau', main=sprintf('%.2f < gamma < %.2f', min(gammas), max(gammas))) lines(gammas, approx_tau_fn(gammas), type='l', col='red') abline(v=gammas[min_tau_indx], col='black') grid() par(mfrow=c(1, 1)) #dev.off()