diff --git a/macros/math/PGnumericalmacros.pl b/macros/math/PGnumericalmacros.pl index 8cbc7eee21..0362818754 100644 --- a/macros/math/PGnumericalmacros.pl +++ b/macros/math/PGnumericalmacros.pl @@ -17,6 +17,12 @@ =head1 NAME Numerical methods for the PG language +=cut + +BEGIN { strict->import; } + +sub _PGnumericalmacros_init { } + =head2 Interpolation methods =head3 Plotting a list of points (piecewise linear interpolation) @@ -34,9 +40,8 @@ =head3 Plotting a list of points (piecewise linear interpolation) =cut BEGIN { strict->import; } - -sub _PGnumericalmacros_init { -} +# TODO: this fails if the $xref are not unique and ordered. Should check for +# these and document. sub plot_list { my ($xref, $yref) = @_; @@ -51,8 +56,7 @@ sub plot_list { my (@x_vals, @y_vals); unless (defined($yref)) { #with only one entry we assume (x0,y0,x1,y1..); if (@$xref % 2 == 1) { - die "ERROR in plot_list -- single array of input has odd number of - elements"; + die "ERROR in plot_list -- single array of input has odd number of elements"; } my @in = @$xref; @@ -64,7 +68,7 @@ sub plot_list { $yref = \@y_vals; } - my $fun = sub { + return sub { my $x = shift; my $y; my ($x0, $x1, $y0, $y1); @@ -85,27 +89,37 @@ sub plot_list { } $y; }; - $fun; } -=head2 Horner polynomial/ Newton polynomial +=head3 Horner polynomial/ Newton polynomial Usage: - $fn = horner([x0,x1,x2],[q0,q1,q2]); + $fn = horner([x0,x1,x2, ...],[q0,q1,q2, ...]); Produces the newton polynomial - &$fn(x) = q0 + q1*(x-x0) +q2*(x-x1)*(x-x0); + &$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. + +=head4 Example -Generates a subroutine which evaluates a polynomial passing through the points C<(x0,q0), (x1,q1), -... > using Horner's method. + $h = horner([0,1,2],[1,-1,2]); + +Then C<&$h(num)> returns the polynomial at the value C. For example, +C<&$h(1.5)=1>. =cut sub horner { my ($xref, $qref) = @_; # get the coefficients - my $fn = sub { + die 'The x inputs and q inputs must be the same length' + unless scalar(@$xref) == scalar(@$qref); + return sub { my $x = shift; my @xvals = @$xref; my @qvals = @$qref; @@ -116,14 +130,13 @@ sub horner { } $y; }; - $fn; } -=head2 Hermite polynomials +=head3 Hermite polynomials Usage: - $poly = hermit([x0,x1...],[y0,y1...],[yp0,yp1,...]); + $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 @@ -133,14 +146,23 @@ =head2 Hermite polynomials The polynomial will be of high degree and may wobble unexpectedly. Use the Hermite splines described below and in Hermite.pm for most graphing purposes. +=head4 Example + + $h = hermite([0,1],[0,0],[1,-1]); + +C<&$h(num)> will evaluate the hermite polynomial at C. For example, C<&$h(0.5)> +returns C<0.25>. + =cut sub hermite { my ($x_ref, $y_ref, $yp_ref) = @_; + die 'The input array refs all must be the same length' + unless scalar(@$x_ref) == scalar(@$y_ref) && scalar(@$y_ref) == scalar(@$yp_ref); my (@zvals, @qvals); my $n = $#{$x_ref}; - my $i = 0; - foreach $i (0 .. $n) { + + for my $i (0 .. $n) { $zvals[ 2 * $i ] = $$x_ref[$i]; $zvals[ 2 * $i + 1 ] = $$x_ref[$i]; $qvals[ 2 * $i ][0] = $$y_ref[$i]; @@ -149,22 +171,22 @@ sub hermite { $qvals[ 2 * $i ][1] = ($qvals[ 2 * $i ][0] - $qvals[ 2 * $i - 1 ][0]) / ($zvals[ 2 * $i ] - $zvals[ 2 * $i - 1 ]) unless $i == 0; - } - my $j; - foreach $i (2 .. (2 * $n + 1)) { - foreach $j (2 .. $i) { + + for my $i (2 .. (2 * $n + 1)) { + for my $j (2 .. $i) { $qvals[$i][$j] = ($qvals[$i][ $j - 1 ] - $qvals[ $i - 1 ][ $j - 1 ]) / ($zvals[$i] - $zvals[ $i - $j ]); } } + my @output; - foreach $i (0 .. 2 * $n + 1) { + for my $i (0 .. 2 * $n + 1) { push(@output, $qvals[$i][$i]); } - horner(\@zvals, \@output); + return horner(\@zvals, \@output); } -=head2 Hermite splines +=head3 Hermite splines Usage @@ -203,15 +225,18 @@ sub hermite_spline { $yp0 = $yp1; } - my $hermite_spline_sub = sub { + return sub { my $x = shift; my $y; my $fun; my @xvals = @$xref; my @fns = @polys; - return $y = &{ $fns[0] }($x) if $x == $xvals[0]; #handle left most endpoint - while (@xvals && $x > $xvals[0]) { # find the function for this range of x + #handle left most endpoint + return $y = &{ $fns[0] }($x) if $x == $xvals[0]; + + # find the function for this range of x + while (@xvals && $x > $xvals[0]) { shift(@xvals); $fun = shift(@fns); } @@ -221,12 +246,11 @@ sub hermite_spline { if (@xvals && defined($fun)) { $y = &$fun($x); } - $y; + return $y; }; - $hermite_spline_sub; } -=head2 Cubic spline approximation +=head3 Cubic spline approximation Usage: @@ -263,7 +287,7 @@ =head2 Cubic spline approximation sub cubic_spline { my ($t_ref, $y_ref) = @_; my @spline_coeff = create_cubic_spline($t_ref, $y_ref); - sub { + return sub { my $x = shift; eval_cubic_spline($x, @spline_coeff); } @@ -274,54 +298,49 @@ sub eval_cubic_spline { # The knot points given by $t_ref should be in order. my $i = 0; my $out = 0; - while (defined($t_ref->[ $i + 1 ]) && $x > $t_ref->[ $i + 1 ]) { - - $i++; - } + $i++ while (defined($t_ref->[ $i + 1 ]) && $x > $t_ref->[ $i + 1 ]); unless (defined($t_ref->[$i]) && ($t_ref->[$i] <= $x) && ($x <= $t_ref->[ $i + 1 ])) { $out = undef; # input value is not in the range defined by the function. } else { - $out = ($t_ref->[ $i + 1 ] - $x) * (($d_ref->[$i]) + ($a_ref->[$i]) * ($t_ref->[ $i + 1 ] - $x)**2) + + $out = + ($t_ref->[ $i + 1 ] - $x) * (($d_ref->[$i]) + ($a_ref->[$i]) * ($t_ref->[ $i + 1 ] - $x)**2) + ($x - $t_ref->[$i]) * (($b_ref->[$i]) * ($x - $t_ref->[$i])**2 + ($c_ref->[$i])); - } - $out; + return $out; } -########################################################################## -# Cubic spline algorithm adapted from p319 of Kincaid and Cheney's Numerical Analysis -########################################################################## - +# Cubic spline algorithm adapted from p319 of Kincaid and Cheney's Numerical Analysis. sub create_cubic_spline { my ($t_ref, $y_ref) = @_; # The knot points are given by $t_ref (in order) and the function values by $y_ref my $num = $#{$t_ref}; # number of knots - my $i; + my (@h, @b, @u, @v, @z); - foreach $i (0 .. $num - 1) { + for my $i (0 .. $num - 1) { $h[$i] = $t_ref->[ $i + 1 ] - $t_ref->[$i]; $b[$i] = 6 * ($y_ref->[ $i + 1 ] - $y_ref->[$i]) / $h[$i]; } $u[1] = 2 * ($h[0] + $h[1]); $v[1] = $b[1] - $b[0]; - foreach $i (2 .. $num - 1) { + for my $i (2 .. $num - 1) { $u[$i] = 2 * ($h[$i] + $h[ $i - 1 ]) - ($h[ $i - 1 ])**2 / $u[ $i - 1 ]; - $v[$i] = $b[$i] - $b[ $i - 1 ] - $h[ $i - 1 ] * $v[ $i - 1 ] / $u[ $i - 1 ]; + $v[$i] = + $b[$i] - $b[ $i - 1 ] - $h[ $i - 1 ] * $v[ $i - 1 ] / $u[ $i - 1 ]; } $z[$num] = 0; - for ($i = $num - 1; $i > 0; $i--) { + for (my $i = $num - 1; $i > 0; $i--) { $z[$i] = ($v[$i] - $h[$i] * $z[ $i + 1 ]) / $u[$i]; } $z[0] = 0; my ($a_ref, $b_ref, $c_ref, $d_ref); - foreach $i (0 .. $num - 1) { + for my $i (0 .. $num - 1) { $a_ref->[$i] = $z[$i] / (6 * $h[$i]); $b_ref->[$i] = $z[ $i + 1 ] / (6 * $h[$i]); $c_ref->[$i] = (($y_ref->[ $i + 1 ]) / $h[$i] - $z[ $i + 1 ] * $h[$i] / 6); $d_ref->[$i] = (($y_ref->[$i]) / $h[$i] - $z[$i] * $h[$i] / 6); } - ($t_ref, $a_ref, $b_ref, $c_ref, $d_ref); + return ($t_ref, $a_ref, $b_ref, $c_ref, $d_ref); } sub javaScript_cubic_spline { @@ -387,9 +406,71 @@ sub javaScript_cubic_spline { $output_str; } +=head3 Newton Divided Difference + +Computes the newton divided difference table with the function C. + +=head4 Arguments + +=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 + +=head4 Ouput + +An arrayref of arrayrefs of Divided Differences. + +=head4 Examples + + $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; +} + =head2 Numerical Integration methods -=head3 Integration by Left Hand Sum +=head3 Left Hand Riemann Sum Usage: @@ -402,24 +483,20 @@ =head3 Integration by Left Hand Sum =cut sub lefthandsum { - my $fn_ref = shift; - my $x0 = shift; - my $x1 = shift; - my %options = @_; + my ($fn_ref, $x0, $x1, %options) = @_; assign_option_aliases(\%options, intervals => 'steps',); - set_default_options(\%options, steps => 30,); + set_default_options(\%options, steps => 30); my $steps = $options{steps}; my $delta = ($x1 - $x0) / $steps; - my $i; - my $sum = 0; + my $sum = 0; - foreach $i (0 .. ($steps - 1)) { - $sum = $sum + &$fn_ref($x0 + ($i) * $delta); + for my $i (0 .. ($steps - 1)) { + $sum += &$fn_ref($x0 + $i * $delta); } - $sum * $delta; + return $sum * $delta; } -=head3 Integration by Right Hand Sum +=head3 Right Hand Riemann Sum Usage: @@ -432,89 +509,75 @@ =head3 Integration by Right Hand Sum =cut sub righthandsum { - my $fn_ref = shift; - my $x0 = shift; - my $x1 = shift; - my %options = @_; + my ($fn_ref, $x0, $x1, %options) = @_; assign_option_aliases(\%options, intervals => 'steps',); - set_default_options(\%options, steps => 30,); + set_default_options(\%options, steps => 30); my $steps = $options{steps}; my $delta = ($x1 - $x0) / $steps; - my $i; - my $sum = 0; + my $sum = 0; - foreach $i (1 .. ($steps)) { - $sum = $sum + &$fn_ref($x0 + ($i) * $delta); + for my $i (1 .. $steps) { + $sum += &$fn_ref($x0 + $i * $delta); } - $sum * $delta; + return $sum * $delta; } -=head3 Integration by Midpoint rule +=head3 Midpoint rule Usage: midpoint(function_reference, start, end, steps=>30 ); -Implements the Midpoint rule using 30 intervals between 'start' and 'end'. +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. =cut sub midpoint { - my $fn_ref = shift; - my $x0 = shift; - my $x1 = shift; - my %options = @_; + my ($fn_ref, $x0, $x1, %options) = @_; assign_option_aliases(\%options, intervals => 'steps',); set_default_options(\%options, steps => 30,); my $steps = $options{steps}; my $delta = ($x1 - $x0) / $steps; - my $i; - my $sum = 0; + my $sum = 0; - foreach $i (0 .. ($steps - 1)) { - $sum = $sum + &$fn_ref($x0 + ($i + 1 / 2) * $delta); + for my $i (0 .. ($steps - 1)) { + $sum += &$fn_ref($x0 + ($i + 1 / 2) * $delta); } - $sum * $delta; + return $sum * $delta; } -=head3 Integration by Simpson's rule +=head3 Simpson's rule Usage: simpson(function_reference, start, end, steps=>30 ); -Implements Simpson's rule using 30 intervals between 'start' and 'end'. +Implements Simpson's rule between 'start' and 'end'. The first three arguments are required. The final argument (number of steps) is optional and defaults to 30, but must be even. =cut sub simpson { - my $fn_ref = shift; - my $x0 = shift; - my $x1 = shift; - my %options = @_; + my ($fn_ref, $x0, $x1, %options) = @_; assign_option_aliases(\%options, intervals => 'steps',); - set_default_options(\%options, steps => 30,); + set_default_options(\%options, steps => 30); my $steps = $options{steps}; - unless ($steps % 2 == 0) { - die "Error: Simpson's rule requires an even number of steps."; - } + die "Error: Simpson's rule requires an even number of steps." unless $steps % 2 == 0; my $delta = ($x1 - $x0) / $steps; - my $i; - my $sum = 0; - for ($i = 1; $i < $steps; $i = $i + 2) { # look this up - loop by two. - $sum = $sum + 4 * &$fn_ref($x0 + $i * $delta); + my $sum = 0; + for (my $i = 1; $i < $steps; $i = $i + 2) { # look this up - loop by two. + $sum += 4 * &$fn_ref($x0 + $i * $delta); } - for ($i = 2; $i < $steps - 1; $i = $i + 2) { # ditto - $sum = $sum + 2 * &$fn_ref($x0 + $i * $delta); + for (my $i = 2; $i < $steps - 1; $i = $i + 2) { # ditto + $sum += 2 * &$fn_ref($x0 + $i * $delta); } - $sum = $sum + &$fn_ref($x0) + &$fn_ref($x1); - $sum * $delta / 3; + $sum += &$fn_ref($x0) + &$fn_ref($x1); + return $sum * $delta / 3; } -=head3 Integration by trapezoid rule +=head3 trapezoid rule Usage: @@ -527,21 +590,17 @@ =head3 Integration by trapezoid rule =cut sub trapezoid { - my $fn_ref = shift; - my $x0 = shift; - my $x1 = shift; - my %options = @_; + my ($fn_ref, $x0, $x1, %options) = @_; assign_option_aliases(\%options, intervals => 'steps',); - set_default_options(\%options, steps => 30,); + set_default_options(\%options, steps => 30); my $steps = $options{steps}; my $delta = ($x1 - $x0) / $steps; - my $i; - my $sum = 0; + my $sum = 0; - foreach $i (1 .. ($steps - 1)) { - $sum = $sum + &$fn_ref($x0 + $i * $delta); + for my $i (1 .. ($steps - 1)) { + $sum += &$fn_ref($x0 + $i * $delta); } - $sum = $sum + 0.5 * (&$fn_ref($x0) + &$fn_ref($x1)); + $sum += 0.5 * (&$fn_ref($x0) + &$fn_ref($x1)); $sum * $delta; } @@ -556,32 +615,18 @@ =head3 Romberg method of integration =cut sub romberg { - my $fn_ref = shift; - my $x0 = shift; - my $x1 = shift; - my %options = @_; - set_default_options(\%options, level => 6,); + my ($fn_ref, $x0, $x1, %options) = @_; + set_default_options(\%options, level => 6); my $level = $options{level}; romberg_iter($fn_ref, $x0, $x1, $level, $level); } sub romberg_iter { - my $fn_ref = shift; - my $x0 = shift; - my $x1 = shift; - my $j = shift; - my $k = shift; - my $out; - if ($k == 1) { - $out = trapezoid($fn_ref, $x0, $x1, steps => 2**($j - 1)); - } else { - - $out = - (4**($k - 1) * romberg_iter($fn_ref, $x0, $x1, $j, $k - 1) - - romberg_iter($fn_ref, $x0, $x1, $j - 1, $k - 1)) / - (4**($k - 1) - 1); - } - $out; + my ($fn_ref, $x0, $x1, $j, $k) = @_; + return $k == 1 + ? trapezoid($fn_ref, $x0, $x1, steps => 2**($j - 1)) + : (4**($k - 1) * romberg_iter($fn_ref, $x0, $x1, $j, $k - 1) - romberg_iter($fn_ref, $x0, $x1, $j - 1, $k - 1)) + / (4**($k - 1) - 1); } =head3 Inverse Romberg @@ -594,19 +639,23 @@ =head3 Inverse Romberg 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. +=head4 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. + + $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. + =cut sub inv_romberg { - my $fn_ref = shift; - my $a = shift; - my $value = shift; - my $b = $a; - my $delta = 1; - my $i = 0; - my $funct; - my $deriv; - - while (abs($delta) > 0.000001 && $i < 5000) { + my ($fn_ref, $a, $value) = @_; + my ($b, $delta, $i, $funct, $deriv) = ($a, 1, 0); + + while (abs($delta) > 0.000001 && $i++ < 5000) { $funct = romberg($fn_ref, $a, $b) - $value; $deriv = &$fn_ref($b); if ($deriv == 0) { @@ -614,35 +663,259 @@ sub inv_romberg { return; } $delta = $funct / $deriv; - $b = $b - $delta; - $i++; + $b -= $delta; } if ($i == 5000) { warn 'Newtons method does not converge.'; return; } - $b; + return $b; +} + +=head3 Newton Cotes functions. + +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]>. + +=head4 options + +=over + +=item method + +The method options are either open or closed methods. The closed newton-cotes formula methods +are trapezoid, simpson, three-eighths, boole. The open newton-cotes formula methods are +open1, open2, open3, open4, 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 Legendre Polynomials + +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 derivative of Legendre Polynomials + +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 Nodes and Weights of Legendre Polynomial + +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) + +=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); +} + +=head3 Gaussian Quadrature + +Compute the integral of a function C<$f> on an interval C<[a,b]> using Gassian +Quadrature. + +Usage: + + gauss_quad($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 gauss_quad { + 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 +=head3 4th-order Runge-Kutta -Finds integral curve of a vector field using the 4th order Runge Kutta method. +Finds integral curve of a vector field using the 4th order Runge Kutta method by +providing the function C Usage: rungeKutta4( &vectorField(t,x),%options); - Returns: \@array of points [t,y}) + Returns: \@array of points [t,y] Default %options: - 'initial_t' => 1, - 'initial_y' => 1, - 'dt' => .01, - 'num_of_points' => 10, #number of reported points - 'interior_points' => 5, # number of 'interior' steps between reported points - 'debug' + 'initial_t' =>1, + 'initial_y' => 1, + 'dt' => .01, + 'num_of_points' => 10, #number of reported points + 'interior_points' => 5, # number of 'interior' steps between reported points + 'debug' =cut @@ -668,22 +941,21 @@ sub rungeKutta4 { my $rf_rhs = sub { my @in = @_; 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); + $errors .= " $err at ( " . join(" , ", @in) . " )
\n" + if defined($err); + $out = 'NaN' if defined($err) and not is_a_number($out); $out; }; my @output = ([ $t, $y ]); - my ($i, $j, $K1, $K2, $K3, $K4); - - for ($j = 0; $j < $num; $j++) { - for ($i = 0; $i < $num2; $i++) { - $K1 = $dt * &$rf_rhs($t, $y); - $K2 = $dt * &$rf_rhs($t + $dt / 2, $y + $K1 / 2); - $K3 = $dt * &$rf_rhs($t + $dt / 2, $y + $K2 / 2); - $K4 = $dt * &$rf_rhs($t + $dt, $y + $K3); - $y = $y + ($K1 + 2 * $K2 + 2 * $K3 + $K4) / 6; - $t = $t + $dt; + for my $j (0 .. $num - 1) { + for my $i (0 .. $num2 - 1) { + my $K1 = $dt * &$rf_rhs($t, $y); + my $K2 = $dt * &$rf_rhs($t + $dt / 2, $y + $K1 / 2); + my $K3 = $dt * &$rf_rhs($t + $dt / 2, $y + $K2 / 2); + my $K4 = $dt * &$rf_rhs($t + $dt, $y + $K3); + $y += ($K1 + 2 * $K2 + 2 * $K3 + $K4) / 6; + $t += $dt; } push(@output, [ $t, $y ]); } @@ -694,5 +966,390 @@ sub rungeKutta4 { } } +=head3 Robust Differential Equation Solver + +Produces a numerical solution to the differential equation y'=f(x,y) with the +function C. + +=head4 Arguments + +=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 + +=head4 Options + +=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 + +=head4 Output + +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 + +=head4 Examples + +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. + +=head4 Arguments + +=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 + +=head4 Options + +=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 + +=head4 Output + +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. + +=head4 Arguments + +=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 + +=head4 Options + +=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 + +=head4 Output + +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. + +=head4 Arguments + +=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 + +=head4 Options + +=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 + +=head4 Output + +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.t b/t/macros/numerical.t new file mode 100644 index 0000000000..9cfd7c7e2c --- /dev/null +++ b/t/macros/numerical.t @@ -0,0 +1,496 @@ +#!/usr/bin/env perl + +# Tests subroutines in the PGnumericamacros.pl macro. + +use Test2::V0 '!E', { E => 'EXISTS' }; + +die "PG_ROOT not found in environment.\n" unless $ENV{PG_ROOT}; +do "$ENV{PG_ROOT}/t/build_PG_envir.pl"; + +loadMacros('PGnumericalmacros.pl', 'MathObjects.pl', 'PGauxiliaryFunctions.pl'); + +subtest 'plot_list' => sub { + ok my $p1 = plot_list([ 0, 0, 1, 2 ]); + is &$p1(0.75), 1.5, 'linear interpolation at $x=0.75'; + is &$p1(0.25), 0.5, 'linear interpolation at $x=0.25'; + + ok my $p2 = plot_list([ (0, 0), (1, 2) ]); + is &$p2(0.75), 1.5, 'linear interpolation at $x=0.75'; + is &$p2(0.25), 0.5, 'linear interpolation at $x=0.25'; + + ok my $p3 = plot_list([ 0, 3 ], [ 4, 0 ]); + is &$p3(1.5), 2, 'linear interpolation at $x=0.75'; + is &$p3(2), 4 / 3, 'linear interpolation at $x=0.25'; + + like dies { plot_list([ 0, 1, 3, 4, 5 ]) }, + qr/single array of input has odd number/, + 'Input of odd number of values'; + like dies { plot_list(0, 1, 3, 4, 5) }, + qr/Error in plot_list:X values must be given as an array reference./, + 'Values are not given as an array reference'; +}; + +subtest 'horner' => sub { + ok my $h1 = horner([ 0, 1, 2 ], [ 1, -1, 2 ]); # 1-1*(x-0)+2(x-0)*(x-1) + is &$h1(0.5), 0, 'h1(0.5)=0'; #1-1*0.5+2*(0.5)*(-0.5) = 0 + is &$h1(1.5), 1, 'h1(1.5)=1'; # 1-1*(1.5)+2*(1.5)*(0.5)= 1 + + ok my $h2 = horner([ -1, 1, 2, 5 ], [ 2, 0, -2, 1 ]); # 2+0(x+1)-2(x+1)(x-1)+(x+1)(x-1)(x-2) + is &$h2(0), 6, 'h2(0)=6'; # 2-2(1)(-1)+(1)(-1)(-2) = 6 + is &$h2(3), -6, 'h2(3)=-6'; # 2-2(4)(2)+(4)(2)(1) = -6 + + like dies { horner([ 0, 1, 2 ], [ -1, 0, 2, 3 ]); }, + qr/The x inputs and q inputs must be the same length/, + 'Input array refs are different lengths.'; +}; + +subtest 'hermite' => sub { + ok my $h1 = hermite([ 0, 1 ], [ 0, 0 ], [ 1, -1 ]); # x-x^2 + is &$h1(0), 0, 'h1(0)=0'; + is &$h1(1), 0, 'h1(1)=0'; + is &$h1(0.5), 0.25, 'h1(0.5)=0.25'; + is &$h1(0.25), 0.1875, 'h1(0.25)=0.1875'; + + ok my $h2 = hermite([ 0, 1, 3 ], [ 2, 0, 1 ], [ 1, 0, -1 ]); + is &$h2(0), 2, 'h2(0)=2'; + is &$h2(1), 0, 'h2(1)=0'; + is Round(&$h2(3), 10), 1, 'h2(3)=1'; + is Round(&$h2(0.5), 10), Round(1573 / 1728, 10), 'h2(1/2)=1573/1728'; + is Round(&$h2(2), 10), Round(55 / 27, 10), 'h2(2)=55/27'; + + like dies { hermite([ 0, 1, 2 ], [ 1, 1, 1 ], [ 0, 2 ]) }, + qr/The input array refs all must be the same length/, + 'Input array refs are different lengths.'; +}; + +subtest 'hermite spline' => sub { + ok my $h = hermite_spline([ 0, 1, 3 ], [ 3, 1, -5 ], [ 1, -2, 0 ]); + is &$h(0), 3, 'h(0)=3'; + is &$h(1), 1, 'h(1)=1'; + is &$h(3), -5, 'h(3)=-5'; + is &$h(0.5), 19 / 8, 'h(1/2)=19/8'; + is &$h(2), -2.5, 'h(2)=-2.5'; + is &$h(1.3), 0.202, 'h(1.3)=0.202'; +}; + +subtest 'cubic spline' => sub { + ok my $s = cubic_spline([ 0, 1, 2 ], [ 0, 1, 0 ]); + is &$s(0), 0, 's(0)=0'; + is &$s(1), 1, 's(1)=1'; + is &$s(2), 0, 's(2)=0'; + # not sure how the spline is constructed to find intermediate points. + # is &$s(0.25), 0.4375; + # is &$s(0.5), 0.75; +}; + +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]'; + is righthandsum($f, 0, 2, steps => 4), 3.75, 'right hand sum of x^2 on [0,2]'; + is midpoint($f, 0, 2, steps => 4), 2.625, 'midpoint rule of x^2 on [0,2]'; +}; + +subtest 'Quadrature' => sub { + my $f = sub { my $x = shift; return $x * $x; }; + my $g = sub { my $x = shift; return exp($x); }; + 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); }, + qr /Error: Simpson's rule requires an even number of steps./, + 'Check for odd number of steps'; + + is trapezoid($f, 0, 2, steps => 4), 2.75, 'Trapezoid rule of x^2 on [0,2]'; + + is romberg($f, 0, 2), 8 / 3, 'Romberg interation for x^2 on [0,2]'; + 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(gauss_quad($f), 15), 0, 'gauss_quad(x^3) on [-1,1]'; + is Round(gauss_quad($f, a => 0, b => 1), 15), 0.25, 'gauss_quad(x^3) on [0,1]'; + + is gauss_quad($f, n => 2, a => 0, b => 1), 0.25, 'gauss_quad(x^3, n=>2) on [0,1]'; + + my $g = sub { my $x = shift; return $x**6; }; + is gauss_quad($g), 2 / 7, 'gauss_quad(x^6) on [-1,1]'; + is Round(gauss_quad($g, n => 2), 15), Round(2 * (1 / sqrt(3))**6, 15), 'gauss_quad(x^6) on [-1,1]'; + + my $e_x = sub { my $x = shift; return exp($x); }; + is Round(gauss_quad($e_x, n => 3), 15), Round(5 * (exp(-sqrt(3 / 5)) + exp(sqrt(3 / 5))) / 9 + 8 / 9, 15), + 'gauss_quad(x^6) on [-1,1]'; + is Round(gauss_quad($e_x, n => 15, a => 0, b => 1), 14), Round(exp(1) - 1, 14), + 'gauss_quad(e^x,n=>15) on [-1,1]'; + + my ($nodes, $weights) = legendreP_nodes_weights(14); + is Round(gauss_quad($e_x, a => 0, b => 1, nodes => $nodes, weights => $weights), 14), Round(exp(1) - 1, 14), + 'gauss_quad(e^x,n=>15) on [-1,1]'; + +}; + +subtest 'Runge Kutta 4th order' => sub { + my $f = sub { + my ($t, $y) = @_; + return $t * $t + $y * $y; + }; + my $rk4 = rungeKutta4( + $f, + initial_t => 0, + initial_y => 1, + dt => 0.2, + num_of_points => 5, + interior_points => 1 + ); + is [ map { $_->[0] } @$rk4 ], [ 0, 0.2, 0.4, 0.6, 0.8, 1.0 ], 'returns correct x values'; + is roundArray([ map { $_->[1] } @$rk4 ]), + roundArray([ 1, 1.25299088, 1.6959198, 2.6421097, 5.7854627, 99.9653469 ]), + '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); + return [ map { defined($_) ? Round($_, $options{digits}) : $_ } @$arr ]; +} + +done_testing;