#-----------------#
#      week3      #
#-----------------#


# a)

eps = 2
x0  = 1
T   = 50

dt = 0.01
tt = seq(from=0,to=T,by=dt)
xt = x0*cos(eps*tt)

plot(tt,xt)


# b)

T  = 50
dt = 0.01
tt = seq(from=0,to=T,by=dt)
M  = length(tt)

xt = rep(0,M)
vt = rep(0,M)
xt[1] = x0
vt[1] = 0

t0 = Sys.time()            # wir messen die Rechenzeit
for(k in 2:M)
{
  xt[k] = xt[k-1] + vt[k-1]*dt
  vt[k] = vt[k-1] - eps^2*xt[k-1]*dt
}
Sys.time() - t0

xt_exact = x0*cos(eps*tt)
plot(tt,xt)
points(tt,xt_exact,col="red")


# c)

y = c(x0,0)
xt = rep(0,M)
xt[1] = x0
 
r1 = c(1,dt)
r2 = c(-eps^2*dt,1)
R = rbind(r1,r2)
R
for(k in 2:M)
{
  y = R %*% y
  xt[k] = y[1] 
}

xt_exact = x0*cos(eps*tt)
plot(tt,xt)
points(tt,xt_exact,col="red")


# d)

T  = 50
dt = 0.01
tt = seq(from=0,to=T,by=dt)
M  = length(tt)

A = rbind( c(0,1) , c(-eps^2,0) )
Id = diag(rep(1,2))
Id
A2 = A %*% A
A3 = A2 %*% A
A4 = A2 %*% A2
expdtA2 = Id + dt*A + dt^2/2 * A2
expdtA3 = expdtA2 + dt^3/6 * A3
expdtA4 = expdtA3 + dt^4/24 * A4


y = c(x0,0)
xt = rep(0,M)
xt[1] = x0
 
for(k in 2:M)
{
  y = expdtA3 %*% y      # das ist also numerisch deutlich stabiler
  xt[k] = y[1] 
}

xt_exact = x0*cos(eps*tt)
plot(tt,xt)
points(tt,xt_exact,col="red")



# e)

y = c(x0,0)
xt = rep(0,M)
xt[1] = x0

u1 = c( cos(eps*dt) , sin(eps*dt)/eps )
u2 = c( -eps*sin(eps*dt) , cos(eps*dt) )
U = rbind( u1 , u2 )


for(k in 2:M)
{
  y = U %*% y
  xt[k] = y[1] 
}

xt_exact = x0*cos(eps*tt)
plot(tt,xt)
points(tt,xt_exact,col="red")