# 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()