examples/sf/erfc_check

   1 #!/usr/bin/perl
   2 # interfc
   3 #   v0.2   Keith Lofstrom   2008 Sept 18
   4 #         with tweaks by Jonathan Leto
   5 #
   6 # Compare the results of numerical integration of
   7 # (2.0/sqrt(PI)*integral[x,inf](e^(-x^2) ) with the
   8 # erfc() function.
   9 #
  10 # The results are pretty good, within 1e-14 or better
  11 # for the examples tried
  12
  13 use Math::GSL::SF qw/:all/;
  14 use Math::GSL::Integration qw/:all/;
  15 use Math::GSL::Const qw/:all/;
  16 use strict;
  17
  18 my $workspace  = gsl_integration_workspace_alloc(101);
  19 my $maxdiff = 0;
  20
  21 printf " arg         erfc             int()             difference\n";
  22
  23 for( my $i=-5.2 ; $i<=5.201 ; $i += 0.1 ) {
  24    my $e  = gsl_sf_erfc($i);
  25    my ($status, $result, $abserr ) = gsl_integration_qagiu(
  26                      sub { exp( -($_[0]**2) ) },
  27                      $i, 0, 1e-9, 100, $workspace );
  28
  29    my $yy = $M_2_SQRTPI * $result ;
  30    my $diff=$yy-$e ;
  31    $maxdiff = abs($diff) > $maxdiff ? abs($diff) : $maxdiff;
  32
  33    printf "%4.1f%18.12f%18.12f %.18g\n", $i, $e, $yy, $diff ;
  34 }
  35 printf "Maximum local error: %.12g\n", $maxdiff;