[system] / trunk / pg / macros / PGnumericalmacros.pl Repository:
ViewVC logotype

View of /trunk/pg/macros/PGnumericalmacros.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1050 - (download) (as text) (annotate)
Fri Jun 6 21:39:42 2003 UTC (16 years, 8 months ago) by sh002i
File size: 12263 byte(s)
moved PG modules and macro files from webwork-modperl to pg
-sam

    1 #!/usr/local/bin/webwork-perl
    2 #use strict;
    3 
    4 ###########
    5 #use Carp;
    6 
    7 =head1 NAME
    8 
    9   Numerical methods for the PG language
   10 
   11 =head1 SYNPOSIS
   12 
   13 
   14 
   15 =head1 DESCRIPTION
   16 
   17 =cut
   18 
   19 =head2 Interpolation  methods
   20 
   21 =head3 Plotting a list of points (piecewise linear interpolation)
   22 
   23 
   24     Usage:  plot_list([x0,y0,x1,y1,...]);
   25         plot_list([(x0,y0),(x1,y1),...]);
   26         plot_list(\x_y_array);
   27 
   28             plot_list([x0,x1,x2...], [y0,y1,y2,...]);
   29             plot_list(\@xarray,\@yarray);
   30 
   31 
   32 =cut
   33 
   34 BEGIN {
   35   be_strict();
   36 }
   37 sub _PGnumericalmacros_init {
   38 }
   39 
   40 sub plot_list {
   41     my($xref,$yref) = @_;
   42     unless( defined($xref) && ref($xref) =~/ARRAY/  ) {
   43       die "Error in plot_list:X values must be given as an array reference.
   44          Remember to use ~~\@array to reference an array in the PG language.";
   45     }
   46     if (defined($yref) && ! ( ref($yref) =~/ARRAY/ ) ) {
   47       die "Error in plot_list:Y values must be given as an array reference.
   48          Remember to use ~~\@array to reference an array in the PG language.";
   49     }
   50     my (@x_vals, @y_vals);
   51     unless( defined($yref) ) { #with only one entry we assume (x0,y0,x1,y1..);
   52         if ( @$xref % 2 ==1) {
   53           die "ERROR in plot_list -- single array of input has odd number of
   54           elements";
   55         }
   56 
   57      my  @in = @$xref;
   58      while (@in) {
   59         push(@x_vals, shift(@in));
   60         push(@y_vals, shift(@in));
   61       }
   62       $xref = \@x_vals;
   63       $yref = \@y_vals;
   64   }
   65 
   66   my $fun =sub {
   67     my $x = shift;
   68     my $y;
   69     my( $x0,$x1,$y0,$y1);
   70     my @x_values = @$xref;
   71     my @y_values = @$yref;
   72     while (( @x_values and $x > $x_values[0]) or
   73            ( @x_values > 0 and $x >= $x_values[0] ) ) {
   74       $x0 = shift(@x_values);
   75       $y0 = shift(@y_values);
   76     }
   77     # now that we have the left hand of the input
   78     #check first that x isn't out of range to the left or right
   79     if (@x_values  && defined($x0) )  {
   80       $x1= shift(@x_values);
   81       $y1=shift(@y_values);
   82       $y = $y0 + ($y1-$y0)*($x-$x0)/($x1-$x0);
   83     }
   84     $y;
   85   };
   86   $fun;
   87 }
   88 
   89 =head3 Horner polynomial/ Newton polynomial
   90 
   91 
   92 
   93   Usege:  $fn = horner([x0,x1,x2],[q0,q1,q2]);
   94     Produces the newton polynomial
   95     &$fn(x) = q0 + q1*(x-x0) +q2*(x-x1)*(x-x0);
   96 
   97 Generates a subroutine which evaluates a polynomial passing through the points C<(x0,q0), (x1,q1),
   98 ... > using Horner's method.
   99 
  100 =cut
  101 
  102 sub horner {
  103   my ($xref,$qref) = @_;  # get the coefficients
  104   my $fn = sub {
  105     my $x = shift;
  106     my @xvals = @$xref;
  107     my @qvals = @$qref;
  108     my $y = pop(@qvals);
  109             pop(@xvals);
  110       while (@qvals) {
  111         $y = $y * ($x - pop(@xvals) )  + pop(@qvals);
  112       }
  113       $y;
  114    };
  115    $fn;
  116 }
  117 
  118 
  119 =head3 Hermite polynomials
  120 
  121 
  122 =pod
  123 
  124   Usage:  $poly = hermit([x0,x1...],[y0,y1...],[yp0,yp1,...]);
  125     Produces a reference to polynomial function
  126     with the specified values and first derivatives
  127     at (x0,x1,...).
  128     &$poly(34) gives a number
  129 
  130 Generates a subroutine which evaluates a polynomial passing through the specified points
  131 with the specified derivatives: (x0,y0,yp0) ...
  132 The polynomial will be of high degree and may wobble unexpectedly.  Use the Hermite splines
  133 described below and in Hermite.pm for  most graphing purposes.
  134 
  135 =cut
  136 
  137 
  138 sub hermite {
  139   my($x_ref,$y_ref,$yp_ref) = @_;
  140   my (@zvals,@qvals);
  141   my $n = $#{$x_ref};
  142   my $i =0;
  143   foreach $i (0..$n ) {
  144     $zvals[2*$i] = $$x_ref[$i];
  145     $zvals[2*$i+1] = $$x_ref[$i];
  146     $qvals[2*$i][0] = $$y_ref[$i];
  147     $qvals[2*$i+1][0] = $$y_ref[$i];
  148     $qvals[2*$i+1][1] = $$yp_ref[$i];
  149     $qvals[2*$i][1]  = ( $qvals[2*$i][0] - $qvals[2*$i-1][0] )
  150                       / ( $zvals[2*$i]- $zvals[2*$i-1] ) unless $i ==0;
  151 
  152   }
  153   my $j;
  154   foreach $i ( 2..(2*$n+1) )  {
  155     foreach $j (2 .. $i) {
  156       $qvals[$i][$j] = ($qvals[$i][$j-1] - $qvals[$i-1][$j-1]) /
  157                         ($zvals[$i] - $zvals[$i-$j]);
  158     }
  159   }
  160   my @output;
  161   foreach $i (0..2*$n+1) {
  162     push(@output,$qvals[$i][$i]);
  163   }
  164   horner(\@zvals,\@output);
  165 }
  166 
  167 
  168 =head3 Hermite splines
  169 
  170 
  171   Usage:  $spline = hermit_spline([x0,x1...],[y0,y1...],[yp0,yp1,...]);
  172     Produces a reference to a piecewise cubic hermit spline
  173     with the specified values and first derivatives
  174     at (x0,x1,...).
  175 
  176     &$spline(45) evaluates to a number.
  177 
  178 Generates a subroutine which evaluates a piecewise cubic polynomial
  179 passing through the specified points
  180 with the specified derivatives: (x0,y0,yp0) ...
  181 
  182 An object oriented version of this is defined in Hermite.pm
  183 
  184 =cut
  185 
  186 
  187 sub hermite_spline {
  188   my ($xref, $yref, $ypref) = @_;
  189   my @xvals  = @$xref;
  190   my @yvals  = @$yref;
  191   my @ypvals = @$ypref;
  192   my $x0 = shift @xvals;
  193   my $y0 = shift @yvals;
  194   my $yp0 = shift @ypvals;
  195   my ($x1,$y1,$yp1);
  196   my @polys;  #calculate a hermite polynomial evaluator for each region
  197   while (@xvals) {
  198     $x1 = shift @xvals;
  199     $y1 = shift @yvals;
  200     $yp1 = shift @ypvals;
  201     push @polys, hermite([$x0,$x1],[$y0,$y1],[$yp0,$yp1]);
  202     $x0  = $x1;
  203     $y0  = $y1;
  204     $yp0 = $yp1;
  205   }
  206 
  207 
  208   my $hermite_spline_sub = sub {
  209     my $x = shift;
  210     my $y;
  211     my $fun;
  212     my @xvals = @$xref;
  213     my @fns = @polys;
  214     return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
  215 
  216     while (@xvals && $x > $xvals[0]) {  # find the function for this range of x
  217       shift(@xvals);
  218       $fun = shift(@fns);
  219     }
  220 
  221     # now that we have the left hand of the input
  222     #check first that x isn't out of range to the left or right
  223     if (@xvals  && defined($fun) )  {
  224       $y =&$fun($x);
  225     }
  226     $y;
  227   };
  228   $hermite_spline_sub;
  229 }
  230 
  231 =head3 Cubic spline approximation
  232 
  233 
  234 
  235     Usage:
  236          $fun_ref = cubic_spline(~~@x_values, ~~@y_values);
  237 
  238 Where the x and y value arrays come from the function to be approximated.
  239 The function reference will take a single value x and produce value y.
  240 
  241       $y = &$fun_ref($x);
  242 
  243 You can also generate javaScript which defines a cubic spline:
  244 
  245         $function_string = javaScript_cubic_spline(~~@_x_values, ~~@y_values,
  246                   name =>  'myfunction1',
  247                   llimit => -3,
  248                   rlimit => 3,
  249                   );
  250 
  251 The string contains
  252 
  253   <SCRIPT LANGUAGE="JavaScript">
  254   <!-- Begin
  255   function myfunction1(x) {
  256   ...etc...
  257   }
  258   </SCRIPT>
  259 
  260 and can be placed in the header of the HTML output using
  261 
  262   HEADER_TEXT($function_string);
  263 
  264 =cut
  265 
  266 sub cubic_spline {
  267   my($t_ref, $y_ref) = @_;
  268   my @spline_coeff = create_cubic_spline($t_ref, $y_ref);
  269   sub {
  270     my $x = shift;
  271     eval_cubic_spline($x,@spline_coeff);
  272   }
  273 }
  274 sub eval_cubic_spline {
  275   my ($x, $t_ref,$a_ref,$b_ref,$c_ref,$d_ref ) = @_;
  276 # The knot points given by $t_ref should be in order.
  277   my $i=0;
  278   my $out =0;
  279   while (defined($t_ref->[$i+1] ) && $x  >  $t_ref->[$i+1] ) {
  280 
  281     $i++;
  282   }
  283   unless (defined($t_ref->[$i]) && ( $t_ref->[$i] <= $x ) && ($x <= $t_ref->[$i+1] ) ) {
  284     $out = undef;
  285     # input value is not in the range defined by the function.
  286   } else {
  287     $out = ( $t_ref->[$i+1]   - $x )* ( ($d_ref->[$i]) +($a_ref->[$i])*( $t_ref->[$i+1] - $x )**2 )
  288            +
  289            ( $x  -  $t_ref->[$i] ) * ( ($b_ref->[$i])*( $x  -  $t_ref->[$i] )**2  + ($c_ref->[$i]) )
  290 
  291   }
  292   $out;
  293 }
  294 
  295 ##########################################################################
  296 # Cubic spline algorithm adapted from p319 of Kincaid and Cheney's Numerical Analysis
  297 ##########################################################################
  298 
  299 sub create_cubic_spline {
  300   my($t_ref, $y_ref) = @_;
  301 # The knot points are given by $t_ref (in order) and the function values by $y_ref
  302   my $num =$#{$t_ref};  # number of knots
  303   my $i;
  304   my (@h,@b,@u,@v,@z);
  305   foreach $i (0..$num-1) {
  306     $h[$i] = $t_ref->[$i+1] - $t_ref->[$i];
  307     $b[$i] = 6*( $y_ref->[$i+1] - $y_ref->[$i] )/$h[$i];
  308   }
  309   $u[1] = 2*($h[0] +$h[1] );
  310   $v[1] = $b[1] - $b[0];
  311   foreach $i (2..$num-1)    {
  312     $u[$i] = 2*( $h[$i] + $h[$i-1] ) - ($h[$i-1])**2/$u[$i-1];
  313     $v[$i] = $b[$i] - $b[$i-1] - $h[$i-1]*$v[$i-1]/$u[$i-1]
  314   }
  315   $z[$num] = 0;
  316   for ($i = $num-1; $i>0; $i--)  {
  317     $z[$i] = ( $v[$i] - $h[$i]*$z[$i+1] )/$u[$i];
  318   }
  319   $z[0] = 0;
  320   my ($a_ref, $b_ref, $c_ref, $d_ref);
  321   foreach $i (0..$num-1) {
  322     $a_ref->[$i] = $z[$i]/(6*$h[$i]);
  323     $b_ref->[$i] = $z[$i+1]/(6*$h[$i]);
  324     $c_ref->[$i] = ( ($y_ref->[$i+1])/$h[$i] - $z[$i+1]*$h[$i]/6 );
  325     $d_ref->[$i] = ( ($y_ref->[$i])/$h[$i] - $z[$i]*$h[$i]/6  );
  326   }
  327   ($t_ref, $a_ref, $b_ref, $c_ref, $d_ref);
  328 }
  329 
  330 
  331 sub javaScript_cubic_spline {
  332   my $x_ref = shift;
  333   my $y_ref = shift;
  334   my %options = @_;
  335   assign_option_aliases(\%options,
  336 
  337   );
  338   set_default_options(\%options,
  339             name      =>  'func',
  340             llimit      =>  $x_ref->[0],
  341             rlimit      =>  $x_ref->[$#$x_ref],
  342   );
  343 
  344   my ($t_ref, $a_ref, $b_ref, $c_ref, $d_ref) = create_cubic_spline ($x_ref, $y_ref);
  345 
  346 
  347 my $str_t_array = "t = new Array(" . join(",",@$t_ref) . ");";
  348 my $str_a_array = "a = new Array(" . join(",",@$a_ref) . ");";
  349 my $str_b_array = "b = new Array(" . join(",",@$b_ref) . ");";
  350 my $str_c_array = "c = new Array(" . join(",",@$c_ref) . ");";
  351 my $str_d_array = "d = new Array(" . join(",",@$d_ref) . ");";
  352 
  353 my $output_str = <<END_OF_JAVA_TEXT;
  354 <SCRIPT LANGUAGE="JavaScript">
  355 <!-- Begin
  356 
  357 
  358 
  359 function $options{name}(x) {
  360   $str_t_array
  361   $str_a_array
  362   $str_b_array
  363   $str_c_array
  364   $str_d_array
  365 
  366   // Evaluate a cubic spline defined by the vectors above
  367   i = 0;
  368   while (x > t[i+1] ) {
  369     i++
  370   }
  371 
  372   if ( t[i] <= x && x <= t[i+1]  && $options{llimit} <= x && x <= $options{rlimit} ) {
  373     return (   ( t[i+1] - x )*( d[i] +a[i]*( t[i+1] - x )*( t[i+1] - x ) )
  374              + ( x -   t[i] )*( b[i]*( x - t[i])*( x - t[i] ) +c[i] )
  375            );
  376 
  377   } else {
  378     return( "undefined" ) ;
  379   }
  380 
  381 }
  382 // End
  383  -->
  384 </SCRIPT>
  385 <NOSCRIPT>
  386 This problem requires a browser capable of processing javaScript
  387 </NOSCRIPT>
  388 END_OF_JAVA_TEXT
  389 
  390 $output_str;
  391 }
  392 
  393 
  394 =head2 Numerical Integration methods
  395 
  396 =head3 Integration by trapezoid rule
  397 
  398 =pod
  399 
  400   Usage:  trapezoid(function_reference, start, end, steps=>30 );
  401 
  402 Implements the trapezoid rule using 30 intervals between 'start' and 'end'.
  403 The first three arguments are required.  The final argument (number of steps)
  404 is optional and defaults to 30.
  405 
  406 =cut
  407 
  408 
  409 sub trapezoid {
  410   my $fn_ref = shift;
  411   my $x0 = shift;
  412   my $x1 = shift;
  413   my %options = @_;
  414   assign_option_aliases(\%options,
  415                          intervals  =>  'steps',
  416   );
  417   set_default_options(\%options,
  418             steps     =>  30,
  419   );
  420   my $steps = $options{steps};
  421   my $delta =($x1-$x0)/$steps;
  422         my $i;
  423   my $sum=0;
  424   foreach $i (1..($steps-1) ) {
  425     $sum =$sum + &$fn_ref( $x0 + $i*$delta );
  426   }
  427   $sum = $sum + 0.5*( &$fn_ref($x0) + &$fn_ref($x1) );
  428   $sum*$delta;
  429 }
  430 
  431 =head3  Romberg method of integration
  432 
  433 =pod
  434 
  435         Usage:  romberg(function_reference, x0, x1, level);
  436 
  437 Implements the Romberg integration routine through 'level' recursive steps.  Level defaults to 6.
  438 
  439 =cut
  440 
  441 
  442 sub romberg {
  443   my $fn_ref = shift;
  444   my $x0 = shift;
  445   my $x1 = shift;
  446   my %options = @_;
  447   set_default_options(\%options,
  448             level     =>  6,
  449   );
  450   my $level = $options{level};
  451   romberg_iter($fn_ref, $x0,$x1, $level,$level);
  452 }
  453 
  454 sub romberg_iter {
  455   my $fn_ref = shift;
  456   my $x0 = shift;
  457   my $x1 = shift;
  458   my $j  = shift;
  459   my $k  = shift;
  460   my $out;
  461   if ($k == 1 ) {
  462     $out = trapezoid($fn_ref, $x0,$x1,steps => 2**($j-1) );
  463   } else {
  464 
  465     $out = (  4**($k-1) * romberg_iter($fn_ref, $x0,$x1,$j,$k-1) -
  466             romberg_iter($fn_ref, $x0,$x1,$j-1,$k-1) ) / ( 4**($k-1) -1) ;
  467   }
  468   $out;
  469 }
  470 
  471 =head3 Inverse Romberg
  472 
  473 =pod
  474 
  475         Usage: inv_romberg(function_reference, a, value);
  476 
  477 Finds b such that the integral of the function from a to b is equal to value.
  478 Assumes that the function is continuous and doesn't take on the zero value.
  479 Uses Newton's method of approximating roots of equations, and Romberg to evaluate definite integrals.
  480 
  481 =cut
  482 
  483 sub inv_romberg {
  484         my $fn_ref = shift;
  485         my $a = shift;
  486         my $value = shift;
  487         my $b = $a;
  488         my $delta = 1;
  489         my $i=0;
  490         my $funct;
  491         my $deriv;
  492         while (abs($delta) > 0.000001 && $i < 5000) {
  493                 $funct = romberg($fn_ref,$a,$b)-$value;
  494                 $deriv = &$fn_ref ( $b );
  495           if ($deriv == 0) {
  496       warn 'Values of the function are too close to 0.';
  497                         return;
  498                 }
  499                 $delta = $funct/$deriv;
  500                 $b = $b - $delta;
  501                 $i++;
  502         }
  503         if ($i == 5000) {
  504                 warn 'Newtons method does not converge.';
  505                 return;
  506   }
  507         $b;
  508 }
  509 
  510 #########################################################
  511 
  512 1;
  513 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9