g = 9.8 m = 0.15 y0 = 30 y0p = 20 t = seq( 0, 5.5, length.out=1000 ) # Position and velocity for problem 20: # v_20 = - g * t + y0p y_20 = -(g/2) * t^2 + y0p * t + y0 # Position and velocity for problem 21: # # t_21 = t v_21 = -30 * m * g + ( y0p + 30 * m * g ) * exp( -t/(30*m) ) y_21 = -30 * m * g * t + 30 * m * ( y0p + 30 * m * g ) * ( 1 - exp( -t/(30*m) ) ) + y0 y_21_max_height = max(y_21) print(sprintf('Prob 21: maximum height %f', y_21_max_height)) mi = which.min(abs(y_21)) t_hit_ground = t_21[mi] print( sprintf('Prob 21: t when hit ground %f', t_hit_ground) ) # Position and velocity for problem 22: # library(deSolve) parameters = c( m=m, g=g ) state = c( y=y0, v=y0p ) dydx = function(t, state, parameters){ with(as.list(c(state, parameters)),{ # rate of change dy = v dv = -g - (sign(v)*v^2)/(1325*m) list(c(dy, dv)) }) } res = ode(y=state, times=t, fun=dydx, parms=parameters) t_22 = res[,'time'] y_22 = res[,'y'] v_22 = res[,'v'] y_22_max_height = max(y_22) print(sprintf('Prob 22: maximum height %f', y_22_max_height)) mi = which.min(abs(y_22)) t_hit_ground = t_22[mi] print( sprintf('Prob 22: t when hit ground %f', t_hit_ground) ) #postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_3_probs_20_22_plot.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1,2)) plot( t, y_20, col='black', type='l', xlab='t', ylab='y(t)' ) lines( t_21, y_21, col='blue' ) lines( t_22, y_22, col='green' ) grid() plot( t, v_20, type='l', xlab='t', ylab='v(t)' ) lines( t_21, v_21, col='blue' ) lines( t_22, v_22, col='green' ) grid() legend( 'topright', c('No drag', '|v| drag', 'v^2 drag'), col=c('black', 'blue', 'green'), lty=1, lwd=2 ) par(mfrow=c(1,1)) #dev.off()