#-----------------------#
#   Loesungen Blatt 4   #
#-----------------------#


# a)

eps = 2
omega = 2.1
x0 = 1
a0 = 0.5
T = 200

dt = 0.01
tt = seq(from=0,to=T,by=dt)
term1 = x0*cos(eps*tt)
term2 = a0/(eps^2-omega^2) * ( sin(omega*tt) - omega/eps*sin(eps*tt) )
xt_exact = term1 + term2

plot(tt,xt_exact)


# b)

eps = 2
omega = 2.1
x0 = 1
a0 = 0.5
T = 200

dt = 0.0001
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()
for(k in 2:M)
{
  xt[k] = xt[k-1] + vt[k-1]*dt
  vt[k] = vt[k-1] + ( a0*sin(omega*(k-1)*dt) - eps^2*xt[k-1] )*dt
}
Sys.time() - t0

plot(tt,xt)

term1 = x0*cos(eps*tt)
term2 = a0/(eps^2-omega^2) * ( sin(omega*tt) - omega/eps*sin(eps*tt) )
xt_exact = term1 + term2
points(tt,xt_exact,col="red")



# Wenn die Anzahl der zu plottenden Punkte sehr gross ist, so ab 100000, 
# weil man das dt sehr klein gewaehlt hat, kann das manchmal ein bischen 
# dauern. Fuer die Optik reichen auch 10000 Punkte, kann man dann so 
# machen:

nt = length(tt)
kplot = seq(from=1,to=nt,length=10000)            # also 10000 Stueck
plot(tt[kplot],xt[kplot])

term1 = x0*cos(eps*tt)
term2 = a0/(eps^2-omega^2) * ( sin(omega*tt) - omega/eps*sin(eps*tt) )
xt_exact = term1 + term2
points(tt[kplot],xt_exact[kplot],col="red")