Least Action

Nontrivializing triviality..and vice versa.

Archive for the ‘Computational Physics’ Category

Gnuplot with C/C++

leave a comment »

Sometimes, it is convenient to call gnuplot from within your own C/C++ code rather than having the code create a data file, and manually execute gnuplot every time to re-plot using the contents of the file. This is easy if you use Linux: just use a pipe. For instance:


#include <iostream>
#include <stdio.h>
#include <cstdlib>

using namespace std;

int main(){
// Code for gnuplot pipe
FILE *pipe = popen("gnuplot -persist", "w");
fprintf(pipe, "\n");

// Your code goes here

// Don't forget to close the pipe
fclose(pipe);
return 0;
}

Note that we have had to use C routines to call gnuplot. But there’s a problem with this code: gnuplot will generate the plot only when the pipe is closed. So if you have multiple plot statements, or if you want to refresh the same plot (as for instance, in an animation scenario) you will find that the above approach will result in the plots getting asynchronously updated right at the end.

So how does one ensure that after every plot command, a plot is actually generated (or an existing one refreshed) Immediately? The key is to flush the buffer using fflush. Simply insert


fflush(pipe);

below every plot call to gnuplot. If you have an animation, then in order to make the transitions smoother, you might want to include a delay subroutine. A simple one can be built using the clock routines in time.h. For an example, click here.

Thanks to Vidushi Sharma for bringing these issues to my attention.

Written by Vivek

June 22, 2012 at 14:13

Evaluation of Complex Error Functions erf(z) using GSL/C

with 5 comments

For a small QFT calculation, I needed to numerically evaluate the imaginary error function erfi(x) = erf(i x). But it turns out that GSL (and most other numerical recipe code I could find) can only deal with erf(x), where x is real. Here’s a poor man’s implementation of erf(z) through a standard Taylor expansion,

erf(z) = \frac{2}{\sqrt{\pi}}\sum_{n = 0}^{\infty}\frac{(-1)^n z^{2n+1}}{n!(2n+1)}

The catch here is to deal with propagation of errors in the complex Taylor series, and to also somehow benchmark the results. Well, I haven’t been able to think about this yet, but I was able to confirm that for real arguments, my Taylor series code is about as good as the GSL error function gsl_sf_erf(x).

I am working on a CUDA implementation of this now, because in my project, I need to perform a numerical integration over the error function, which is quite intensive even for Mathematica. For now, I’m just sharing the serial implementation. (Hint: if you can write wrapper for each of the gsl functions inside my Taylor series calculating method, you can embed them within a __global__ kernel call through a struct. Alternatively — and less fun — you can just parallelize the for loop.)

/* Function to compute erf(z) using a Taylor series expansion
/* Author: Vivek Saxena
/* Last updated: January 28, 2011 21:11 hrs */
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
#include <stdio.h>

#define PI 3.1415926543
double cz = 2/sqrt(PI);
const int TERMS = 10; // no of terms to use in the Taylor series

gsl_complex erfTaylor( gsl_complex z, int trunc ){
  gsl_complex res = gsl_complex_rect(0,0),
  num  = gsl_complex_rect(0,0),
  den  = gsl_complex_rect(1,0),
  snum = gsl_complex_rect(1,0),
  temp = gsl_complex_rect(0,0);
  for(int i = 0; i < trunc; i++){
     snum = gsl_complex_rect( cz * gsl_sf_pow_int(-1, i), 0 );
     num  = gsl_complex_mul(snum, gsl_complex_pow_real(z, 2*i+1));
     den  = gsl_complex_rect((2*i + 1)*gsl_sf_fact(i),0);
     temp = gsl_complex_div(num, den);
     res  = gsl_complex_add(res, temp);
  }

 return res;
}

int main ( void ){
   printf( "Real error function\n\n");
   for ( float i = 0; i <= 1; i += 0.01 ){
     float gslerror = gsl_sf_erf(i);
     float taylor   = GSL_REAL( erfTaylor( gsl_complex_rect(i, 0), 10 ) );
     printf("erf(%f): gsl = %f, taylor = %f, mag error = %f\n", i, gslerror, taylor, abs(gslerror-taylor));
   }

   gsl_complex t, arg;
   printf( "\n\nImaginary error function\n\n");
   for (float i = 0; i <= 1; i += 0.01 ){
      /* this would be your generic argument z in erf(z).
       * so if z = x + iy, then
       * arg = gsl_complex_rect(x, y);
      */
      arg = gsl_complex_rect(0, i);
      t   = erfTaylor( arg, TERMS );
      printf("erf(%f + i %f) = %f + i %f\n", GSL_REAL(arg), GSL_IMAG(arg), GSL_REAL(t), GSL_IMAG(t));
   }
   return 0;
}

Written by Vivek

January 28, 2011 at 21:12

Cadabra – a Computer Algebra System for field theory

leave a comment »

This seems very exciting:

Cadabra is a computer algebra system (CAS) designed specifically for the solution of problems encountered in field theory. It has extensive functionality for tensor computer algebra, tensor polynomial simplification including multi-term symmetries, fermions and anti-commuting variables, Clifford algebras and Fierz transformations, implicit coordinate dependence, multiple index types and many more. The input format is a subset of TeX. Both a command-line and a graphical interface are available.

Source: http://cadabra.phi-sci.com/

There are two interesting papers on this. The first is a semi-technical overview, and the other (hep-th/0701238) is a more comprehensive one geared towards an audience familiar with various problems in modern field theory. The abstract of the second paper reads:

Introducing Cadabra: a symbolic computer algebra system for field theory problems

(Submitted on 25 Jan 2007 (v1), last revised 14 Jun 2007 (this version, v2))

Abstract: Cadabra is a new computer algebra system designed specifically for the solution of problems encountered in field theory. It has extensive functionality for tensor polynomial simplification taking care of Bianchi and Schouten identities, for fermions and anti-commuting variables, Clifford algebras and Fierz transformations, implicit coordinate dependence, multiple index types and many other field theory related concepts. The input format is a subset of TeX and thus easy to learn. Both a command-line and a graphical interface are available. The present paper is an introduction to the program using several concrete problems from gravity, supergravity and quantum field theory.

Source = http://arxiv.org/abs/hep-th/0701238

Written by Vivek

November 8, 2010 at 00:35