Code-Bin #879

untitled C++ Code
Posted by: Billy | October 6, 2010 @ 7:04pm

C++ Code

#include <iostream>
//include for random numbers
#include <cstdlib>
//include for pow function
#include <cmath>
//include to write to file
#include <fstream>
//include to get the time for seed
#include <time.h>

using namespace std;

//number of time steps N, total time T = Na
const unsigned int N=400;

//time step a
const float a=0.25;

const int m=1;
const int w=1;
const int hbar=1;
const float delta=1.5;

float du, k, g;
float accepted = 0;

//generates a random number between -k and k
float randomNumber(int k) {
  return (rand()/(float)RAND_MAX)*2*k -k;
}

//method to write to file
int writeOutputToFile(char * filename, float x[]) {
  ofstream f;
  f.open(filename);
  for(int i=0; i<N+1; i++) {
    f << x[i] << " " << i << endl;
  }
  f.close();
  return 0;
}

//this infact computes partial S/a. Leave x a till the end to improve speed slightly?
/*
float partialS(float xi1, float xi) {
  float s;
  s = (m/2)* pow( (xi1 - xi)/a, 2 ) + 0.5 * m * pow(w*(xi1+xi)/2, 2);
  return s;
}

//calculate the total action
float totalS(float x[]) {
  float action = 0;

//loop through all sites to get the total action by summing each partial action
  for (int i=0; i<N+1; i++) {
    action = action + partialS(x[i+1], x[i]);
  }
  
  action = a * action;
  return action;
}
*/

//calculate the change in the action/hbar
//x1 = u_{i-1}, x2= u_{i}, x3 = u_{i+1}
float dAction(float x1, float x2, float x3) {
  return du*(du + 2*x2 - g*(x1 + x3));
}

//metropolis algorithm
int metropolis(float x[]) {
  float newx, newPartialS, oldPartialS, r;
  float dS;
  //I initially had i<N+1 but this changes x[0] = x[N+1] and I need end points fixed
  
  for(int i=1; i<N; i++) {
   
    //newx is the oldx +- randomNumber between delta and -delta
    newx = x[i] + randomNumber(delta);

//calculate the change in action/hbar
    dS = dAction(x[i-1], x[i], x[i+1]);

//if the action is reduced with the new x, then set x[i] = newx;
    if (dS<0) {
      x[i] = newx;
      accepted++;
      //otherwise...
    } else {

//create a random number between 0..1 with uniform probability
      r = rand()/(float)RAND_MAX;
    
      
      //check to see if we accept newx based on the random number generated and the change in action
      if (exp(-1 * dS ) > r) {
	//if true, ACCEPT!
	x[i] = newx;
	accepted++;
      } else {
	//cout << "here" << endl;

}
    //cout << (accepted/(float)N) << endl;
    accepted = 0;
  }
  
  return 0;
  
}

float monteCarlo(float x[]) {
  float integral =0;

//number of monte carlo interations
  int iter = 300;
  writeOutputToFile("output1.txt", x);
  //300 monte carlo iterations
  for(int i=0; i<iter; i++) {

//calculate new path/configuration, a few times
    for(int j=0; j<50; j++) {
      metropolis(x);
    }

//compute correlation function
    integral = integral + x[N-1] * x[1];

}
  
  writeOutputToFile("output.txt", x);
  return integral/(float)iter;
 
}

//the main method!
int main() {

//set up the variables a little

k = 0.25 * a*a * w*w;
  g = (1-k)/(1+k);

//using float because RAND_MAX is only 10^9, no point in using double
  //create an array of size N, since j=0..N-1, the Markov chain
  float u[N+1];

//set the seed for rand()
  srand(time(0));

//inputs random variables into sites, being between -5 and 5
  u[0] = 0;
  u[N] = 0;
  for(int i=1; i<N; i++) {
    u[i] = randomNumber(5);
  }

//monteCarlo(u);
  //cout <<  monteCarlo(u) << endl;
 
  return 0;

}

Syntax Highlighting