diff --git a/macros/math/PGnumericalmacros.pl b/macros/math/PGnumericalmacros.pl index 87e6791fd..ecab029ad 100644 --- a/macros/math/PGnumericalmacros.pl +++ b/macros/math/PGnumericalmacros.pl @@ -1,7 +1,11 @@ =head1 NAME -Numerical methods for the PG language +PGnumericalmacros.pl - Numerical methods for the PG language + +=head1 DESCRIPTION + +The functions in this macro are perl subroutines of standard numerical methods. =cut @@ -9,18 +13,17 @@ =head1 NAME sub _PGnumericalmacros_init { } -=head2 Interpolation methods +=head1 INTERPOLATION METHODS -=head3 plot_list +=head2 plot_list Usage: - plot_list([x0,y0,x1,y1,...]); - plot_list([(x0,y0),(x1,y1),...]); - plot_list(\x_y_array); + plot_list([x0, y0, x1, y1, ...]); + plot_list([(x0, y0), (x1, y1), ...]); - plot_list([x0,x1,x2...], [y0,y1,y2,...]); - plot_list(\@xarray,\@yarray); + plot_list([x0, x1, x2, ...], [y0, y1, y2, ...]); + plot_list(~~@xarray, ~~@yarray); It is important that the x values in any form are unique or this method fails. There is no check for this however. @@ -68,7 +71,11 @@ sub plot_list { }; } -=head3 horner +=head2 horner + +Generates a subroutine which evaluates a polynomial passing through the points +C<(x0,q0), (x1,q1), (x2, q2)>, ... using Horner's method, which is a form of the polynomial +that is used for evaluation with fewer operations than the standard form. Usage: @@ -78,9 +85,6 @@ =head3 horner &$fn(x) = q0 + q1*(x-x0) +q2*(x-x1)*(x-x0) + ...; -Generates a subroutine which evaluates a polynomial passing through the points -C<(x0,q0), (x1,q1), (x2, q2)>, ... using Horner's method. - The array refs for C and C can be any length but must be the same length. Example @@ -108,15 +112,15 @@ sub horner { }; } -=head3 hermite +=head2 hermite + +Produces a reference to polynomial function with the specified values and first derivatives +at (x0,x1,...). C<&$poly(34)> gives a number. Usage: $poly = hermite([x0,x1...],[y0,y1...],[yp0,yp1,...]); -Produces a reference to polynomial function with the specified values and first derivatives -at (x0,x1,...). C<&$poly(34)> gives a number - Generates a subroutine which evaluates a polynomial passing through the specified points with the specified derivatives: (x0,y0,yp0) ... The polynomial will be of high degree and may wobble unexpectedly. Use the Hermite splines @@ -162,15 +166,15 @@ sub hermite { return horner(\@zvals, \@output); } -=head3 hermite_spline - -Usage - - $spline = hermit_spline([x0,x1...],[y0,y1...],[yp0,yp1,...]); +=head2 hermite_spline Produces a reference to a piecewise cubic hermit spline with the specified values and first derivatives at (x0,x1,...). +Usage: + + $spline = hermit_spline([x0,x1...],[y0,y1...],[yp0,yp1,...]); + C<&$spline(45)> evaluates to a number. Generates a subroutine which evaluates a piecewise cubic polynomial @@ -225,7 +229,10 @@ sub hermite_spline { }; } -=head3 cubic_spline +=head2 cubic_spline + +Generic the standard cubic spline of a pair of array references for the x and y +values. Usage: @@ -379,17 +386,195 @@ sub javaScript_cubic_spline { return $output_str; } -=head2 Numerical Integration methods +=head2 newtonDividedDifference + +Computes the Newton divided difference table. + +B + +=over + +=item * C an array reference for x-values. + +=item * C an array reference for y-values. This is the first row/column in the divided +difference table. + +=back + +B + +An array reference of array referencess of divided differences. + +B + + $x = [0, 1, 3, 6]; + $y = [0 ,1 ,2, 5]; + + $c = newtonDividedDifference($x, $y); + +The result of C<$c> is + + [ + [0, 1, 2, 5], + [1, 0.5, 1], + [-0.1667, 0.1], + [0.0444] + ] + +This is generally laid out in the following way: + + 0 0 + 1 + 1 1 -0.1667 + 0.5 0.04444 + 3 2 0.1 + 1 + 6 5 + +where the first column is C<$x>, the second column is C<$y> and the rest of the table is -=head3 lefthandsum + f[x_i, x_j] = (f[x_j] - f[x_i])/(x_j - x_i) -Left Hand Riemann Sum +This table can then be used to find the interpolated polynomial through the given points. For example in the +above example, the polynomial through C<$x> and C<$y> is + + 0 + 1x -0.1167x(x-1) + 0.04444x(x-1)(x-2). + +=cut + +sub newtonDividedDifference { + my ($x, $y) = @_; + my $a = [ [@$y] ]; + for my $j (0 .. (scalar(@$x) - 2)) { + for my $i (0 .. (scalar(@$x) - ($j + 2))) { + $a->[ $j + 1 ][$i] = ($a->[$j][ $i + 1 ] - $a->[$j][$i]) / ($x->[ $i + $j + 1 ] - $x->[$i]); + } + } + return $a; +} + +=head2 legendreP + +Returns a code reference to the Legendre Polynomial of degree C. Usage: - lefthandsum(function_reference, start, end, steps=>30 ); + $poly = legendreP($n); + +For example, if C<$poly = legendreP(7);>, then C<&$poly(0.5)> returns the value of the 7th degree +Legendre polynomial at x = 0.5. Legendre polynomials have a standard domain of [-1,1]. +This subroutine does not check for that. + +=cut + +# This uses the recurrence formula (n+1)P_{n+1}(x) = (2n+1)P_n(x) - n P_{n-1}(x), with P_0(x) = 1 and P_1(x) = x. +# After testing, this is found to have less round-off error than other formula. + +sub legendreP { + my ($n) = @_; + return sub { + my ($x) = @_; + return 1 if $n == 0; + return $x if $n == 1; + my $P1 = legendreP($n - 1); + my $P2 = legendreP($n - 2); + return ((2 * $n - 1) * $x * &$P1($x) - ($n - 1) * &$P2($x)) / $n; + }; +} + +=head2 diffLegendreP + +Returns a code reference to the derivative of the Legendre polynomial of degree C. + +Usage: + + $dp = diffLegendreP($n); + +If C<$dp = diffLegendreP(5)>, then C<&$dp(0.5)> will find the value of the derivative of the 5th +degree Legendre polynomial at C. + +=cut + +# This uses the recurrence relation P'_{n+1}(x) = (n+1)P_n(x) + x P'_n(x). Like the subroutine +# legendreP, it was found that round off error is smaller for this method than others. +sub diffLegendreP { + my ($n) = @_; + return sub { + my ($x) = @_; + return 0 if $n == 0; + my $P = legendreP($n - 1); + my $dP = diffLegendreP($n - 1); + return $n * &$P($x) + $x * &$dP($x); + }; +} + +=head2 legendreP_nodes_weights + +Finds the nodes (roots) and weights of the Legendre Polynomials of degree C. These are used in +Gaussian Quadrature. + +Usage: -Implements the Left Hand sum using 30 intervals between 'start' and 'end'. + ($nodes, $weights) = legendreP_nodes_weights($n); + +The C<$nodes> and C<$weights> are array references of nodes and weights. + +=cut + +# This calculates the roots and weights of the Legendre polynomial of degree n. The roots +# can be determined exactly for n<=9, due to symmetry, however, this uses newton's method +# to solve them based on an approximate value +# (see https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial ) +# +# The weights can then be calculated based on a formula shown in +# https://en.wikipedia.org/wiki/Gaussian_quadrature . + +sub legendreP_nodes_weights { + my ($n) = @_; + + my $leg = legendreP($n); + my $dleg = diffLegendreP($n); + my $pi = 4 * atan(1.0); + + my @nodes; + my @weights; + my $m; + # If $n is odd, then there is a node at x=0. + if ($n % 2 == 1) { + push(@nodes, 0); + push(@weights, 2 / &$dleg(0)**2); + $m = ($n + 1) / 2 + 1; + } else { + $m = $n / 2 + 1; + } + # Compute only nodes for half of the nodes and use symmetry to fill in the rest. + for my $k ($m .. $n) { + my $node = newton( + $leg, $dleg, + (1 - 1 / (8 * $n**2) + 1 / (8 * $n**3)) * cos($pi * (4 * $k - 1) / (4 * $n + 2)), + feps => 1e-14 + )->{root}; + my $w = 2 / ((1 - $node**2) * &$dleg($node)**2); + unshift(@nodes, $node); + push(@nodes, -$node); + + unshift(@weights, $w); + push(@weights, $w); + } + return (\@nodes, \@weights); +} + +=head1 NUMERICAL INTEGRATION METHODS + +=head2 lefthandsum + +Calculate the Left Hand Riemann Sum + +Usage: + + lefthandsum(function_reference, start, end, steps => 30); + +Implements the Left Hand sum using the given number of intervals between 'start' and 'end'. The first three arguments are required. The final argument (number of steps) is optional and defaults to 30. @@ -409,15 +594,15 @@ sub lefthandsum { return $sum * $delta; } -=head3 righthandsum +=head2 righthandsum -Right Hand Riemann Sum +Calculate the Right Hand Riemann Sum Usage: - righthandsum(function_reference, start, end, steps=>30 ); + righthandsum(function_reference, start, end, steps => 30); -Implements the right hand sum using 30 intervals between 'start' and 'end'. +Implements the right hand sum using the given number of intervals between 'start' and 'end'. The first three arguments are required. The final argument (number of steps) is optional and defaults to 30. @@ -437,15 +622,16 @@ sub righthandsum { return $sum * $delta; } -=head3 midpoint +=head2 midpoint + +Approximate the integral of a function using the midpoint rule. Usage: - midpoint(function_reference, start, end, steps=>30); + midpoint(function_reference, start, end, steps => 30); -Implements the Midpoint rule between 'start' and 'end'. -The first three arguments are required. The final argument (number of steps) -is optional and defaults to 30. +Implements the midpoint rule between 'start' and 'end' using the given number of steps. +The first three arguments are required. The number of steps is optional and defaults to 30. =cut @@ -463,13 +649,15 @@ sub midpoint { return $sum * $delta; } -=head3 simpson +=head2 simpson + +Approximate the integral of a function using Simpson's rule. Usage: - simpson(function_reference, start, end, steps=>30 ); + simpson(function_reference, start, end, steps => 30); -Implements Simpson's rule between 'start' and 'end'. +Implements Simpson's rule between 'start' and 'end' using the given number of steps. The first three arguments are required. The final argument (number of steps) is optional and defaults to 30, but must be even. @@ -494,13 +682,15 @@ sub simpson { return $sum * $delta / 3; } -=head3 trapezoid +=head2 trapezoid + +Approximate the integral of a function using the Trapezoid rule. Usage: - trapezoid(function_reference, start, end, steps=>30); + trapezoid(function_reference, start, end, steps => 30); -Implements the trapezoid rule using 30 intervals between 'start' and 'end'. +Implements the trapezoid rule between 'start' and 'end' using the given number of steps. The first three arguments are required. The final argument (number of steps) is optional and defaults to 30. @@ -521,13 +711,15 @@ sub trapezoid { return $sum * $delta; } -=head3 romberg +=head2 romberg + +Implements the Romberg integration routine through 'level' recursive steps. Level defaults to 6. Usage: - romberg(function_reference, x0, x1, level); + romberg(function_reference, x0, x1, level => n); + -Implements the Romberg integration routine through 'level' recursive steps. Level defaults to 6. =cut @@ -545,24 +737,24 @@ sub romberg_iter { / (4**($k - 1) - 1); } -=head3 inv_romberg +=head2 inv_romberg -Inverse Romberg +Finds the upper limit of an integral using Romberg integration and newton's method. Usage: inv_romberg(function_reference, a, value); -Finds b such that the integral of the function from a to b is equal to value. -Assumes that the function is continuous and doesn't take on the zero value. -Uses Newton's method of approximating roots of equations, and Romberg to evaluate definite integrals. +Finds b such that the integral of the function from a to b is equal to value. Assumes that the +function is continuous and doesn't take on the zero value. Uses Newton's method of approximating +roots of equations, and Romberg to evaluate definite integrals. -Example +Example: -Find the value of b such that the integral of e^(-x^2/2)/sqrt(2*pi) from 0 to b is 0.25. +Find the value of b such that the integral of e^(-x^2/2)/sqrt(2*pi) from 0 to b is 0.45. - $f = sub { my $x = shift; return exp(-$x*$x/2)/sqrt(4*acos(0));}; - $b = inv_romberg($f,0,0.45); + $f = sub { my $x = shift; return exp(-$x * $x / 2) / sqrt(4 * acos(0)); }; + $b = inv_romberg($f, 0, 0.45); this returns C<1.64485362695934>. This is the standard normal curve and this value is the z value for the 90th percentile. @@ -590,26 +782,144 @@ sub inv_romberg { return $b; } -=head2 Differential Equation Methods +=head2 newtonCotes + +Perform quadrature (numerical integration) using a newtonCotes composite formula (trapezoid, +Simpson's, the 3/8 rule or Boole's). + +Usage: + + newtonCotes($f, $a, $b, n => 4, method => 'simpson'); + +where C<$f> is a subroutine reference (function that takes a single numerical value and +returns a single value), and C<[$a, $b]> is an interval to integrate over. + +B + +=over + +=item method + +The method options are either open or closed methods. The closed newton-cotes formula methods +are C. The open newton-cotes formula methods are +C, where the number indicates the number of used nodes for the formula. + +=item n + +This is the number of subintervals to use for a composite version of the formula. +If n is set to 1, then this uses the non-composite version of the method. + +=back + +=cut + +sub newtonCotes { + my ($f, $a, $b, @args) = @_; + my %opts = (n => 10, method => 'simpson', @args); + my $h = ($b - $a) / $opts{n}; + my @weights; + my @innernodes; + + if ($opts{method} eq 'trapezoid') { + @weights = (1 / 2, 1 / 2); + @innernodes = (0, 1); + } elsif ($opts{method} eq 'simpson') { + @weights = (1 / 6, 4 / 6, 1 / 6); + @innernodes = (0, 0.5, 1); + } elsif ($opts{method} eq 'three-eighths') { + @weights = (1 / 8, 3 / 8, 3 / 8, 1 / 8); + @innernodes = (0, 1 / 3, 2 / 3, 1); + } elsif ($opts{method} eq 'boole') { + @weights = (7 / 90, 32 / 90, 12 / 90, 32 / 90, 7 / 90); + @innernodes = (0, 1 / 4, 1 / 2, 3 / 4, 1); + } elsif ($opts{method} eq 'open1') { + @weights = (undef, 1); + @innernodes = (undef, 0.5); + } elsif ($opts{method} eq 'open2') { + @weights = (undef, 1 / 2, 1 / 2); + @innernodes = (undef, 1 / 3, 2 / 3); + } elsif ($opts{method} eq 'open3' || $opts{method} eq 'milne') { + @weights = (undef, 2 / 3, -1 / 3, 2 / 3); + @innernodes = (undef, 1 / 4, 1 / 2, 3 / 4); + } elsif ($opts{method} eq 'open4') { + @weights = (undef, 11 / 24, 1 / 24, 1 / 24, 11 / 24); + @innernodes = (undef, 1 / 5, 2 / 5, 3 / 5, 4 / 5); + } + + my $quad = 0; + for my $i (0 .. $opts{n} - 1) { + for my $k (0 .. $#innernodes) { + $quad += &$f($a + ($i + $innernodes[$k]) * $h) * $weights[$k] if $weights[$k]; + } + } + return $h * $quad; +} + +=head2 gaussQuad + +Compute the integral of a function C<$f> on an interval C<[a, b]> using Gassian Quadrature. + +Usage: + + gaussQuad($f, n => 5, a => -1, b => 1, weights => $w, nodes => $nodes); -=head3 rungeKutta4 +where C<$f> is a code reference to a function from R => R, C and C are the endpoints of the +interval, C is the number of nodes to use. The weights and nodes will depend on the value of +C. -Finds integral curve of a vector field using the 4th order Runge Kutta method by -providing the function C +If C or C are included, they must both be used and will override the C option. +These will not be checked and assumed to be correct. These should be used for performance +in that calculating the weights and nodes have some computational time. + +=cut + +sub gaussQuad { + my ($f, %opts) = @_; + # defines default values. + %opts = (n => 5, a => -1, b => 1, %opts); + die 'The optional value n must be an integer >=2' unless $opts{n} =~ /\d+/ && $opts{n} >= 2; + die 'The optional value a must be a number' unless $opts{a} =~ /[+-]?\d*\.?\d+/; + die 'The optional value b must be a number' unless $opts{b} =~ /[+-]?\d*\.?\d+/; + die 'The optional value b must be greater than a' unless $opts{b} > $opts{a}; + die 'The argument f must be a code ref' unless ref($f) eq 'CODE'; + + my ($x, $w) = ($opts{nodes}, $opts{weights}); + if ((!defined($w) && !defined($x))) { + ($x, $w) = legendreP_nodes_weights($opts{n}); + } elsif (!defined($w) || !defined($w)) { + die 'If either option "weights" or "nodes" is used, both must be used.'; + } + die 'The options weights and nodes must be array refs of the same length' + unless ref $w eq 'ARRAY' && ref $x eq 'ARRAY' && scalar($x) == scalar($x); + + my $sum = 0; + $sum += $w->[$_] * &$f(0.5 * ($opts{b} + $opts{a}) + 0.5 * ($opts{b} - $opts{a}) * $x->[$_]) + for (0 .. scalar(@$w) - 1); + return 0.5 * ($opts{b} - $opts{a}) * $sum; +} + +=head1 DIFFERENTIAL EQUATION METHODS + +=head2 rungeKutta4 + +Finds integral curve of a vector field using the 4th order Runge Kutta method. In particular +find an approximate solution to + +dx/dt = F(t,x) Usage: - rungeKutta4( &vectorField(t,x),%options); + rungeKutta4(&F(t, x), %options); -Returns: array ref of points [t,y] +Default %options: + 'initial_t' => 1, + 'initial_y' => 1, + 'dt' => 0.01, + 'num_of_points' => 10, # number of reported points + 'interior_points' => 5, # number of 'interior' steps between reported points + 'debug' - Default %options: - 'initial_t' => 1, - 'initial_y' => 1, - 'dt' => 0.01, - 'num_of_points' => 10, # number of reported points - 'interior_points' => 5, # number of 'interior' steps between reported points - 'debug' +Returns: array ref of points [t,y] =cut @@ -636,7 +946,7 @@ sub rungeKutta4 { my ($out, $err) = &$rf_fun(@in); $errors .= " $err at ( " . join(" , ", @in) . " )
\n" if defined($err); $out = 'NaN' if defined($err) and not is_a_number($out); - $out; + return $out; }; my @output = ([ $t, $y ]); @@ -658,4 +968,376 @@ sub rungeKutta4 { } } +=head2 solveDiffEqn + +Produces a numerical solution to the differential equation y' = f(x,y) using a number of optional +methods. + +B + +=over + +=item * C a subroutine reference that take two inputs (x,y) and returns +a single number. Note: if you use a Formula to generate a function, create a perl +function with the C<<$f->perlFunction>> method. + +=item * C a real-valued number for the initial point + +=back + +B + +=over + +=item * C the initial x-value (defaults to 0) + +=item * C the step size of the numerical method (defaults to 0.25) + +=item * C the number of steps to perform (defaults to 4) + +=item * C one of 'euler', 'improved_euler', 'heun' or 'rk4' (defaults to euler) + +=back + +B + +An hash with the following fields: + +=over + +=item * C an array reference of the x-values which are C + +=item * C an array reference of the y-values (which depend on the method used) + +=item * C the intermediate function values used (which depend on the method). + +=back + +B + +The following performs Euler's method on C using C for C points, +so the last x-value is 5. + + $f = sub { my ($x, $y) = @_; return $x*$y; } + $sol1 = solveDiffEqn($f, 1, x0 => 0, h => 0.5, n => 10, method => 'euler'); + +The output C<$sol> is a hash reference with fields x and y, where each have 11 points. + +The following uses the improved Euler method on C using C for +C points (the last x-value is 1.0). Note, this shows how to pass the perl function to the +method. + + Context()->variables->add(y => 'Real'); + $G = Formula("x^2 + y^2"); + $g = $G->perlFunction; + $sol2 = solveDiffEqn($g, 1, method => 'improved_euler', x0 => 0, h => 0.2,n => 5); + +In this case, C<$sol2> returns both x and y, but also, the values of C and C. + +=cut + +sub solveDiffEqn { + my ($f, $y0, @args) = @_; + my %opts = (x0 => 0, h => 0.25, n => 4, method => 'euler', @args); + + die 'The first argument must be a subroutine reference' unless ref($f) eq 'CODE'; + die 'The option n must be a positive integer' unless $opts{n} =~ /^\d+$/; + die 'The option h must be a positive number' unless $opts{h} > 0; + die 'The option method must be one of euler/improved_euler/heun/rk4' + unless grep { $opts{method} eq $_ } qw/euler improved_euler heun rk4/; + + my $x0 = $opts{x0}; + my $h = $opts{h}; + my @y = ($y0); + my @k1; + my @k2; + my @k3; + my @k4; + my @x = map { $x0 + $_ * $h } (0 .. $opts{n}); + + for my $j (1 .. $opts{n}) { + if ($opts{method} eq 'euler') { + $y[$j] = $y[ $j - 1 ] + $h * &$f($x[ $j - 1 ], $y[ $j - 1 ]); + } elsif ($opts{method} eq 'improved_euler') { + $k1[$j] = &$f($x[ $j - 1 ], $y[ $j - 1 ]); + $k2[$j] = &$f($x[$j], $y[ $j - 1 ] + $h * $k1[$j]); + $y[$j] = $y[ $j - 1 ] + 0.5 * $h * ($k1[$j] + $k2[$j]); + } elsif ($opts{method} eq 'heun') { + $k1[$j] = &$f($x[ $j - 1 ], $y[ $j - 1 ]); + $k2[$j] = &$f($x[ $j - 1 ] + 2 * $h / 3, $y[ $j - 1 ] + 2 * $h / 3 * $k1[$j]); + $y[$j] = $y[ $j - 1 ] + 0.25 * $h * ($k1[$j] + 3 * $k2[$j]); + } elsif ($opts{method} eq 'rk4') { + $k1[$j] = &$f($x[ $j - 1 ], $y[ $j - 1 ]); + $k2[$j] = &$f($x[ $j - 1 ] + 0.5 * $h, $y[ $j - 1 ] + $h * 0.5 * $k1[$j]); + $k3[$j] = &$f($x[ $j - 1 ] + 0.5 * $h, $y[ $j - 1 ] + $h * 0.5 * $k2[$j]); + $k4[$j] = &$f($x[$j], $y[ $j - 1 ] + $h * $k3[$j]); + $y[$j] = $y[ $j - 1 ] + $h / 6 * ($k1[$j] + 2 * $k2[$j] + 2 * $k3[$j] + $k4[$j]); + } + } + if ($opts{method} eq 'euler') { + return { y => \@y, x => \@x }; + } elsif ($opts{method} eq 'improved_euler' || $opts{method} eq 'heun') { + return { k1 => \@k1, k2 => \@k2, y => \@y, x => \@x }; + } elsif ($opts{method} eq 'rk4') { + return { + k1 => \@k1, + k2 => \@k2, + k3 => \@k3, + k4 => \@k4, + y => \@y, + x => \@x + }; + } +} + +=head1 ROOTFINDING + +=head2 bisection + +Performs the bisection method for the function C<$f> and initial interval C<$int> (an array +reference). An example is + + $f = sub { $x = shift; $x**2 - 2; } + $bisect = bisection($f, [1, 2]); + +The result is a hash with fields C (the estimated root) and C (an array reference +of intervals for each step of bisection) or a hash with field C if there is an error with +either the inputs or from the method. + +B + +=over + +=item * C, a reference to a subroutine with a single input number and single output +value. + +=item * C, an array reference for the interval C<[a,b]> where a < b. + +=back + +B + +=over + +=item * C, the maximum error of the root or stopping condition. Default is C<1e-6> + +=item * C, the maximum number of iterations to run the bisection method. Default is C<40>. + +=back + +B + +A hash with the following fields + +=over + +=item * C, the approximate root using bisection. + +=item * C, an array reference of the intervals (each interval also an array reference) + +=item * C, a string specifying the error (either argument argument error or too many steps) + +=back + +=cut + +sub bisection { + my ($f, $int, @args) = @_; + my %opts = (eps => 1e-6, max_iter => 40, @args); + + # Check that the arguments/options are valid. + return { error => 'The function must be a code reference' } unless ref($f) eq 'CODE'; + + return { error => 'The interval must be an array ref of length 2' } + unless ref($int) eq 'ARRAY' && scalar(@$int) == 2; + + return { error => 'The initial interval [a, b] must satisfy a < b' } unless $int->[0] < $int->[1]; + + return { error => 'The function may not have a root on the given interval' } + unless &$f($int->[0]) * &$f($int->[1]) < 0; + + return { error => 'The option eps must be a positive number' } unless $opts{eps} > 0; + + return { error => 'The option max_iter must be a positive integer' } + unless $opts{max_iter} > 0 && int($opts{max_iter}) == $opts{max_iter}; + + # stores the intervals for each step + my $ints = [$int]; + my $i = 0; + do { + my $mid = 0.5 * ($ints->[$i][0] + $ints->[$i][1]); + my $fmid = &$f($mid); + push(@$ints, $fmid * &$f($ints->[$i][0]) < 0 ? [ $ints->[$i][0], $mid ] : [ $mid, $ints->[$i][1] ]); + $i++; + } while ($i < $opts{max_iter} + && ($ints->[$i][1] - $ints->[$i][0]) > $opts{eps}); + + if ($i == $opts{max_iter}) { + return { error => "You have reached the maximum number of iterations: $opts{max_iter} without " + . 'reaching a root.' }; + } + + return { + root => 0.5 * ($ints->[$i][0] + $ints->[$i][1]), + intervals => $ints + }; +} + +=head2 newton + +Performs Newton's method for the function C<$f> and initial point C<$x0>. +An example is + + $f = sub { my $x = shift; return $x**2 - 2; }; + $df = sub { my $x = shift; return 2*$x; }; + $newton = newton($f, $df, 1); + +The result is a hash with fields C (the estimated root) and C (an array reference +of the iterations with the first being C<$x0>. The result hash will contain the field C +if there is an error. + +B + +=over + +=item * C, a reference to a subroutine with a single input number and single output value. + +=item * C, a subroutine reference that is the derivative of f. + +=item * C, a perl number or math object number. + +=back + +B + +=over + +=item * C, the maximum number of iterations to run Newton's method. Default is C<15>. + +=item * C, the cutoff value in the C direction or stopping condition. The default is C<1e-8> + +=item * C, the allowed functional value for the stopping condition. The default value is C<1e-10>. + +=back + +B + +A hash with the following fields + +=over + +=item * C, the approximate root + +=item * C, an array reference of the iterations + +=item * C, a string specifying the error (either argument argument error or too many steps) + +=back + +=cut + +sub newton { + my ($f, $df, $x0, @args) = @_; + my %opts = (eps => 1e-8, feps => 1e-10, max_iter => 15, @args); + + # Check that the arguments/options are valid. + return { error => 'The function must be a code reference' } unless ref($f) eq 'CODE'; + + return { error => 'The option eps must be a positive number' } + unless $opts{eps} > 0; + + return { error => 'The option feps must be a positive number' } + unless $opts{feps} > 0; + + return { error => 'The option max_iter must be a positive integer' } + unless $opts{max_iter} > 0; + + my @iter = ($x0); + my $i = 0; + do { + $iter[ $i + 1 ] = $iter[$i] - &$f($iter[$i]) / &$df($iter[$i]); + $i++; + return { error => "Newton's method did not converge in $opts{max_iter} steps" } + if $i > $opts{max_iter}; + } while abs($iter[$i] - $iter[ $i - 1 ]) > $opts{eps} || &$f($iter[$i]) > $opts{feps}; + + return { root => $iter[$i], iterations => \@iter }; +} + +=head2 secant + +Performs the secant method for finding a root of the function C<$f> with initial points C<$x0> and +C<$x1>. An example is + + $f = sub { my $x = shift; return $x**2 - 2; }; + $secant = secant($f, 1, 2); + +The result is a hash with fields C (the estimated root) and C (an array reference +of the iterations with the first two being C<$x0> and C<$x1>. The result hash will contain +the field C if there is an error. + +B + +=over + +=item * C, a reference to a subroutine with a single input number and single output value + +=item * C, a number + +=item * C, a number + +=back + +B + +=over + +=item * C, the maximum number of iterations to run the Secant method. The default value is C<20>. + +=item * C, the cutoff value in the C direction or stopping condition. The default value is C<1e-8>. + +=item * C, the allowed functional value for the stopping condition. The default value is C<1e-10>. + +=back + +B + +A hash with the following fields + +=over + +=item * C, the approximate root. + +=item * C, an array reference to the iterations. + +=item * C, a string specifying the error (either argument argument error or too many steps) + +=back + +=cut + +sub secant { + my ($f, $x0, $x1, @args) = @_; + my %opts = (eps => 1e-8, feps => 1e-10, max_iter => 20, @args); + + # Check that the arguments/options are valid. + return { error => 'The function must be a code reference' } unless ref($f) eq 'CODE'; + return { error => 'The option eps must be a positive number' } unless $opts{eps} > 0; + return { error => 'The option feps must be a positive number' } unless $opts{feps} > 0; + return { error => 'The option max_iter must be a positive integer' } unless $opts{max_iter} > 0; + + my @iter = ($x0, $x1); + my $i = 1; + do { + my $m = (&$f($iter[$i]) - &$f($iter[ $i - 1 ])) / ($iter[$i] - $iter[ $i - 1 ]); + $iter[ $i + 1 ] = $iter[$i] - &$f($iter[$i]) / $m; + $i++; + return { error => "The secant method did not converge in $opts{max_iter} steps" } + if $i > $opts{max_iter}; + + } while abs($iter[$i] - $iter[ $i - 1 ]) > $opts{eps}; + + return { root => $iter[$i], iterations => \@iter }; +} + 1; diff --git a/t/macros/numerical_methods.t b/t/macros/numerical_methods.t index b5a00a2a3..bef08fa24 100644 --- a/t/macros/numerical_methods.t +++ b/t/macros/numerical_methods.t @@ -83,6 +83,21 @@ subtest 'cubic spline' => sub { is &$s(0.5), 0.6875, 'check s(0.5)'; }; +subtest 'Newton Divided difference' => sub { + my @x = (0, 1, 3, 6); + my @y = (0, 1, 2, 5); + my $a = + [ [ 0, 1, 2, 5 ], [ 1, 1 / 2, 1 ], [ -1 / 6, 1 / 10 ], [ 2 / 45 ] ]; + + is newtonDividedDifference(\@x, \@y), $a, 'Newton Divided difference, test 1'; + + @x = (5, 6, 9, 11); + @y = (12, 13, 14, 16); + $a = + [ [ 12, 13, 14, 16 ], [ 1, 1 / 3, 1 ], [ -1 / 6, 4 / 30 ], [ 1 / 20 ] ]; + is newtonDividedDifference(\@x, \@y), $a, 'Newton Divided difference, test 2'; +}; + subtest 'Riemann Sums' => sub { my $f = sub { my $x = shift; return $x * $x; }; is lefthandsum($f, 0, 2, steps => 4), 1.75, 'left hand sum of x^2 on [0,2]'; @@ -93,6 +108,7 @@ subtest 'Riemann Sums' => sub { subtest 'Quadrature' => sub { my $f = sub { my $x = shift; return $x * $x; }; my $g = sub { my $x = shift; return exp($x); }; + my $h = sub { my $x = shift; return exp(-$x * $x / 2) / sqrt(4 * acos(0)); }; is simpson($f, 0, 2, steps => 4), 8 / 3, "Simpson's rule of x^2 on [0,2]"; is Round(simpson($g, 0, 1), 7), Round(exp(1) - 1, 7), "Simpson's rule of e^x on [0,1]"; like dies { simpson($f, 0, 2, steps => 5); }, @@ -105,6 +121,124 @@ subtest 'Quadrature' => sub { is romberg($g, 0, 1), exp(1) - 1, 'Romberg interation on e^x on [0,1]'; is inv_romberg($g, 0, exp(1) - 1), 1.0, 'Inverse Romberg to find b with int of e^x on [0,b] returns 1'; + is inv_romberg($h, 0, 0.45), 1.64485362695934, + 'Inverse Romberg to find b with int of the normal curve on [0,b] returns 1.64485362695934'; + + is newtonCotes($f, 0, 2, n => 4, method => 'trapezoid'), 2.75, 'Newton-Cotes (trapezoid) of x^2 on [0,2]'; + is newtonCotes($f, 0, 2, n => 4, method => 'simpson'), 8 / 3, 'Newton-Cotes (simpson) of x^2 on [0,2]'; + is newtonCotes($f, 0, 2, n => 4, method => 'three-eighths'), 8 / 3, 'Newton-Cotes (3/8) of x^2 on [0,2]'; + is newtonCotes($f, 0, 2, n => 4, method => 'boole'), 8 / 3, 'Newton-Cotes (boole) of x^2 on [0,2]'; + + is newtonCotes($g, -1, 1, n => 1, method => 'trapezoid'), 3.0861612696304874, + 'Newton-Cotes (trapezoid) of e^x on [-1,1]'; + is newtonCotes($g, -1, 1, n => 1, method => 'simpson'), 2.362053756543496, + 'Newton-Cotes (simpsons) of e^x on [-1,1]'; + is newtonCotes($g, -1, 1, n => 1, method => 'three-eighths'), 2.355648119152531, + 'Newton-Cotes (3/8) of e^x on [-1,1]'; + is newtonCotes($g, -1, 1, n => 1, method => 'boole'), 2.350470903569373, + 'Newton-Cotes (boole) of e^x on [-1,1]'; + + is newtonCotes($g, -1, 1, n => 4, method => 'trapezoid'), 2.3991662826140026, + 'Newton-Cotes (composite trapezoid, n=4) of e^x on [-1,1]'; + is newtonCotes($g, -1, 1, n => 4, method => 'simpson'), 2.3504530172422795, + 'Newton-Cotes (composite simpson, n=4) of e^x on [-1,1]'; + is newtonCotes($g, -1, 1, n => 4, method => 'three-eighths'), 2.350424908072871, + 'Newton-Cotes (composite 3/8, n=4) of e^x on [-1,1]'; + is newtonCotes($g, -1, 1, n => 4, method => 'boole'), 2.3504024061087962, + 'Newton-Cotes (composite boole, n=4) of e^x on [-1,1]'; +}; + +subtest 'Quadrature - Open Newton-Cotes' => sub { + my $f = sub { my $x = shift; return $x * $x; }; + my $g = sub { my $x = shift; return exp($x); }; + is newtonCotes($f, 0, 2, n => 1, method => 'open1'), 2, 'Newton-Cotes (open, k=1) of x^2 on [0,2]'; + is newtonCotes($f, 0, 2, n => 1, method => 'open2'), 20 / 9, 'Newton-Cotes (open, k=2) of x^2 on [0,2]'; + is newtonCotes($f, 0, 2, n => 1, method => 'open3'), 8 / 3, 'Newton-Cotes (open, k=3) of x^2 on [0,2]'; + is newtonCotes($f, 0, 2, n => 1, method => 'open4'), 8 / 3, 'Newton-Cotes (open, k=4) of x^2 on [0,2]'; +}; + +subtest 'Legendre Polynomial' => sub { + my $leg3 = legendreP(3); + is &$leg3( 0.5), (5 * (0.5)**3 - 3 * (0.5)) / 2.0, 'testing legendreP(3,0.5)'; + is &$leg3(-0.9), (5 * (-0.9)**3 - 3 * (-0.9)) / 2.0, 'testing legendreP(3,0.5)'; + is &$leg3( 1), 1, 'testing legendreP(3,1)'; + is &$leg3(-1), -1, 'testing legendreP(3,-1)'; + + my $leg6 = legendreP(6); + is &$leg6(0.5), (231 * 0.5**6 - 315 * 0.5**4 + 105 * 0.5**2 - 5) / 16.0, 'testing legendreP(6,0.5)'; + is Round(&$leg6(-0.3), 10), Round((231 * (-0.3)**6 - 315 * (-0.3)**4 + 105 * (-0.3)**2 - 5) / 16.0, 10), + 'testing legendreP(6,-0.3)'; + is &$leg6( 1), 1, 'testing legendreP(6,1)'; + is &$leg6(-1), 1, 'testing legendreP(6,-1)'; + + my $leg12 = legendreP(12); + is Round(&$leg12(0.5), 15), Round(980431 / 4194304, 15), 'evaluating legendreP(12,0.5)'; + is Round(&$leg12(-0.9), 15), Round(41726683414959 / 1024000000000000, 15), 'evaluating legendreP(12,-0.9)'; + + my $dleg3 = diffLegendreP(3); + is &$dleg3( 0.5), (15 * (0.5)**2 - 3) / 2.0, 'testing diffLegendreP(3,0.5)'; + is &$dleg3(-0.9), (15 * (0.9)**2 - 3) / 2.0, 'testing diffLegendreP(3,-0.9)'; + + my $dleg10 = diffLegendreP(10); + is &$dleg10(0.4), -2.70832364, 'testing diffLegendreP(10) at x=0.4'; + + my $dleg12 = diffLegendreP(12); + + is &$dleg12(-0.8), -16152097767 / 3125000000, 'testing diffLegendreP(12) at x=-0.8'; + +}; + +subtest 'Legendre Polynomial Roots and Weights' => sub { + my ($roots5, $weights5) = legendreP_nodes_weights(5); + is $roots5, [ -0.906179845938664, -0.5384693101056831, 0.0, 0.5384693101056831, 0.906179845938664 ], + 'roots of LegendreP(5)'; + is $weights5, + [ 0.23692688505618908, 0.47862867049936647, 0.5688888888888889, 0.47862867049936647, 0.23692688505618908 ], + 'weights of LegendreP(5)'; + my ($roots12, $weights12) = legendreP_nodes_weights(12); + is roundArray($roots12, digits => 14), + roundArray( + [ + -0.9815606342467192, -0.9041172563704748, -0.7699026741943047, -0.5873179542866175, + -0.3678314989981802, -0.1252334085114689, 0.1252334085114689, 0.3678314989981802, + 0.5873179542866175, 0.7699026741943047, 0.9041172563704748, 0.9815606342467192 + ], + digits => 14 + ), + 'roots of LegendreP(12)'; + is roundArray($weights12, digits => 14), + roundArray( + [ + 0.04717533638651175, 0.10693932599531826, 0.16007832854334625, 0.20316742672306587, + 0.23349253653835492, 0.24914704581340288, 0.24914704581340288, 0.23349253653835492, + 0.20316742672306587, 0.16007832854334625, 0.10693932599531826, 0.04717533638651175 + ], + digits => 14 + ), + 'weights of LegendreP(12)'; + +}; + +subtest 'Gaussian Quadrature' => sub { + my $f = sub { my $x = shift; return $x**3; }; + is Round(gaussQuad($f), 15), 0, 'gaussQuad(x^3) on [-1,1]'; + is Round(gaussQuad($f, a => 0, b => 1), 15), 0.25, 'gaussQuad(x^3) on [0,1]'; + + is gaussQuad($f, n => 2, a => 0, b => 1), 0.25, 'gaussQuad(x^3, n=>2) on [0,1]'; + + my $g = sub { my $x = shift; return $x**6; }; + is gaussQuad($g), 2 / 7, 'gaussQuad(x^6) on [-1,1]'; + is Round(gaussQuad($g, n => 2), 15), Round(2 * (1 / sqrt(3))**6, 15), 'gaussQuad(x^6) on [-1,1]'; + + my $e_x = sub { my $x = shift; return exp($x); }; + is Round(gaussQuad($e_x, n => 3), 15), Round(5 * (exp(-sqrt(3 / 5)) + exp(sqrt(3 / 5))) / 9 + 8 / 9, 15), + 'gaussQuad(x^6) on [-1,1]'; + is Round(gaussQuad($e_x, n => 15, a => 0, b => 1), 14), Round(exp(1) - 1, 14), 'gaussQuad(e^x,n=>15) on [-1,1]'; + + my ($nodes, $weights) = legendreP_nodes_weights(14); + is Round(gaussQuad($e_x, a => 0, b => 1, nodes => $nodes, weights => $weights), 14), Round(exp(1) - 1, 14), + 'gaussQuad(e^x,n=>15) on [-1,1]'; + }; subtest 'Runge Kutta 4th order' => sub { @@ -126,6 +260,235 @@ subtest 'Runge Kutta 4th order' => sub { 'returns correct y values'; }; +subtest 'Options for solveDiffEqn' => sub { + my $g = sub { + my ($x, $y) = @_; + return $x**2 + $y**2; + }; + + like dies { + Context()->variables->add(y => 'Real'); + my $f = Formula('x^2+y^2'); + solveDiffEqn($f, 1); + }, qr/The first argument must be a subroutine reference/, 'The first argument must be a sub.'; + like dies { solveDiffEqn($g, 1, n => -3) }, qr/The option n must be a positive integer/, + 'The option n is a positive integer'; + like dies { solveDiffEqn($g, 1, h => -0.25) }, qr/The option h must be a positive number/, + 'The option h is a positive number'; + like dies { solveDiffEqn($g, 1, method => 'error') }, + qr/The option method must be one of euler\/improved_euler\/heun\/rk4/, 'Checking for a value method'; +}; + +subtest "Solve an ODE using Euler's method" => sub { + my $g = sub { + my ($x, $y) = @_; + return $x**2 + $y**2; + }; + + my $soln = solveDiffEqn( + $g, 1, + method => 'euler', + h => 0.2, + n => 5 + ); + is $soln->{x}, [ 0, 0.2, 0.4, 0.6, 0.8, 1.0 ], 'returns correct x'; + is roundArray($soln->{y}), + roundArray([ 1, 1.2, 1.496, 1.9756032, 2.8282048008, 4.5559532799 ]), + 'returns correct y'; +}; + +subtest 'Solve an ODE using improved Euler\'s method ' => sub { + my $g = sub { + my ($x, $y) = @_; + return $x**2 + $y**2; + }; + + my $soln = solveDiffEqn( + $g, 1, + x0 => 0, + method => 'improved_euler', + h => 0.2, + n => 5 + ); + is $soln->{x}, [ 0, 0.2, 0.4, 0.6, 0.8, 1.0 ], 'returns correct x'; + # check the following to 6 digits. + is roundArray($soln->{k1}), + roundArray([ undef, 1, 1.597504, 2.947084257, 6.662185892, 22.89372811 ]), + 'returns correct k1'; + is roundArray($soln->{k2}), + roundArray([ undef, 1.48, 2.617058758, 5.462507804, 15.40751657, 87.41805808 ]), + 'returns correct k2'; + is roundArray($soln->{y}), + roundArray([ 1, 1.248, 1.669456276, 2.510415482, 4.717385728, 15.74856435 ]), + 'returns correct y'; +}; + +subtest "Solve an ODE using Heun's method" => sub { + my $g = sub { + my ($x, $y) = @_; + return $x**2 + $y**2; + }; + + my $soln = solveDiffEqn( + $g, 1, + x0 => 0, + method => 'heun', + h => 0.2, + n => 5 + ); + is $soln->{x}, [ 0, 0.2, 0.4, 0.6, 0.8, 1.0 ], 'returns correct x'; + # check the following to 6 digits. + is roundArray($soln->{k1}), + roundArray([ undef, 1.0, 1.5908551111111113, 2.9161500566582608, 6.502422880077087, 21.460193376361623 ]), + 'returns correct k1'; + is roundArray($soln->{k2}), + roundArray([ + undef, 1.302222222222222, 2.235263883735181, 4.482786757206292, 11.72935117869894, 55.9909574019759 ]), + 'returns correct k2'; + is roundArray($soln->{y}), + roundArray([ + 1, 1.2453333333333334, 1.6601656714491662, 2.478391187863023, 4.562915008671718, 14.034568287786184 ]), + 'returns correct y'; +}; + +subtest 'Solve an ODE using 4th order Runge-Kutta ' => sub { + my $g = sub { + my ($x, $y) = @_; + return $x**2 + $y**2; + }; + + my $soln = solveDiffEqn($g, 1, method => 'rk4', h => 0.2, n => 5); + is $soln->{x}, [ 0, 0.2, 0.4, 0.6, 0.8, 1.0 ], 'returns correct x'; + # check the following to 6 digits. + is roundArray($soln->{k1}), + roundArray([ undef, 1, 1.6099859, 3.0361440, 7.3407438, 34.1115788 ]), + 'returns correct k1'; + is roundArray($soln->{k2}), + roundArray([ undef, 1.22000, 2.0893660, 4.2481371, 11.8886191, 85.3878304 ]), + 'returns correct k2'; + is roundArray($soln->{k3}), + roundArray([ undef, 1.2688840, 2.2272318, 4.7475107, 15.1663436, 205.9940166 ]), + 'returns correct k3'; + is roundArray($soln->{k4}), + roundArray([ undef, 1.6119563, 3.0446888, 7.3582574, 32.8499206, 2208.5212543 ]), + 'returns correct k4'; + is roundArray($soln->{y}), + roundArray([ 1, 1.25299088, 1.6959198, 2.6421097, 5.7854627, 99.9653469 ]), + 'returns correct y'; +}; + +subtest 'Test that errors of the bisection method are returned correctly' => sub { + my $bisect = bisection(Formula('x^2+2'), [ 0, 1 ]); + like $bisect->{error}, qr/The function must be a code reference/, 'The function is not a code reference'; + + my $g = sub { return (shift)**2 - 2; }; + + $bisect = bisection($g, [ 0, 1 ]); + like $bisect->{error}, qr/The function may not have a root/, 'The function may not have a root'; + + $bisect = bisection($g, [ 0, 1, 2 ]); + is $bisect->{error}, 'The interval must be an array ref of length 2', 'The interval must be an array ref'; + + $bisect = bisection($g, [ 1, 0 ]); + is $bisect->{error}, 'The initial interval [a, b] must satisfy a < b', 'Check the initial interval for a < b'; + + $bisect = bisection($g, [ 0, 2 ], eps => -1); + is $bisect->{error}, 'The option eps must be a positive number', 'The option eps must be a positive number'; + + $bisect = bisection($g, [ 0, 2 ], max_iter => -1); + is $bisect->{error}, 'The option max_iter must be a positive integer', + 'The option max_iter must be a positive integer'; + + $bisect = bisection($g, [ 0, 2 ], max_iter => 1.5); + is $bisect->{error}, 'The option max_iter must be a positive integer', + 'The option max_iter must be a positive integer'; + + $bisect = bisection(sub { (shift)**2 - 19 }, [ 0, 100 ], max_iter => 20); + like $bisect->{error}, qr/You have reached the maximum/, 'Reached the maximum number of iterations.'; +}; + +subtest 'Find a root via bisection' => sub { + my $g = sub { return (shift)**2 - 2; }; + + my $bisect = bisection($g, [ 0, 2 ]); + is roundArray([ map { $_->[0] } @{ $bisect->{intervals} }[ 0 .. 10 ] ]), + roundArray([ 0.0, 1.0, 1.0, 1.25, 1.375, 1.375, 1.40625, 1.40625, 1.4140625, 1.4140625, 1.4140625 ]), + 'left endpoints of the bisection method'; + is roundArray([ map { $_->[1] } @{ $bisect->{intervals} }[ 0 .. 10 ] ]), + roundArray([ 2.0, 2.0, 1.5, 1.5, 1.5, 1.4375, 1.4375, 1.421875, 1.421875, 1.41796875, 1.416015625 ]), + 'right endpoints of the bisection method'; + is sqrt(2), float($bisect->{root}, precision => 6), 'The root was found successfully.'; +}; + +subtest "Test that the errors from Newton's method" => sub { + my $newton = newton(Formula('x^2+2'), 1); + like $newton->{error}, qr/The function must be a code reference/, 'The function is not a code reference'; + + my $g = sub { return (shift)**2 - 2; }; + my $dg = sub { return 2 * (shift); }; + + $newton = newton($g, $dg, 1, eps => -1e-8); + like $newton->{error}, qr/The option eps must be a positive number/, 'The option eps must be a positive number'; + + $newton = newton($g, $dg, 1, feps => -1e-8); + like $newton->{error}, qr/The option feps must be a positive number/, + 'The option feps must be a positive number'; + + $newton = newton($g, $dg, 1, max_iter => -10); + like $newton->{error}, qr/The option max_iter must be a positive integer/, + 'The option max_iter must be a positive number'; + + $newton = newton(sub { my $x = shift; ($x)**2 + 2 }, sub { my $x = shift; 2 * $x; }, 1); + like $newton->{error}, qr/Newton's method did not converge in \d+ steps/, "Newton's method did not converge."; +}; + +subtest "Find a root using Newton's method" => sub { + my $g = sub { return (shift)**2 - 2; }; + my $dg = sub { return 2 * (shift); }; + + my $newton = newton($g, $dg, 10); + is sqrt(2), float($newton->{root}), 'The root was found successfully.'; + + is roundArray([ @{ $newton->{iterations} }[ 0 .. 5 ] ]), + roundArray([ 10.0, 5.1, 2.7460784313725486, 1.7371948743795982, 1.444238094866232, 1.4145256551487377 ]), + "iterations of newton's method"; + +}; + +subtest 'Test that the errors from the Secant method' => sub { + + my $secant = secant(Formula('x^2+2'), 1, 2); + like $secant->{error}, qr/The function must be a code reference/, 'The function is not a code reference'; + + my $g = sub { return (shift)**2 - 2; }; + + $secant = secant($g, 1, 2, eps => -1e-8); + like $secant->{error}, qr/The option eps must be a positive number/, 'The option eps must be a positive number'; + + $secant = secant($g, 1, 2, feps => -1e-8); + like $secant->{error}, qr/The option feps must be a positive number/, + 'The option feps must be a positive number'; + + $secant = secant($g, 1, 2, max_iter => -10); + like $secant->{error}, qr/The option max_iter must be a positive integer/, + 'The option max_iter must be a positive number'; + + $secant = secant(sub { return (shift)**2 + 2; }, 1, 2); + like $secant->{error}, qr/The secant method did not converge in \d+ steps/, + 'The secant method did not converge.'; +}; + +subtest 'Find a root using the Secant method' => sub { + my $g = sub { return (shift)**2 - 2; }; + my $secant = secant($g, 1, 2); + is sqrt(2), float($secant->{root}), 'The root was found successfully.'; + + is roundArray([ @{ $secant->{iterations} }[ 0 .. 6 ] ]), + roundArray([ 1.0, 2.0, 1.3333333333333335, 1.4, 1.4146341463414633, 1.41421143847487, 1.4142135620573204 ]), + 'iterations of the secant method'; + +}; + sub roundArray { my ($arr, %options) = @_; %options = (digits => 6, %options);