Perl Monk, Perl Meditation PerlMonks

### Numerical integration

by neniro (Priest)
 on Dec 24, 2005 at 09:43 UTC Need Help??
 Description: That code was my answer to a question on how to do numerical integration using perl:
```#!/usr/bin/perl

use strict;
use warnings;

package Range;

sub new {
my \$class       = shift;
my \$lower       = shift;
my \$upper       = shift;
my \$delta       = shift || 0.0001;

my \$self = {
lower => \$lower,
upper => \$upper,
delta => \$delta,
current => \$lower,
};
return bless \$self, \$class;
}

sub next {
my \$self = shift;
if (\$self->{current} + \$self->{delta} < \$self->{upper}) {
\$self->{current} += \$self->{delta};
return 1;
} else {
return undef;
}
}

package main;

sub integrate {
my \$function    = shift;
my \$iter        = shift;
my \$sum;

do {
\$sum += \$function->(\$iter->{current}) * \$iter->{delta};
} while (\$iter->next());
return \$sum;
}

print sprintf("%.2f", integrate(sub {(shift)**2}, new Range (2, 5))),
+"\n";
```
Replies are listed 'Best First'.
Re: Numerical integration
by educated_foo (Vicar) on Dec 27, 2005 at 13:44 UTC
Okay, I wasted way too much time playing with this. First, it always bugs me when people reinvent perl's bulitins, and perl already has ranges:
```sub integrate(&@)
{
local \$x;
my \$sum;
my \$f = shift;
my %o = (from => 0, by => 0.01, @_);
for (\$o{from} / \$o{by} .. \$o{to} / \$o{by}) {
\$x = \$_ * \$o{by};
\$sum += &\$f * \$o{by};
}
\$sum;
}

#> integrate { \$x } from => 0, to => 1
## 0.505
And then I thought "why not make this multidimensional?
```sub integrate(&@)
{
my \$f = shift;
## Set up some sensical defaults for ranges and vars:
my (\$from, \$to, \$by, \$vars) = do {
my %o = @_;
for (qw(from to by vars)) {
\$o{\$_} = [\$o{\$_}] if exists \$o{\$_} && ! ref \$o{\$_};
}
my \$dim = @{\$o{from}};
\$o{vars} = [qw(x y z w)[0..\$dim-1]] unless \$o{vars};
\$o{from} = [(0) x \$dim] unless \$o{from};
\$o{by} = [(0.01) x \$dim] unless \$o{by};
@o{qw(from to by vars)};
};
my \$vol = 1;
\$vol *= \$_ for @\$by;
## Generate nested evaluation loops:
local *intgen = sub {
my (\$n, \$body) = @_;
if (\$n < 0) {
'\$sum += &\$f * ' . \$vol;
} else {
my (\$by, \$v) = (\$by->[\$n], \$vars->[\$n]);
my (\$lo,\$hi) = (\$from->[\$n] / \$by, \$to->[\$n] / \$by);
"for \\$\$v (\$lo .. \$hi) { \\$\$v *= \$by ;\n"
. intgen(\$n-1)
. "\n}\n";
}
};

## Do it:
(eval 'sub { my \$f = shift; my \$sum = 0;'.intgen(\$#{\$from}).' \$sum
+ }')
->(\$f);
}

#> integrate { \$x * \$y } from => [0,0], to => [1,1]
## 0.255025
Lightly tested, and does no error checking, but it was fun to build.
Re: Numerical integration
by blazar (Canon) on Dec 27, 2005 at 12:09 UTC

Indeed an interesting and clean approach!! As minor remark for a side note, when doing numerical integration one may likely be concerned about speed - not that relevant after all: I know full well that "premature optimization is the root of all evil", and if one is really concerned then he'd probably not doing this in Perl at all...

But all those dereferentiations and method accesses do not make for *speed* and apart a clean syntax I'm not really sure about what this gives you more than a C-style for loop which for once seems appropriate! (Most often in my experience it is not.)

Whatever, I stress once again that especially for educational purposes I find it great. In this sense, just for fun here's my rewrite of your code according to my own stylistic preferences. It should be just the same as yours modulo one tiny change in logic ("adaptive delta")

.
```#!/usr/bin/perl

use strict;
use warnings;

package Range;

sub new {
my (\$class, \$lower, \$upper, \$delta) = @_;
bless {
lower => \$lower,
upper => \$upper,
delta => \$delta || (\$upper-\$lower)/1000,
current => \$lower,
}, \$class;
}

sub next {
my \$self = shift;
(\$self->{current} += \$self->{delta}) < \$self->{upper};
}

package main;

sub integrate {
my (\$function, \$iter) = @_;
my \$sum;
\$sum += \$function->(\$iter->{current}) * \$iter->{delta} while \$iter
+->next;
\$sum;
}

printf "%.2f\n", integrate sub {(shift)**2}, new Range (2, 5);

__END__

Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: snippet [id://518887]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others romping around the Monastery: (4)
As of 2022-12-10 08:58 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?

No recent polls found

Notices?