update week 9

doc/src/week9/programs/automersenne.cpp

@@ -0,0 +1,76 @@
+//  This function computes the autocorrelation function for 
+//  the standard c++ random number generator
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <cmath>
+#include <random>
+
+using namespace std;
+// output file as global variable
+ofstream ofile;  
+
+//     Main function begins here     
+int main(int argc, char* argv[])
+{
+     int n;
+     char *outfilename;
+
+     cin >> n;
+     double MCint = 0.;      double MCintsqr2=0.;
+  // Initialize the seed and call the Mersienne algo
+     std::random_device rd;
+     std::mt19937_64 gen(rd());
+  // Set up the uniform distribution for x \in [[0, 1]
+     std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
+     // Compute the variance and the mean value of the uniform distribution
+     // Compute also the specific values x for each cycle in order to be able to
+     // the covariance and the correlation function  
+     // Read in output file, abort if there are too few command-line arguments
+     if( argc <= 2 ){
+       cout << "Bad Usage: " << argv[0] << 
+	 " read also output file and number of cycles on same line" << endl;
+       exit(1);
+     }
+     else{
+       outfilename=argv[1];
+     }
+     ofile.open(outfilename); 
+     // Get  the number of Monte-Carlo samples
+     n = atoi(argv[2]);
+     double *X;  
+     X = new double[n];
+     for (int i = 0;  i < n; i++){
+           double x = RandomNumberGenerator(gen);
+           X[i] = x;
+           MCint += x;
+           MCintsqr2 += x*x;
+     }
+     double Mean = MCint/((double) n );
+     MCintsqr2 = MCintsqr2/((double) n );
+     double STDev = sqrt(MCintsqr2-Mean*Mean);
+     double Variance = MCintsqr2-Mean*Mean;
+//   Write mean value and standard deviation 
+     cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl;
+
+     // Now we compute the autocorrelation function
+     double *autocor;  autocor = new double[n];
+     for (int j = 0; j < n; j++){
+       double sum = 0.0;
+       for (int k = 0; k < (n-j); k++){
+	 sum  += (X[k]-Mean)*(X[k+j]-Mean); 
+       }
+       autocor[j] = sum/Variance/((double) n );
+       ofile << setiosflags(ios::showpoint | ios::uppercase);
+       ofile << setw(15) << setprecision(8) << j;
+       ofile << setw(15) << setprecision(8) << autocor[j] << endl;
+     }
+     ofile.close();  // close output file
+     return 0;
+}  // end of main program 
+
+
+
+
+

doc/src/week9/programs/mc.py

@@ -0,0 +1,54 @@
+from  matplotlib import pyplot as plt
+from math import exp, acos, log10
+import numpy as np
+from sympy import Symbol, integrate, exp, oo
+import random
+
+
+# function for the trapezoidal rule
+def TrapezoidalRule(a,b,f,n):
+   h = (b-a)/float(n)
+   s = 0
+   x = a
+   for i in range(1,n,1):
+       x = x+h
+       s = s+ f(x)
+   s = 0.5*(f(a)+f(b))+s
+   return h*s
+# function to perform the Monte Carlo calculations
+def MonteCarloIntegration(f,n):
+    sum = 0
+# Define the seed for the rng
+    random.seed()    
+    for i in range (1, n, 1):
+        x = random.random()
+        sum = sum +f(x)
+    return sum/n
+
+#  function to compute
+def function(x):
+    return 4/(1+x*x)
+
+# Integration limits for the Trapezoidal rule
+a = 0.0; b = 1.0
+# define x as a symbol to be used by sympy
+x = Symbol('x')
+# find result from sympy
+#exact = integrate(function(x), (x, a, b))
+exact = acos(-1.0)
+# set up the arrays for plotting the relative error
+log10n = np.zeros(6); Trapez = np.zeros(6); MCint = np.zeros(6);
+# find the relative error as function of integration points
+for i in range(1, 6):
+    npts = 10**(i+1)
+    log10n[i] = log10(npts)
+    Trapez[i] = log10(abs((TrapezoidalRule(a,b,function,npts)-exact)/exact))
+    MCint[i] = log10(abs((MonteCarloIntegration(function,npts)-exact)/exact))
+plt.plot(log10n, Trapez ,'b-',log10n, MCint,'g-')
+plt.axis([1,6,-14.0, 0.0])
+plt.xlabel('$\log_{10}(n)$')
+plt.ylabel('Relative error')
+plt.title('Relative errors for Monte Carlo integration and Trapezoidal rule')
+plt.legend(['Trapezoidal rule', 'Brute force Monte Carlo integration'], loc='best') 
+plt.savefig('mcintegration.pdf')
+plt.show()

doc/src/week9/programs/plot.py

@@ -0,0 +1,16 @@
+import numpy as np
+from  matplotlib import pyplot as plt
+# Load in data file
+data = np.loadtxt("autocor.dat")
+data1 = np.loadtxt("automersenne.dat")
+# Make arrays containing x-axis and binding energies as function of A
+x = data[:,0]
+corr = data[:,1]
+corr2 = data1[:,1]
+plt.plot(x, corr ,'ro', x, corr2, 'b')
+plt.axis([0,1000,-0.2, 1.1])
+plt.xlabel(r'$d$')
+plt.ylabel(r'$C_d$')
+plt.title(r'autocorrelation function for RNG')
+plt.savefig('autocorr.pdf')
+plt.show()

doc/src/week9/programs/uniformhisto.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python
+import numpy as np
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+import random
+
+# initialize the rng with a seed
+random.seed() 
+counts = 10000
+values = np.zeros(counts)   
+for i in range (1, counts, 1):
+    values[i] = random.random()
+
+# the histogram of the data
+n, bins, patches = plt.hist(values, 10, facecolor='green')
+
+plt.xlabel('$x$')
+plt.ylabel('Number of counts')
+plt.title(r'Test of uniform distribution')
+plt.axis([0, 1, 0, 1100])
+plt.grid(True)
+
+plt.show()