added exercise 3 solution

examples/lanczos/Lanczos.ipynb

@@ -116,16 +116,8 @@
    "name": "python2"
   },
   "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 2
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
    "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.8"
+   "version": "2.7.6"
   }
  },
  "nbformat": 4,

exercises/ex03_solution/TFIM_time_evolution/cpp/Makefile

+
+magnetization.png:  magnetization.txt
+	python3 plotmagn.py
+
+magnetization.txt: trFieldIsingTimeEvolution.x
+	./trFieldIsingTimeEvolution.x > magnetization.txt
+
+trFieldIsingTimeEvolution.x: trFieldIsingTimeEvolution.cpp
+	g++ -O2 -o trFieldIsingTimeEvolution.x trFieldIsingTimeEvolution.cpp
+

exercises/ex03_solution/TFIM_time_evolution/cpp/plotmagn.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+magn = np.loadtxt('magnetization.txt')
+Nt   = magn.shape[0]
+N    = magn.shape[1]
+dt   = 0.5
+
+fig = plt.figure()
+plt.pcolor(range(N+1),np.arange(Nt+1)*dt,magn)
+plt.xlim([0,N])
+plt.ylim([0,Nt*dt])
+plt.xlabel('site')
+plt.ylabel('time')
+plt.colorbar()
+plt.show()
+fig.savefig('magnetization.png')

exercises/ex03_solution/TFIM_time_evolution/cpp/trFieldIsingTimeEvolution.cpp

+#include <valarray>
+#include <vector>
+#include <iostream>
+#include <cmath>
+#include <complex>
+
+typedef std::valarray< std::complex<double> > state_vector_t;
+typedef std::vector<double> param_vector_t;
+typedef unsigned state_t;
+const std::complex<double> I(0.,1.);
+
+
+template <class Vector>
+void tstep_transv(unsigned l, const param_vector_t& cosinesPrecomputed, const Vector& sinesPrecomputed, Vector& x, Vector& tmp)
+{
+    state_t dim = 1 << l;
+    tmp *= 0;
+    for( int r = 0; r < l ; ++r )
+    {
+        // diagonal
+        for( state_t s = 0; s < dim; ++s )
+            tmp[s] =  cosinesPrecomputed[r] * x[s];
+        // off-diagonal
+        for( state_t s = 0; s < dim; ++s )
+            tmp[s^(1<<r)] += sinesPrecomputed[r] * x[s];
+        std::swap(x,tmp);
+    }
+}
+
+param_vector_t precomputeCosines(const param_vector_t &hx, double dt){
+    param_vector_t res(hx.size(),double(0.0));
+    for(int r = 0; r<hx.size(); ++r){
+        res[r] = std::cos(dt*hx[r]);
+    }
+    
+    return res;
+}
+
+template <class Vector>
+Vector precomputeSines(const param_vector_t &hx, double dt){
+    Vector res(std::complex<double>(0.0,0.0), hx.size());
+    for(int r = 0; r<hx.size(); ++r){
+        res[r] = I * std::sin(dt*hx[r]);
+    }
+    
+    return res;
+}
+
+template <class Vector>
+void tstep_diag(unsigned l, const param_vector_t &j, double dt, Vector& x)
+{
+  state_t dim = 1 << l;
+  for( state_t s = 0; s < dim; ++s )
+    {
+      double jtotal = 0.;
+      //  +J for parallel, -J for antiparallel neighbors
+      for( int r = 0; r < l - 1; ++r )
+        {
+            jtotal +=  ((s >> r)^(s >> (r+1)))&1 ? -j[r] : j[r]; // check if spins are parallel at site r and r+1
+        }        
+      x[s] *= std::exp(I*jtotal*dt);
+    }
+}
+
+template <class Vector>
+void evolve(unsigned l, const param_vector_t &j, const param_vector_t &hx, double dt, unsigned n, Vector& x)
+{
+  state_t dim = 1 << l;
+  state_vector_t tmp(x.size()); // memory for tstep_transv
+  param_vector_t cosinesPrecomputed = precomputeCosines(hx,dt);
+  Vector sinesPrecomputed = precomputeSines<Vector>(hx,dt);
+  tstep_diag(l,j,dt/2.,x);
+  for(int i=0 ; i<n-1 ; ++i){
+      tstep_transv(l,cosinesPrecomputed,sinesPrecomputed,x,tmp);
+      tstep_diag(l,j,dt,x);
+  }
+  tstep_transv(l,cosinesPrecomputed,sinesPrecomputed,x,tmp);
+  tstep_diag(l,j,dt/2.,x);
+}
+
+template <class Vector>
+void printMagnetization(unsigned l, const Vector& v, std::ostream& outputStream)
+{    
+  state_t dim = 1 << l;
+  for(int r=0; r<l; r++){
+    double m(0);
+    for(long s=0; s<dim; s++){
+      m += std::norm(v[s]) * 2. * ( (bool)(s&(1 << r)) -0.5);
+    }
+    outputStream << m << "\t\t";
+  }
+  outputStream << std::endl;
+}
+
+int main(){
+  // ======== input
+  unsigned l = 11; // length of chain
+  param_vector_t j(l-1,1.0); // Ising coupling
+  param_vector_t hx(l,0.4); // transverse magnetic field
+  double tmax = 40.; // maximum time reached in time evolutin
+  double tmeas = 0.5; // time after which a measurment is performed
+  double dt = 0.1; // time step used for the evolution
+  unsigned nmeas = tmax/tmeas; // number of measurements
+  unsigned nstep = tmeas/dt; // number of steps between measurements
+  // ======== input
+
+  // Starting configuration
+  state_t initialInd = 1<<5;
+  state_vector_t x(0., 1 << l);
+  x[initialInd] = 1;   
+  printMagnetization(l,x,std::cout);
+
+  // Do time evolution
+  for(int n=1; n<=nmeas; ++n){    
+    evolve(l, j, hx, dt, nstep, x);
+    printMagnetization(l,x,std::cout);    
+  }
+  
+  return 0;
+}

