#-------------#
#    week4    #
#-------------#   


# a)

eps = 2
mu  = 0.1
x0  = 1
v0  = 0
T   = 50
om  = sqrt(eps^2-mu^2)

dt  = 0.01
tt  = seq(0,T,by=dt)

xexact = x0*exp(-mu*tt)*( cos(om*tt) + mu/om*sin(om*tt) )

plot(tt,xexact)


# b)

nt = length(tt)
x  = rep(0,nt)
v  = rep(0,nt)
x[1] = x0
v[1] = v0

for(k in 2:nt)
{
  x[k] = x[k-1] + v[k-1]*dt
  v[k] = v[k-1] -( eps^2*x[k-1] + 2*mu*v[k-1] )*dt
}

plot(tt,x)
points(tt,xexact,col="red")


# c)

Id = diag(c(1,1))
A = rbind( c(0,1) , c(-eps^2,-2*mu) )
R = Id + A*dt

y = c(x0,v0)
x[1] = y[1]
for(k in 2:nt)
{
  y = R %*% y
  x[k] = y[1]
}

plot(tt,x)
points(tt,xexact,col="red")


# c2)

Id = diag(c(1,1))
A = rbind( c(0,1) , c(-eps^2,-2*mu) )
R = Id + A*dt

y = matrix(0,2,nt)
y[,1] = c(x0,v0)

for(k in 2:nt)
{
  y[,k] = R %*% y[,k-1]
}

plot(tt,y[1,])
points(tt,xexact,col="red")


# d)  Matrix-Exponential machen wir im week5