[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 6538 - (download) (as text) (annotate)
Sat Nov 20 23:46:23 2010 UTC (9 years, 1 month ago) by gage
File size: 17324 byte(s)
added routines provided by Steven McKay



    1 
    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 Left Hand Sum
  397 
  398 =pod
  399 
  400   Usage:  lefthandsum(function_reference, start, end, steps=>30 );
  401 
  402 Implements the Left Hand sum using 30 intervals between 'start' and 'end'.
  403 The first three arguments are required.  The final argument (number of steps) is optional and defaults to 30.
  404 
  405 =cut
  406 
  407 sub lefthandsum {
  408   my $fn_ref = shift;
  409         my $x0 = shift;
  410         my $x1 = shift;
  411         my %options = @_;
  412         assign_option_aliases(\%options,
  413                               intervals   =>  'steps',
  414   );
  415   set_default_options(\%options,
  416             steps     =>  30,
  417         );
  418   my $steps = $options{steps};
  419   my $delta = ($x1-$x0)/$steps;
  420         my $i;
  421   my $sum=0;
  422   foreach $i (0..($steps-1) ) {
  423     $sum =$sum + &$fn_ref( $x0 + ($i)*$delta );
  424   }
  425   $sum*$delta;
  426 }
  427 
  428 =head3 Integration by Right Hand Sum
  429 
  430 =pod
  431 
  432   Usage:  righthandsum(function_reference, start, end, steps=>30 );
  433 
  434 Implements the right hand sum using 30 intervals between 'start' and 'end'.
  435 The first three arguments are required.  The final argument (number of steps) is optional and defaults to 30.
  436 
  437 =cut
  438 
  439 sub righthandsum {
  440   my $fn_ref = shift;
  441         my $x0 = shift;
  442         my $x1 = shift;
  443         my %options = @_;
  444         assign_option_aliases(\%options,
  445                               intervals   =>  'steps',
  446   );
  447   set_default_options(\%options,
  448             steps     =>  30,
  449         );
  450   my $steps = $options{steps};
  451   my $delta = ($x1-$x0)/$steps;
  452         my $i;
  453   my $sum=0;
  454   foreach $i (1..($steps) ) {
  455     $sum =$sum + &$fn_ref( $x0 + ($i)*$delta );
  456   }
  457   $sum*$delta;
  458 }
  459 
  460 
  461 =head3 Integration by Midpoint rule
  462 
  463 =pod
  464 
  465   Usage:  midpoint(function_reference, start, end, steps=>30 );
  466 
  467 Implements the Midpoint rule using 30 intervals between 'start' and 'end'.
  468 The first three arguments are required.  The final argument (number of steps) is optional and defaults to 30.
  469 
  470 =cut
  471 
  472 sub midpoint {
  473   my $fn_ref = shift;
  474         my $x0 = shift;
  475         my $x1 = shift;
  476         my %options = @_;
  477         assign_option_aliases(\%options,
  478                               intervals   =>  'steps',
  479   );
  480   set_default_options(\%options,
  481             steps     =>  30,
  482         );
  483   my $steps = $options{steps};
  484   my $delta = ($x1-$x0)/$steps;
  485         my $i;
  486   my $sum=0;
  487   foreach $i (0..($steps-1) ) {
  488     $sum =$sum + &$fn_ref( $x0 + ($i+1/2)*$delta );
  489   }
  490   $sum*$delta;
  491 }
  492 
  493 =head3 Integration by Simpson's rule
  494 
  495 =pod
  496 
  497   Usage:  simpson(function_reference, start, end, steps=>30 );
  498 
  499 Implements Simpson's rule using 30 intervals between 'start' and 'end'.
  500 The first three arguments are required.  The final argument (number of steps) is optional and defaults to 30,
  501 but must be even.
  502 
  503 =cut
  504 
  505 
  506 sub simpson {
  507   my $fn_ref = shift;
  508         my $x0 = shift;
  509         my $x1 = shift;
  510         my %options = @_;
  511         assign_option_aliases(\%options,
  512                               intervals   =>  'steps',
  513   );
  514   set_default_options(\%options,
  515             steps     =>  30,
  516         );
  517   my $steps = $options{steps};
  518       unless( $steps % 2 == 0  ) {
  519       die "Error: Simpson's rule requires an even number of steps.";
  520     }
  521 
  522   my $delta = ($x1-$x0)/$steps;
  523         my $i;
  524   my $sum=0;
  525   for ($i=1;$i<$steps;$i=$i+2) { # look this up - loop by two.
  526     $sum =$sum + 4*&$fn_ref( $x0 + $i*$delta );
  527   }
  528   for ($i=2;$i<$steps-1;$i=$i+2) { # ditto
  529     $sum = $sum + 2*&$fn_ref( $x0 + $i*$delta);
  530   }
  531   $sum = $sum + &$fn_ref($x0) + &$fn_ref($x1);
  532   $sum*$delta/3;
  533 }
  534 
  535 
  536 =head3 Integration by trapezoid rule
  537 
  538 =pod
  539 
  540   Usage:  trapezoid(function_reference, start, end, steps=>30 );
  541 
  542 Implements the trapezoid rule using 30 intervals between 'start' and 'end'.
  543 The first three arguments are required.  The final argument (number of steps)
  544 is optional and defaults to 30.
  545 
  546 =cut
  547 
  548 
  549 sub trapezoid {
  550   my $fn_ref = shift;
  551   my $x0 = shift;
  552   my $x1 = shift;
  553   my %options = @_;
  554   assign_option_aliases(\%options,
  555                          intervals  =>  'steps',
  556   );
  557   set_default_options(\%options,
  558             steps     =>  30,
  559   );
  560   my $steps = $options{steps};
  561   my $delta =($x1-$x0)/$steps;
  562         my $i;
  563   my $sum=0;
  564   foreach $i (1..($steps-1) ) {
  565     $sum =$sum + &$fn_ref( $x0 + $i*$delta );
  566   }
  567   $sum = $sum + 0.5*( &$fn_ref($x0) + &$fn_ref($x1) );
  568   $sum*$delta;
  569 }
  570 
  571 =head3  Romberg method of integration
  572 
  573 =pod
  574 
  575         Usage:  romberg(function_reference, x0, x1, level);
  576 
  577 Implements the Romberg integration routine through 'level' recursive steps.  Level defaults to 6.
  578 
  579 =cut
  580 
  581 
  582 sub romberg {
  583   my $fn_ref = shift;
  584   my $x0 = shift;
  585   my $x1 = shift;
  586   my %options = @_;
  587   set_default_options(\%options,
  588             level     =>  6,
  589   );
  590   my $level = $options{level};
  591   romberg_iter($fn_ref, $x0,$x1, $level,$level);
  592 }
  593 
  594 sub romberg_iter {
  595   my $fn_ref = shift;
  596   my $x0 = shift;
  597   my $x1 = shift;
  598   my $j  = shift;
  599   my $k  = shift;
  600   my $out;
  601   if ($k == 1 ) {
  602     $out = trapezoid($fn_ref, $x0,$x1,steps => 2**($j-1) );
  603   } else {
  604 
  605     $out = (  4**($k-1) * romberg_iter($fn_ref, $x0,$x1,$j,$k-1) -
  606             romberg_iter($fn_ref, $x0,$x1,$j-1,$k-1) ) / ( 4**($k-1) -1) ;
  607   }
  608   $out;
  609 }
  610 
  611 =head3 Inverse Romberg
  612 
  613 =pod
  614 
  615         Usage: inv_romberg(function_reference, a, value);
  616 
  617 Finds b such that the integral of the function from a to b is equal to value.
  618 Assumes that the function is continuous and doesn't take on the zero value.
  619 Uses Newton's method of approximating roots of equations, and Romberg to evaluate definite integrals.
  620 
  621 =cut
  622 
  623 sub inv_romberg {
  624         my $fn_ref = shift;
  625         my $a = shift;
  626         my $value = shift;
  627         my $b = $a;
  628         my $delta = 1;
  629         my $i=0;
  630         my $funct;
  631         my $deriv;
  632         while (abs($delta) > 0.000001 && $i < 5000) {
  633                 $funct = romberg($fn_ref,$a,$b)-$value;
  634                 $deriv = &$fn_ref ( $b );
  635           if ($deriv == 0) {
  636       warn 'Values of the function are too close to 0.';
  637                         return;
  638                 }
  639                 $delta = $funct/$deriv;
  640                 $b = $b - $delta;
  641                 $i++;
  642         }
  643         if ($i == 5000) {
  644                 warn 'Newtons method does not converge.';
  645                 return;
  646   }
  647         $b;
  648 }
  649 
  650 #########################################################
  651 =pod
  652 
  653     rungeKutta4
  654 
  655     Finds integral curve of a vector field using the 4th order Runge Kutta method.
  656 
  657   Useage:  rungeKutta4( &vectorField(t,x),%options);
  658 
  659     Returns:  \@array of points [t,y})
  660 
  661     Default %options:
  662           'initial_t'         =>  1,
  663           'initial_y'         =>  1,
  664           'dt'            =>  .01,
  665           'num_of_points'       =>  10,     #number of reported points
  666           'interior_points'     =>  5,      # number of 'interior' steps between reported points
  667           'debug'
  668 
  669 =cut
  670 
  671 sub rungeKutta4 {
  672   my $rf_fun = shift;
  673   my %options = @_;
  674   set_default_options(  \%options,
  675           'initial_t'         =>  1,
  676           'initial_y'         =>  1,
  677           'dt'            =>  .01,
  678           'num_of_points'       =>  10,     #number of reported points
  679           'interior_points'     =>  5,      # number of 'interior' steps between reported points
  680           'debug'           =>  1,      # remind programmers to always pass the debug parameter
  681           );
  682   my $t = $options{initial_t};
  683   my $y = $options{initial_y};
  684 
  685   my $num = $options{'num_of_points'};  # number of points
  686   my $num2 = $options{'interior_points'};  # number of steps between points.
  687   my $dt   = $options{'dt'};
  688   my $errors = undef;
  689   my $rf_rhs = sub {  my @in = @_;
  690         my ( $out, $err) = &$rf_fun(@in);
  691         $errors .= " $err at ( ".join(" , ", @in) . " )<br>\n" if defined($err);
  692         $out = 'NaN' if defined($err) and not is_a_number($out);
  693         $out;
  694           };
  695 
  696   my @output = ([$t, $y]);
  697   my ($i, $j, $K1,$K2,$K3,$K4);
  698 
  699   for ($j=0; $j<$num; $j++) {
  700       for ($i=0; $i<$num2; $i++) {
  701     $K1 = $dt*&$rf_rhs($t, $y);
  702     $K2 = $dt*&$rf_rhs($t+$dt/2,$y+$K1/2);
  703     $K3 = $dt*&$rf_rhs($t+$dt/2, $y+$K2/2);
  704     $K4 = $dt*&$rf_rhs($t+$dt, $y+$K3);
  705     $y = $y + ($K1 + 2*$K2 + 2*$K3 + $K4)/6;
  706     $t = $t + $dt;
  707       }
  708       push(@output, [$t, $y]);
  709   }
  710   if (defined $errors) {
  711     return $errors;
  712   } else {
  713     return \@output;
  714   }
  715 }
  716 
  717 
  718 
  719 1;
  720 

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9