t/Integration.t

   1 package Math::GSL::Integration::Test;
   2 use base 'Test::Class';
   3 use Test::More tests => 29;
   4 use Test::Exception;
   5 use Math::GSL              qw/:all/;
   6 use Math::GSL::Test        qw/:all/;
   7 use Math::GSL::Errno       qw/:all/;
   8 use Math::GSL::Const       qw/:all/;
   9 use Math::GSL::Integration 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 sub TEST_CONSTANTS : Tests {
  26     ok(defined $GSL_INTEG_COSINE  , '$GSL_INTEG_COSINE');
  27     ok(defined $GSL_INTEG_SINE    , '$GSL_INTEG_SINE');
  28     ok(defined $GSL_INTEG_GAUSS15 , '$GSL_INTEG_GAUSS15');
  29     ok(defined $GSL_INTEG_GAUSS21 , '$GSL_INTEG_GAUSS21');
  30     ok(defined $GSL_INTEG_GAUSS31 , '$GSL_INTEG_GAUSS31');
  31     ok(defined $GSL_INTEG_GAUSS41 , '$GSL_INTEG_GAUSS41');
  32     ok(defined $GSL_INTEG_GAUSS51 , '$GSL_INTEG_GAUSS51');
  33     ok(defined $GSL_INTEG_GAUSS61 , '$GSL_INTEG_GAUSS61');
  34 }
  35
  36 sub TEST_QAG : Tests {
  37     my $self = shift;
  38     my ($status, $result, $abserr) = gsl_integration_qag (
  39                                 sub { $_[0]**2} ,
  40                                 0, 1, 0, 1e-7, 1000,
  41                                 6, $self->{wspace},
  42                            );
  43     ok_status($status);
  44     my $res = abs($result - 1/3);
  45     ok_similar(
  46         [$result], [1/3],
  47         sprintf('gsl_integration_qag: res=%.18f, abserr=%.18f',$res,$abserr),
  48         $abserr
  49     );
  50 }
  51
  52 sub TEST_QAGI : Tests {
  53     my $integrator = 'gsl_integration_qagi';
  54     my $func = sub { my $x=shift; exp( -$x **2) };
  55     verify_integral($integrator, $func, sqrt($M_PI), 0, 1e-7);
  56     $func = sub { my $x = shift; 1/($x**2 + 1) };
  57     verify_integral($integrator, $func, $M_PI, 0, 1e-7);
  58 }
  59
  60 sub TEST_QAGIU : Tests {
  61     my $self = shift;
  62     my ($status, $result, $abserr) = gsl_integration_qagiu (
  63                                 sub { my $x=shift; log($x)/(1+100+$x**2) } ,
  64                                     0.0, 0.0, 1.0E-3, 1000,
  65                                 $self->{wspace}
  66                             );
  67     ok_status($status);
  68     ok_similar([$abserr],[ 3.016716913328831851E-06], "gsl_integration_qagiu absolute error",1e-5);
  69
  70     my $integrator = 'gsl_integration_qagiu';
  71     my $func = sub { my $x=shift; log($x)/(1+100+$x**2) };
  72     verify_integral($integrator, $func, 3.616892186127022568E-01, 0, 2e-3, 0);
  73
  74     ok_similar([$result],[ 3.616892186127022568E-01], "gsl_integration_qagiu",2e-3);
  75
  76 }
  77
  78 sub verify_integral {
  79     my ($integrator,$func,$actual,$epsabs,$epsrel,$lower,$upper,$key) = @_;
  80     my $wspace = gsl_integration_workspace_alloc(1000);
  81     my $format = "$integrator: actual=%.9f relerr=%.18f abserr=%.18f";
  82     my ($status, $result, $abserr, @params);
  83     if ( $integrator eq 'gsl_integration_qag' ) {
  84          push @params, $func, $lower, $upper, $epsabs, $epsrel, 1000, $key, $wspace;
  85     } elsif ($integrator eq 'gsl_integration_qags' ) {
  86          push @params, $func, $lower, $upper, $epsabs, $epsrel, 1000, $wspace;
  87     } elsif ($integrator eq 'gsl_integration_qagi' ) {
  88          push @params, $func, $epsabs, $epsrel, 1000, $wspace;
  89     } elsif ($integrator eq 'gsl_integration_qagiu' ) {
  90          push @params, $func, $lower, $epsabs, $epsrel, 1000, $wspace;
  91     }
  92     { no strict 'refs'; ($status, $result, $abserr) = $integrator->( @params ) }
  93     ok_status($status);
  94     ok_similar( [$result], [$actual],
  95         sprintf( $format,$actual,
  96             abs($result-$actual)/abs($actual),$abserr
  97         ), $epsrel
  98     );
  99 }
 100
 101
 102 sub TEST_QAG2 : Tests {
 103     my $self = shift;
 104     my ($status, $result, $abserr) = gsl_integration_qag (
 105                                 sub { 1/$_[0] } ,
 106                                 1, 10, 0, 1e-7, 1000,
 107                                 5, $self->{wspace},
 108                            );
 109     ok_status($status);
 110     my $res = abs($result - log 10);
 111     ok_similar(
 112         [$result],[log 10],
 113         sprintf('gsl_integration_qags: res=%.18f, abserr=%.18f',$res,$abserr),
 114         $abserr
 115     );
 116 }
 117 # f1(x) = x^alpha * log(1/x) */
 118 # integ(f1,x,0,1) = 1/(alpha + 1)^2 */
 119
 120 sub TEST_QAGS3 : Tests {
 121     my $self = shift;
 122     my ($status, $result, $abserr) = gsl_integration_qags (
 123                                 sub { my $x=shift; $x ** 2 * log (1/$x) } ,
 124                                 0, 1, 0, 1e-7, 1000,
 125                                 $self->{wspace},
 126                            );
 127     ok_status($status);
 128     my $res = abs($result - 1/9);
 129     ok_similar(
 130         [$result],[1/9],
 131         sprintf('gsl_integration_qags: res=%.18f, abserr=%.18f',$res,$abserr),
 132         $abserr
 133     );
 134 }
 135
 136
 137 sub TEST_WORKSPACE_ALLOC : Tests {
 138     my $self = shift;
 139     isa_ok($self->{wspace}, 'Math::GSL::Integration');
 140 }
 141
 142 sub TEST_QNG : Tests {
 143     my ($status, $result, $abserr, $neval) = gsl_integration_qng (
 144                                 sub { $_[0]**2} ,
 145                                 0, 1, 0, 1e-7,
 146                            );
 147     ok_status($status);
 148     my $res = abs($result - 1/3);
 149     ok( $neval > 0, 'returned number of evaluations');
 150     ok_similar(
 151         [$result], [1/3],
 152         sprintf('gsl_integration_qng: res=%.18f, abserr=%.18f',$res,$abserr),
 153        $abserr
 154     );
 155 }
 156
 157 sub QAWS_ALLOC : Tests {
 158    my $table = gsl_integration_qaws_table_alloc(0, 0, 1, 0);
 159    isa_ok($table, 'Math::GSL::Integration');
 160 }
 161
 162 sub QAWO_ALLOC : Tests {
 163    my $table = gsl_integration_qawo_table_alloc(10.0 * $M_PI, 1.0,$GSL_INTEG_SINE, 1000);
 164    isa_ok($table, 'Math::GSL::Integration');
 165 }
 166 Test::Class->runtests;