added ex03 solution

exercises/ex03_solution/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/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/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;
+}