Homework "plots"

Prerequisites: install pyxplot and/or gnuplot and/or plotutils on your box:

• Pyxplot:

  sudo apt install pyxplot
  

  Alternatively, you can install the latest version of Pyxplot directly from its github repository. First, install the needed libraries,

  sudo apt-get install libfftw3-dev gcc ghostscript gv imagemagick \
         libc6-dev libcfitsio-dev libgsl-dev libkpathsea-dev \
         libpng-dev libreadline-dev libxml2-dev make \
         texlive-latex-extra texlive-latex-recommended \
         texlive-fonts-extra texlive-fonts-recommended wget \
         zlib1g-dev \
  

  (perhaps all these "extra" LaTeX things (there are about a gigabyte of them) are not strictly needed.)

  Then install Pyxplot:

  mkdir ~/local
  cd ~/local
  git clone https://github.com/dcf21/pyxplot9.git
  cd pyxplot9
  ./configure # and install missing libraries, if anything is missing
  make -j4    # or whatever number of processors you've got
  make -j4 test
  sudo make install
  

• Gnuplot:

  sudo apt --yes install gnuplot
  

  MacOS with homebrew:

  brew install gnuplot [--with-qt|--with-x11] [--with-cairo]
  

• Plotutils:

  sudo apt --yes install plotutils
  

The homework:

1. Plot the error-function together with several of its tabulated values (as a test). You can use the following approximation,

   static double erf(double x){
   /// single precision error function (Abramowitz and Stegun, from Wikipedia)
   if(x<0) return -erf(-x);
   double[] a={0.254829592,-0.284496736,1.421413741,-1.453152027,1.061405429};
   double t=1/(1+0.3275911*x);
   double sum=t*(a[0]+t*(a[1]+t*(a[2]+t*(a[3]+t*a[4]))));/* the right thing */
   return 1-sum*Exp(-x*x);
   } 

   It should produce a curve passing through the tabulated points.

2. Plot the gamma-function together with several of its tabulated values (factorials) as a test. Try to reproduce the plot from the Wikipedia article. You can use the following [Stirling approximation],

   static double gamma(double x){
   ///single precision gamma function (Gergo Nemes, from Wikipedia)
   if(x<0)return PI/sin(PI*x)/gamma(1-x);
   if(x<9)return gamma(x+1)/x;
   double lngamma=x*Log(x+1/(12*x-1/x/10))-x+Log(2*PI/x)/2;
   return Exp(lngamma);
   }

   The gamma-function overflows very easily, so the logarithm of the gamma function, lngamma is often a more useful function. Figure out how to modify the above formula to calculate lngamma. Make a similar plot as above this time for the logarithm of the gamma-function.

   Please not that this is an asymptotic formula so that the shift of the argument to larger values is indeed necessary.

3. Implement the complex function complex G(complex z) that calculates Γ(z) for complex arguments: use the (suitable modified) formula from wikipedia above, it works for complex arguments as well, I believe. Try to reproduce the 3d plot of the absolute value of the Γ-function in the complex plane from Wikipedia article.