#-------------#
#    week5    #
#-------------#   


# a)

eps = 2
z0  = 1 + 0*1i
T   = 15
dt  = 0.01
tt  = seq(0,T,by=dt)

zexact = z0*exp(1i*eps*tt)

nt = length(tt)
z  = rep(0+0*1i,nt)
z[1] = z0

for(k in 2:nt)
{
  z[k] = z[k-1]*( 1 + 1i*eps*dt )
}

plot(z)
points(zexact,col="red")


# b) 

n = 100
sn = matrix(0,n,nt)
iepst = 1i*eps*tt

sn[1,] = 1 + iepst

for(k in 2:n)
{
  sn[k,] = sn[k-1,] + iepst^k/factorial(k) 
  if(k>70)
  {
    #plot(sn[k,]) 
    #readline("press enter for next plot")
  }
}

plot(sn[25,])
plot(sn[50,])
plot(sn[75,])
plot(sn[100,])


# c) 

T   = 15
dt  = 0.1
tt  = seq(0,T,by=dt)
nt = length(tt)

en = matrix(0,4,nt)
iepst = 1i*eps*tt
k = 1
for( n in c(100,1000,10000,100000) )
{
  en[k,] = ( 1 + iepst/n )^n 
  plot(en[k,])
  readline("press enter for next plot")
  k = k+1
}


# d) 

T   = 15
dt  = 0.1
tt  = seq(0,T,by=dt)
nt = length(tt)

n  = 100000
J  = rbind(c(0,1),c(-1,0))
Id = diag(c(1,1))
x = rep(0,nt)
y = rep(0,nt)

for(k in 1:nt)
{
  tk = tt[k]
  En = Id + eps*tk/n * J
  res = Id
  for(j in 1:n)
  {
    res = En %*% res
  }
  x[k] = res[1,1]
  y[k] = res[1,2]
}

plot(x,y)

k=4
points(Re(en[k,]),Im(en[k,]),col="red")


# e)

MyExp = function( A , m )
{
  n  = 2^m
  d  = length(A[1,])
  Id = diag(rep(1,d))
  B  = Id + A/n
  for(j in 1:m)
  {
    B = B %*% B
  }
  res = B
  return(res)
}

x  = rep(0,nt)
y  = rep(0,nt)
for(k in 1:nt)
{
  tk = tt[k]
  res = MyExp(eps*tk*J,20)
  x[k] = res[1,1]
  y[k] = res[1,2]
}

plot(x,y)