Parent Directory
| 
Revision Log
Changed input to rungeKutta4 so that it matches the pod description -- i.e. it takes as its first argument a reference to a subroutine defining a vectorField as opposed to an answer hash.
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 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 =pod 512 513 rungeKutta4 514 515 Finds integral curve of a vector field using the 4th order Runge Kutta method. 516 517 Useage: rungeKutta4( &vectorField(t,x),%options); 518 519 Returns: \@array of points [t,y}) 520 521 Default %options: 522 'initial_t' => 1, 523 'initial_y' => 1, 524 'dt' => .01, 525 'num_of_points' => 10, #number of reported points 526 'interior_points' => 5, # number of 'interior' steps between reported points 527 'debug' 528 529 =cut 530 531 sub rungeKutta4 { 532 my $rf_fun = shift; 533 my %options = @_; 534 set_default_options( \%options, 535 'initial_t' => 1, 536 'initial_y' => 1, 537 'dt' => .01, 538 'num_of_points' => 10, #number of reported points 539 'interior_points' => 5, # number of 'interior' steps between reported points 540 'debug' => 1, # remind programmers to always pass the debug parameter 541 ); 542 my $t = $options{initial_t}; 543 my $y = $options{initial_y}; 544 545 my $num = $options{'num_of_points'}; # number of points 546 my $num2 = $options{'interior_points'}; # number of steps between points. 547 my $dt = $options{'dt'}; 548 my $errors = undef; 549 my $rf_rhs = sub { my @in = @_; 550 my ( $out, $err) = &$rf_fun(@in); 551 $errors .= " $err at ( ".join(" , ", @in) . " )<br>\n" if defined($err); 552 $out = 'NaN' if defined($err) and not is_a_number($out); 553 $out; 554 }; 555 556 my @output = ([$t, $y]); 557 my ($i, $j, $K1,$K2,$K3,$K4); 558 559 for ($j=0; $j<$num; $j++) { 560 for ($i=0; $i<$num2; $i++) { 561 $K1 = $dt*&$rf_rhs($t, $y); 562 $K2 = $dt*&$rf_rhs($t+$dt/2,$y+$K1/2); 563 $K3 = $dt*&$rf_rhs($t+$dt/2, $y+$K2/2); 564 $K4 = $dt*&$rf_rhs($t+$dt, $y+$K3); 565 $y = $y + ($K1 + 2*$K2 + 2*$K3 + $K4)/6; 566 $t = $t + $dt; 567 } 568 push(@output, [$t, $y]); 569 } 570 if (defined $errors) { 571 return $errors; 572 } else { 573 return \@output; 574 } 575 } 576 577 578 579 1; 580
| aubreyja at gmail dot com | ViewVC Help |
| Powered by ViewVC 1.0.9 |