Permalink
Browse files

Fiddling with Roots and adding erfc() check from Keith Lofstrom

  • Loading branch information...
1 parent f4061eb commit ee7c8bcbf7b72583a1802a1260fee20bf222bcc6 Jonathan Leto committed Sep 29, 2008
Showing with 56 additions and 2 deletions.
  1. +1 −0 Roots.i
  2. +35 −0 examples/sf/erfc_check
  3. +4 −0 gsl_typemaps.i
  4. +16 −2 t/Roots.t
View
@@ -1,5 +1,6 @@
%module "Math::GSL::Roots"
%include "gsl_typemaps.i"
+%include "typemaps.i"
%{
#include "gsl/gsl_types.h"
#include "gsl/gsl_roots.h"
View
@@ -0,0 +1,35 @@
+#!/usr/bin/perl
+# interfc
+# v0.2 Keith Lofstrom 2008 Sept 18
+# with tweaks by Jonathan Leto
+#
+# Compare the results of numerical integration of
+# (2.0/sqrt(PI)*integral[x,inf](e^(-x^2) ) with the
+# erfc() function.
+#
+# The results are pretty good, within 1e-14 or better
+# for the examples tried
+
+use Math::GSL::SF qw/:all/;
+use Math::GSL::Integration qw/:all/;
+use Math::GSL::Const qw/:all/;
+use strict;
+
+my $workspace = gsl_integration_workspace_alloc(101);
+my $maxdiff = 0;
+
+printf " arg erfc int() difference\n";
+
+for( my $i=-5.2 ; $i<=5.201 ; $i += 0.1 ) {
+ my $e = gsl_sf_erfc($i);
+ my ($status, $result, $abserr ) = gsl_integration_qagiu(
+ sub { exp( -($_[0]**2) ) },
+ $i, 0, 1e-9, 100, $workspace );
+
+ my $yy = $M_2_SQRTPI * $result ;
+ my $diff=$yy-$e ;
+ $maxdiff = abs($diff) > $maxdiff ? abs($diff) : $maxdiff;
+
+ printf "%4.1f%18.12f%18.12f %.18g\n", $i, $e, $yy, $diff ;
+}
+printf "Maximum local error: %.12g\n", $maxdiff;
View
@@ -78,3 +78,7 @@
$1 = &F;
};
+%typemap(in) gsl_function_fdf * {
+ fprintf(stderr, 'FDF_FUNC');
+
+}
View
@@ -5,8 +5,11 @@ use Test::More;
use Math::GSL qw/:all/;
use Math::GSL::Roots qw/:all/;
use Math::GSL::Test qw/:all/;
+use Math::GSL::Errno qw/:all/;
use Data::Dumper;
+BEGIN { gsl_set_error_handler_off(); }
+
sub make_fixture : Test(setup) {
my $self = shift;
$self->{solver} = gsl_root_fsolver_alloc($gsl_root_fsolver_bisection);
@@ -17,6 +20,15 @@ sub teardown : Test(teardown) {
my $self = shift;
}
+sub GSL_FDFSOLVER_BASIC : Tests {
+ my $self = shift;
+ local $TODO = 'this blows up';
+ #ok_status(gsl_root_fdfsolver_set($self->{fdfsolver},
+ # sub { my $x=shift; ($x-3.2)**3 },
+ # 5
+ #));
+
+}
sub GSL_ROOTS_ALLOC_FREE : Tests {
my $self = shift;
my $x = $self->{solver};
@@ -37,14 +49,16 @@ sub GSL_ROOTS_SET : Tests {
sub GSL_ROOT_ITERATE : Tests {
my $self = shift;
- local $TODO = q{???};
my $solver = gsl_root_fsolver_alloc($gsl_root_fsolver_brent);
ok_status(gsl_root_fsolver_set($solver,
sub { my $x=shift; ($x-3.2)**3 },
0, 5
));
# This currently blows up
- #ok_status( gsl_root_fsolver_iterate($solver) );
+ #local $TODO = q{???};
+ #ok_status( gsl_root_fsolver_iterate($solver));
+ my $root = gsl_root_fsolver_root($solver);
+ ok_similar([$root], [2.5], 'gsl_root_fsolver_root');
}
sub SOlVER_TYPES : Tests {

0 comments on commit ee7c8bc

Please sign in to comment.