#-----------------------#
#   Blatt8, Aufgabe 1   #
#-----------------------#

# R-Code:

alpha = 4
lam = 1
T = 30

dt = 0.0005
tt = seq(from=0, to=T, by=dt)
nt = length(tt)

x1 = rep(0,nt)
p1 = rep(0,nt)
x2 = rep(0,nt)
p2 = rep(0,nt)
x3 = rep(0,nt)
p3 = rep(0,nt)
x1[1] = 0.5
x2[1] = 1 
x3[1] = 2

for(k in 2:nt)
{
  x1[k] = x1[k-1] + p1[k-1]*dt
  x2[k] = x2[k-1] + p2[k-1]*dt
  x3[k] = x3[k-1] + p3[k-1]*dt
  p1[k] = p1[k-1] - lam * sign(x1[k-1]) * abs(x1[k-1])^(alpha-1) * dt
  p2[k] = p2[k-1] - lam * sign(x2[k-1]) * abs(x2[k-1])^(alpha-1) * dt
  p3[k] = p3[k-1] - lam * sign(x3[k-1]) * abs(x3[k-1])^(alpha-1) * dt
}

info = "x0 = 1 in black,   x0 = 1/2 in green,   x0 = 2 in blue 
   alpha = 4,  quartic oscillator"

plot(tt,x3,col="blue",main=info,xlab="time t",ylab="x(t)")
points(tt,x1,col="green")
points(tt,x2,col="black")