#------------------------------------#
#                                    #
#    Week3: Rechnen mit Matrizen,    #
#              Teil 2                #
#------------------------------------#


# 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.001
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")



# 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])
points(tt[kplot],xt_exact[kplot],col="red")