[system] / trunk / webwork / system / courseScripts / PGnumericalmacros.pl Repository:
ViewVC logotype

View of /trunk/webwork/system/courseScripts/PGnumericalmacros.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 23 - (download) (as text) (annotate)
Tue Jun 19 15:37:34 2001 UTC (18 years, 5 months ago) by gage
File size: 11392 byte(s)
Modified the  documentation for the trapezoid function.

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9