t/Integration.t

   1 package Math::GSL::Integration::Test;
   2 use Math::GSL::Test qw/:all/;
   3 use base 'Test::Class';
   4 use Test::More 'no_plan';
   5 use Math::GSL qw/:all/;
   6 use Math::GSL::Integration qw/:all/;
   7 use Math::GSL::Errno qw/:all/;
   8 use Test::Exception;
   9 use Math::GSL::Const qw/:all/;
  10 use Data::Dumper;
  11 use strict;
  12 use warnings;
  13
  14 BEGIN{ gsl_set_error_handler_off() };
  15
  16 sub make_fixture : Test(setup) {
  17     my $self = shift;
  18     $self->{wspace} = gsl_integration_workspace_alloc(1000);
  19 }
  20
  21 sub teardown : Test(teardown) {
  22     my $self = shift;
  23     gsl_integration_workspace_free($self->{wspace});
  24 }
  25
  26 sub TEST_QAG : Tests {
  27     my $self = shift;
  28     my ($status, $result, $abserr) = gsl_integration_qag (
  29                                 sub { $_[0]**2} ,
  30                                 0, 1, 0, 1e-7, 1000,
  31                                 6, $self->{wspace},
  32                            );
  33     ok_status($status);
  34     my $res = abs($result - 1/3);
  35     ok_similar(
  36         [$result], [1/3],
  37         sprintf('gsl_integration_qag: res=%.18f, abserr=%.18f',$res,$abserr),
  38         $abserr
  39     );
  40 }
  41
  42 sub verify_integral {
  43     my ($integrator,$func,$actual,$epsabs,$epsrel,$lower,$upper,$key) = @_;
  44     my $wspace = gsl_integration_workspace_alloc(1000);
  45     my $format = "$integrator: actual=%.9f res=%.18f abserr=%.18f";
  46     my ($status, $result, $abserr, @params);
  47     if ( $integrator eq 'gsl_integration_qag' ) {
  48          push @params, $func, $lower, $upper, $epsabs, $epsrel, 1000, $key, $wspace;
  49     } elsif ($integrator eq 'gsl_integration_qags' ) {
  50          push @params, $func, $lower, $upper, $epsabs, $epsrel, 1000, $wspace;
  51     } elsif ($integrator eq 'gsl_integration_qagi' ) {
  52          push @params, $func, $epsabs, $epsrel, 1000, $wspace;
  53     } elsif ($integrator eq 'gsl_integration_qagiu' ) {
  54          push @params, $func, $lower, $epsabs, $epsrel, 1000, $wspace;
  55     }
  56     { no strict 'refs'; ($status, $result, $abserr) = $integrator->( @params ) }
  57     ok_status($status);
  58     ok_similar( [$result], [$actual],
  59         sprintf($format,abs($result-$actual),$actual,$abserr), $abserr
  60     );
  61 }
  62
  63 sub TEST_QAGI : Tests {
  64     my $integrator = 'gsl_integration_qagi';
  65     my $func = sub { my $x=shift; exp( -$x **2) };
  66     verify_integral($integrator, $func, sqrt($M_PI), 0, 1e-7);
  67     $func = sub { my $x = shift; 1/($x**2 + 1) };
  68     verify_integral($integrator, $func, $M_PI, 0, 1e-7);
  69 }
  70
  71 sub TEST_QAGIU : Tests {
  72     my $self = shift;
  73     my ($status, $result, $abserr) = gsl_integration_qagiu (
  74                                 sub { my $x=shift; log($x)/(1+100+$x**2) } ,
  75                                  0.0, 0.0, 1.0E-3, 1000,
  76                                 $self->{wspace}
  77                            );
  78    ok_status($status);
  79    local $TODO = "Testing with the data from gsl tests, don't know where's the error";
  80 #   ok_similar_relative($result, -3.616892186127022568E-01, "gsl_integration_qagiu",1e-14);
  81 #   ok_similar_relative($abserr, 3.016716913328831851E-06, "gsl_integration_qagiu absolute error",1e-5);
  82
  83
  84
  85
  86 #    my $integrator = 'gsl_integration_qagiu';
  87 #    my $func = sub { my $x=shift; log($x)/(1+100+$x**2) };
  88 #    verify_integral($integrator, $func, -3.616892186127022568E-01, 0, 1.0e-3, 0);
  89 }
  90
  91 sub TEST_QAG2 : Tests {
  92     my $self = shift;
  93     my ($status, $result, $abserr) = gsl_integration_qag (
  94                                 sub { 1/$_[0] } ,
  95                                 1, 10, 0, 1e-7, 1000,
  96                                 5, $self->{wspace},
  97                            );
  98     ok_status($status);
  99     my $res = abs($result - log 10);
 100     ok_similar(
 101         [$result],[log 10],
 102         sprintf('gsl_integration_qags: res=%.18f, abserr=%.18f',$res,$abserr),
 103         $abserr
 104     );
 105 }
 106 # f1(x) = x^alpha * log(1/x) */
 107 # integ(f1,x,0,1) = 1/(alpha + 1)^2 */
 108
 109 sub TEST_QAGS3 : Tests {
 110     my $self = shift;
 111     my ($status, $result, $abserr) = gsl_integration_qags (
 112                                 sub { my $x=shift; $x ** 2 * log (1/$x) } ,
 113                                 0, 1, 0, 1e-7, 1000,
 114                                 $self->{wspace},
 115                            );
 116     ok_status($status);
 117     my $res = abs($result - 1/9);
 118     ok_similar(
 119         [$result],[1/9],
 120         sprintf('gsl_integration_qags: res=%.18f, abserr=%.18f',$res,$abserr),
 121         $abserr
 122     );
 123 }
 124
 125
 126 sub TEST_WORKSPACE_ALLOC : Tests {
 127     my $self = shift;
 128     isa_ok($self->{wspace}, 'Math::GSL::Integration');
 129 }
 130
 131 sub TEST_QNG : Tests {
 132     my ($status, $result, $abserr, $neval) = gsl_integration_qng (
 133                                 sub { $_[0]**2} ,
 134                                 0, 1, 0, 1e-7,
 135                            );
 136     ok_status($status);
 137     my $res = abs($result - 1/3);
 138     ok_similar(
 139         [$result], [1/3],
 140         sprintf('gsl_integration_qng: res=%.18f, abserr=%.18f',$res,$abserr),
 141        $abserr
 142     );
 143 }
 144
 145 sub QAWS_ALLOC : Tests {
 146    my $table = gsl_integration_qaws_table_alloc(0, 0, 1, 0);
 147    isa_ok($table, 'Math::GSL::Integration');
 148 }
 149
 150 sub QAWO_ALLOC : Tests {
 151    my $table = gsl_integration_qawo_table_alloc(10.0 * $M_PI, 1.0,$GSL_INTEG_SINE, 1000);
 152    isa_ok($table, 'Math::GSL::Integration');
 153 }
 154 Test::Class->runtests;