From b42c85c02c7edc37aa0faec207bbf54e46d22e40 Mon Sep 17 00:00:00 2001 From: Peter Staab Date: Wed, 6 Mar 2024 18:57:08 -0500 Subject: [PATCH] Add additional methods and tests for numerical methods. This adds methods for newtonDividedDifference, Legendre Polynomials, Newton-Cotes quadrature, Gaussian Quadrature, numerical differential equation solvers, and rootfinding techniques. u --- macros/math/PGnumericalmacros.pl | 666 ++++++++++++++++++++++++++++++- t/macros/numerical_methods.t | 360 +++++++++++++++++ 2 files changed, 1025 insertions(+), 1 deletion(-) diff --git a/macros/math/PGnumericalmacros.pl b/macros/math/PGnumericalmacros.pl index 921fae7f28..a3ddb023e5 100644 --- a/macros/math/PGnumericalmacros.pl +++ b/macros/math/PGnumericalmacros.pl @@ -393,6 +393,177 @@ sub javaScript_cubic_spline { return $output_str; } +=head3 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 arrayref of arrayrefs 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 + + f[x_i,x_j] = (f[x_j]-f[x_i])/(x_j - x_i) + +=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; +} + +=head3 legendreP + +Returns a code reference to the Legendre Polynomial of degree C. + +Usage: + + $poly = legendreP($n) + +And then evaluations can be found with C<&$poly(0.5)> for example to evaluate the polynomial at +C. Even though this is a polynomial, the standard domain of these are [-1,1], although 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; + }; +} + +=head3 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); + }; +} + +=head3 legendreP_nodes_weights + +Finds the nodes (roots) and weights of the Legendre Polynomials of degree C. These are used in +Gaussian Quadrature. + +Usage: + + ($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); +} + =head2 Numerical Integration methods =head3 lefthandsum @@ -604,6 +775,123 @@ sub inv_romberg { return $b; } +=head3 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), C<$a> and C<$b> is the interval C<[$a,$b]>. + +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, the number indicates the number of used nodes for the formula. + +=item n + +This number 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; +} + +=head3 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) + +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. + +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; +} + =head2 Differential Equation Methods =head3 rungeKutta4 @@ -650,7 +938,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 ]); @@ -672,4 +960,380 @@ sub rungeKutta4 { } } +=head3 solveDiffEqn + +Produces a numerical solution to the differential equation y'=f(x,y) using a number of optional methods. + +B + +=over + +=item * C an 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-values number for the initial point + +=back + +B + +=over + +=item * C the initial x value (defaults to 0) + +=item * C the stepsize 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 ref of the x values which are C + +=item * C an array ref of the y values (depending on the method used) + +=item * C the intermediate function values used (depending 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 ref 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 + }; + } +} + +=head2 Rootfinding + +=head3 bisection + +Performs the bisection method for the function C<$f> and initial interval C<$int> (arrayref). +An example is + + $f = sub { $x = shift; $x**2-2;} + $bisect = bisection($f, [1, 2]); + +The result is a hash with fields root (the estimated root), intervals (an array ref or +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 ref of 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 arrayref of the intervals (each interval also an array ref) + +=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 + }; +} + +=head3 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 arrayref +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 arrayref 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 }; +} + +=head3 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 arrayref +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. Default is C<20>. + +=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 arrayref of 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 e0e36f4f43..0cef9b66f7 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]'; @@ -105,6 +120,122 @@ 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 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 +257,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);