exercises/ex03_solution/TFIM_time_evolution/python/trising.py

+import numpy as np
+import matplotlib.pyplot as plt
+from cmath import exp
+from math import sin,  cos
+
+def evolve(l, j, hx, v, dt, n):
+    """ Evolve the state v in time by n * dt """
+    dim = 1 << l
+    def tstep_diag(x, dt):
+        # Apply diagonal part
+        for s in range(dim):
+            jtotal = 0.
+            for r in range(l-1):
+                jtotal += -j[r] if ((s >> r)^(s >> (r+1)))&1 else j[r]
+            x[s] = exp(-1.j*jtotal*dt)*x[s]
+        return x
+    def tstep_transv(x, dt):
+        # Apply transverse part succesively for each site
+        y =  np.zeros_like(x)
+        for r in range(l):
+            cos_v = cos(dt*hx[r])
+            sin_v = sin(dt*hx[r])
+            # diagonal
+            for s in range(dim):
+                y[s] = cos_v * x[s]
+            # off-diagonal
+            for s in range(dim):
+                y[s^(1<<r)] += 1.j*sin_v * x[s]
+            # swap vectors
+            tmp = x
+            x   = y
+            y   = tmp
+        return x
+    v = tstep_diag(v, dt/2.)
+    for i in range(int(n-1)):
+        v = tstep_transv(v, dt)
+        v = tstep_diag(v, dt)
+    v = tstep_transv(v, dt)
+    v = tstep_diag(v, dt/2.)
+    return v
+    
+def magnetization(v, l):
+    """ Measure the magnetization per site for the state v """
+    dim = 1 << l
+    m = np.zeros((l))
+    for r in range(l):
+        for s in range(dim):
+            m[r] += abs(v[s])**2 * 2. * ( bool(s&(1 << r)) -0.5)
+    return m
+
+
+l   =  9
+dim =  1 << l
+j   =  1.00*np.ones((l-1)) # Ising coupling
+hx  =  0.40*np.ones((l)) # transverse field strength Gamma
+
+# Starting configuration
+ind = int('0000001000000',base=2) # specify spin configuration: 0 = down, 1 = up
+v   = np.zeros((dim),dtype=complex)
+v[ind] = 1.
+
+# Parameters for time evolution
+tend   = 20. # duration of the time evolution
+tmeas  = 0.5 # time after which a measurment is performed
+dt     = 0.1 # time step used for the evolution
+nmeas  = int((tend)/tmeas) # number of measurements
+nstep  = int((tmeas)/dt) # number of steps between measurments
+
+# measure magnetization
+m=np.zeros((nmeas,l));
+m[0,:] = magnetization(v, l)
+
+# Do time evolution
+for t in range(1,nmeas):
+    print(t*tmeas)
+    v=evolve(l,j,hx,v,dt,nstep)
+    m[t,:] = magnetization(v, l)
+
+# plot magnetization
+plt.pcolor(np.array(range(l+1)),np.linspace(0.,tend,nmeas),m,
+           vmin=-1.0, vmax=1.0)
+plt.xlabel('site')
+plt.ylabel('time')
+plt.axis([0,l,0,tend])
+plt.title("Magnetisation, J = "+str(j[0])+" , hx = "+str(hx[0]))
+plt.colorbar()
+plt.show()

