Math-Function-Roots

 view release on metacpan or  search on metacpan

lib/Math/Function/Roots.pm  view on Meta::CPAN

=head2 I<guess or min/max> Parameters

Most algorithms require an initial range or guesses. If referred to as
'guesses' then the solution (root) need not be in the range
[guess1,guess2]. If a range or min and max are required, then to
solution B<must> lie within [min,max].

=head2 I<epsilon> and I<max_iter> Parameters

Epsilon (I<e>) is used to set the desired accuracy. Less accurate
answers take fewer iterations and are therefore quicker to compute. In
general I<e> referres to the maximum distance from the given solution
to the actual solution. If you need 6 decimals of accuracy, then I<e>
= .000_000_1 is appropriate, this is the default. Calculating a few
decimals beyond what you need is generally recommended to prevent
compounding rounding errors. I<epsilon> is a named parameter to set
I<e> for that particular run of the algorithm, it should always follow
mandatory parameters:

    my $result = bisection( sub{...}, -10, 10, epsilon => .01 );

lib/Math/Function/Roots.pm  view on Meta::CPAN

sub epsilon(;$){
    if( @_ ){
	$E = shift;
	$E = 0 if $E < 0;
    }
    return $E;
}

=pod

Similar to epsilon, the maximum number of iterations an algorithm
should run for may be set with the I<max_iter> named parameter, or
globally with I<max_iter>(i). This maximum is normally used to catch
errors, i.e. when a given function doesn't converge, or when there is
a bug (nah...). Do not use this to control run-time, if you need an
answer faster, use a larger epsilon. The only reason to change this
would be if you had a slowly converging function, and you were willing
to wait for a good answer, then you could raise the maximum to allow
the algorithm to continue working. Default is 5,000.

=cut

lib/Math/Function/Roots.pm  view on Meta::CPAN

	$Max_Iter = shift;
	$Max_Iter = 1 if ($Max_Iter < 1);
    }	
    return $Max_Iter;
}

=head1 PERFORMANCE CHECKING

=head2 last_iter( )

This will return the number of iterations used to find the last
result. This might help to give an indication on how an algorithm
performs on your data. 

=cut

sub last_iter{
    return $Last_Iter;
}

=head1 ALGORITHMS

lib/Math/Function/Roots.pm  view on Meta::CPAN

	    $by = $py;
	}
	else{
	    # If $py != 0, $ay and $by have opposite signs, $py and $ay have same sign
	    # then $py and $by must have opposite signs
	    $a = $p;
	    $ay = $py;
	}
    }
    
    carp "Maximum iterations: possible bad solution";
    return $p;
}

=head2 fixed_point( I<fixed point function, guess> )

The Fixed-Point Iteration algorithm is a fast robust method which,
unfortunately, works on a limited domain of problems, and requires
some algebra. The benefits are that it can converge rapidly, and the
range the root is in does not need to be known, any guess will
converge, eventually.

lib/Math/Function/Roots.pm  view on Meta::CPAN

    for (1..$Max_Iter){
	$Last_Iter = $_;

	# Each iter we compute p = g(p')
	$p = &$g($last_p);
	if( abs( $p - $last_p ) <= $E ){
	    return $p;
	}
	$last_p = $p;
    }
    carp "Maximum iterations: divergence likely";
    return undef;
}

=head2 secant( I<function>, guess1, guess2 >)

The secant method is a simplification of the Newton method, which uses
the derivitive of the function to better predict the root of the
function. The secant method uses a secant (line between two points on
the function) as a substitute for knowing or calculating the
derivative of the function.

lib/Math/Function/Roots.pm  view on Meta::CPAN

	$p0 = $p1;
	$q0 = $q1;
	$q1 = &$f( $p );
	
	if( $q1 eq 0 || abs( $p - $p1 ) <= $E ){
	    return $p;
	}
	
	$p1 = $p;	
    }
    carp "Maximum iterations: divergence likely";
    return undef;
}

=head2 false_position( I<function, min, max> )

False Position is an algorithm similar to Secant, it uses secants
of the function to pick better guesses. The difference is that this
method incorporates the bracketing of the Bisection method, with the
speed of the Secant method.

lib/Math/Function/Roots.pm  view on Meta::CPAN

	    $a = $p;
	    $ay = $py;
	}
	else{
	    # root is in range [a..p]
	    $b = $p;
	    $by = $py;
	}
	$last_py = $py;
    }
    carp "Maximum iterations: possible bad solution";
    return $p;
}

=head2 find()

This a hybrid function which uses a combination of algorithms to find
the root of the given function. Both I<guess1> and I<guess2> are
optional. If one is provided, it is used as an approximate starting
point. If both are given, then they are taken as a range, the root
B<must> be within this range.

It will most likely return the root nearest your guess, but no
guarantees. Don't provide a range with more than one root in it, you
might find one, you might not. More information will give higher
performance and more control over which root is being found, but if
you don't know anything about the function, give it a try without a
guess. Settings from epsilon and maximum iterations apply as normal.

=cut

sub find(&;$$%){
    my $f = shift;
    my ($a,$b);

    # This is totally wrong, need to not assign to $a and $b when no 
    # arguments

t/01.bisection.t  view on Meta::CPAN

is( last_iter(), 5 );

# shouldn't execute with bad a/b
eval { bisection( sub {shift()}, 1, 2 ) };
like( $@, qr/^Bad range/, "bisection: bad range" );

# test finding perfect exit value
is( bisection( sub{ shift() - 1}, 0, 4 ), 1, "bisection: perfect exit");
is( last_iter(), 2, "perfect exit takes 2 iterations");


# Test standard use with normal epsilon
max_iter(50_000);
my $e_test = .0001;
cmp_ok( abs( 
	 bisection( sub{ 2*shift()+2 }, -10, 10, epsilon => $e_test ) 
	 + 1 ), '<=', $e_test, "normal bisection" );
# These two test should have the same answer, testing the two ways
# of setting epsilon, global and local