exercises/ex03_solution/TFIM_time_evolution/python/trising_skeleton.py

+# ---- CQP FS16, Ex3.1 skeleton ----
+# ----------------------------------
+import numpy as np
+import matplotlib.pyplot as plt
+
+def evolve(l, j, hx, v, dt, n):
+    """ Evolve the state v in time by n * dt 
+    
+    Parameters:
+    l  length of chain
+    j  Ising coupling for each neighbors, 1-d numpy array with length at least l-1
+    hx transverse field for each site, 1-d numpy array with length at least l
+    v  the state at time t, complex 1-d numpy array with length at least 1 << l
+    dt timestep
+    n  number of consecutive timesteps
+
+    Return:
+    v  the state at time t+n*dt
+    """
+    # Implement the time evolution of the state v here.
+    # Use the split-operator scheme we discussed in the
+    # exercise class.
+    return v
+    
+    
+def magnetization(v, l):
+    """ Measure the magnetization per site for the state v 
+    
+    Parameters:
+    v  the state to be measured, complex 1-d numpy array with length at least 1 << l
+    l  length of chain
+
+    Return:
+    m  the magnetization at each site, 1-d numpy array with length l
+    """
+    # Implement the measurement of the magnetization here.
+    return m
+    
+
+
+l   =  9 # length of chain
+dim =  1 << l
+j   =  1.00*np.ones((l)) # Ising coupling
+hx  =  0.40*np.ones((l-1)) # transverse field strength Gamma
+
+# Starting configuration
+ind = int('000010000',base=2) # specify spin configuration: 0 = down, 1 = up
+v   = np.zeros((dim),dtype=complex)
+v[ind] = 1. # start time evolution with a single basis state
+
+# Parameters for time evolution
+tend   = 10. # duration of the time evolution
+tmeas  = 0.5 # time after which a measurment is performed
+dt     = 0.1 # time step used for the evolution
+nmeas  = int((tend)/tmeas) # number of measurements
+nstep  = int((tmeas)/dt) # number of steps between measurements
+
+# measure magnetization
+m=np.zeros((nmeas,l));
+m[0,:] = magnetization(v, l)
+
+# Do time evolution
+for t in range(1,nmeas):
+    print(t*tmeas)
+    v=evolve(l,j,hx,v,dt,nstep)
+    m[t,:] = magnetization(v, l)
+
+# plot magnetization
+plt.pcolor(np.array(range(l+1)),np.linspace(0.,tend,nmeas),m,
+           vmin=-1.0, vmax=1.0)
+plt.xlabel('site')
+plt.ylabel('time')
+plt.axis([0,l,0,tend])
+plt.title("Magnetisation, J = "+str(j[0])+" , hx = "+str(hx[0]))
+plt.colorbar()
+plt.show()