Parent Directory
| 
Revision Log
dev-1-7-01
1 #!/usr/local/bin/perl 2 #use strict; 3 4 ########### 5 6 =head2 Numerical Methods for WeBWorK 7 8 9 =head3 Plotting a list of points 10 11 =pod 12 13 14 15 Usage: plot_list([x0,y0,x1,y1,...]); 16 plot_list([(x0,y0),(x1,y1),...]); 17 plot_list(\x_y_array); 18 19 plot_list([x0,x1,x2...], [y0,y1,y2,...]); 20 plot_list(\@xarray,\@yarray); 21 22 =cut 23 BEGIN { 24 be_strict(); 25 } 26 27 sub plot_list { 28 my($xref,$yref) = @_; 29 unless( defined($xref) && ref($xref) =~/ARRAY/ ) { 30 die "Error in plot_list:X values must be given as an array reference. 31 Remember to use ~~\@array to reference an array in the PG language."; 32 } 33 if (defined($yref) && ! ( ref($yref) =~/ARRAY/ ) ) { 34 die "Error in plot_list:Y values must be given as an array reference. 35 Remember to use ~~\@array to reference an array in the PG language."; 36 } 37 my (@x_vals, @y_vals); 38 unless( defined($yref) ) { #with only one entry we assume (x0,y0,x1,y1..); 39 if ( @$xref % 2 ==1) { 40 die "ERROR in plot_list -- single array of input has odd number of 41 elements"; 42 } 43 44 my @in = @$xref; 45 while (@in) { 46 push(@x_vals, shift(@in)); 47 push(@y_vals, shift(@in)); 48 } 49 $xref = \@x_vals; 50 $yref = \@y_vals; 51 } 52 53 my $fun =sub { 54 my $x = shift; 55 my $y; 56 my( $x0,$x1,$y0,$y1); 57 my @x_values = @$xref; 58 my @y_values = @$yref; 59 while (( @x_values and $x > $x_values[0]) or 60 ( @x_values > 0 and $x >= $x_values[0] ) ) { 61 $x0 = shift(@x_values); 62 $y0 = shift(@y_values); 63 } 64 # now that we have the left hand of the input 65 #check first that x isn't out of range to the left or right 66 if (@x_values && defined($x0) ) { 67 $x1= shift(@x_values); 68 $y1=shift(@y_values); 69 $y = $y0 + ($y1-$y0)*($x-$x0)/($x1-$x0); 70 } 71 $y; 72 }; 73 $fun; 74 } 75 76 =head4 Horner polynomial 77 78 79 =pod 80 81 82 Useage: $fn = horner([x0,x1,x2],[q0,q1,q2]); 83 Produces the newton polynomial 84 &$fn(x) = q0 + q1*(x-x0) +q2*(x-x1)*(x-x0); 85 86 =cut 87 88 sub horner { 89 my ($xref,$qref) = @_; # get the coefficients 90 my $fn = sub { 91 my $x = shift; 92 my @xvals = @$xref; 93 my @qvals = @$qref; 94 my $y = pop(@qvals); 95 pop(@xvals); 96 while (@qvals) { 97 $y = $y * ($x - pop(@xvals) ) + pop(@qvals); 98 } 99 $y; 100 }; 101 $fn; 102 } 103 104 105 =head4 hermite 106 107 108 =pod 109 110 Useage: hermit([x0,x1...],[y0,y1...],[yp0,yp1,...]); 111 Produces a reference to polynomial function 112 with the specified values and first derivatives 113 at (x0,x1,...). 114 115 116 =cut 117 118 119 sub hermite { 120 my($x_ref,$y_ref,$yp_ref) = @_; 121 my (@zvals,@qvals); 122 my $n = $#{$x_ref}; 123 my $i =0; 124 foreach $i (0..$n ) { 125 $zvals[2*$i] = $$x_ref[$i]; 126 $zvals[2*$i+1] = $$x_ref[$i]; 127 $qvals[2*$i][0] = $$y_ref[$i]; 128 $qvals[2*$i+1][0] = $$y_ref[$i]; 129 $qvals[2*$i+1][1] = $$yp_ref[$i]; 130 $qvals[2*$i][1] = ( $qvals[2*$i][0] - $qvals[2*$i-1][0] ) 131 / ( $zvals[2*$i]- $zvals[2*$i-1] ) unless $i ==0; 132 133 } 134 my $j; 135 foreach $i ( 2..(2*$n+1) ) { 136 foreach $j (2 .. $i) { 137 $qvals[$i][$j] = ($qvals[$i][$j-1] - $qvals[$i-1][$j-1]) / 138 ($zvals[$i] - $zvals[$i-$j]); 139 } 140 } 141 my @output; 142 foreach $i (0..2*$n+1) { 143 push(@output,$qvals[$i][$i]); 144 } 145 horner(\@zvals,\@output); 146 } 147 148 149 =head4 Hermite splines 150 151 152 =pod 153 Useage: hermit_spline([x0,x1...],[y0,y1...],[yp0,yp1,...]); 154 Produces a reference to a piecewise cubic hermit spline 155 with the specified values and first derivatives 156 at (x0,x1,...). 157 158 =cut 159 160 161 sub hermite_spline { 162 my ($xref, $yref, $ypref) = @_; 163 my @xvals = @$xref; 164 my @yvals = @$yref; 165 my @ypvals = @$ypref; 166 my $x0 = shift @xvals; 167 my $y0 = shift @yvals; 168 my $yp0 = shift @ypvals; 169 my ($x1,$y1,$yp1); 170 my @polys; #calculate a hermite polynomial evaluator for each region 171 while (@xvals) { 172 $x1 = shift @xvals; 173 $y1 = shift @yvals; 174 $yp1 = shift @ypvals; 175 push @polys, hermite([$x0,$x1],[$y0,$y1],[$yp0,$yp1]); 176 $x0 = $x1; 177 $y0 = $y1; 178 $yp0 = $yp1; 179 } 180 181 182 my $hermite_spline = sub { 183 my $x = shift; 184 my $y; 185 my $fun; 186 my @xvals = @$xref; 187 my @fns = @polys; 188 return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint 189 190 while (@xvals && $x > $xvals[0]) { # find the function for this range of x 191 shift(@xvals); 192 $fun = shift(@fns); 193 } 194 195 # now that we have the left hand of the input 196 #check first that x isn't out of range to the left or right 197 if (@xvals && defined($fun) ) { 198 $y =&$fun($x); 199 } 200 $y; 201 }; 202 $hermite_spline; 203 } 204 205 =head4 Cubic spline approximation 206 207 208 209 Useage: 210 $fun_ref = cubic_spline(~~@x_values, ~~@y_values); 211 212 Where the x and y value arrays come from the function to be approximated. 213 214 The function reference will take a single value x and produce value y. 215 216 $y = &$fun_ref($x); 217 218 219 =cut 220 221 sub cubic_spline { 222 my($t_ref, $y_ref) = @_; 223 my @spline_coeff = create_cubic_spline($t_ref, $y_ref); 224 sub { 225 my $x = shift; 226 eval_cubic_spline($x,@spline_coeff); 227 } 228 } 229 sub eval_cubic_spline { 230 my ($x, $t_ref,$a_ref,$b_ref,$c_ref,$d_ref ) = @_; 231 # The knot points given by $t_ref should be in order. 232 my $i=0; 233 my $out =0; 234 while (defined($t_ref->[$i+1] ) && $x > $t_ref->[$i+1] ) { 235 236 $i++; 237 } 238 unless (defined($t_ref->[$i]) && ( $t_ref->[$i] <= $x ) && ($x <= $t_ref->[$i+1] ) ) { 239 $out = undef; 240 # input value is not in the range defined by the function. 241 } else { 242 $out = ( $t_ref->[$i+1] - $x )* ( ($d_ref->[$i]) +($a_ref->[$i])*( $t_ref->[$i+1] - $x )**2 ) 243 + 244 ( $x - $t_ref->[$i] ) * ( ($b_ref->[$i])*( $x - $t_ref->[$i] )**2 + ($c_ref->[$i]) ) 245 246 } 247 $out; 248 } 249 250 sub create_cubic_spline { 251 my($t_ref, $y_ref) = @_; 252 # The knot points are given by $t_ref (in order) and the function values by $y_ref 253 my $num =$#{$t_ref}; # number of knots 254 my $i; 255 my (@h,@b,@u,@v,@z); 256 foreach $i (0..$num-1) { 257 $h[$i] = $t_ref->[$i+1] - $t_ref->[$i]; 258 $b[$i] = 6*( $y_ref->[$i+1] - $y_ref->[$i] )/$h[$i]; 259 } 260 $u[1] = 2*($h[0] +$h[1] ); 261 $v[1] = $b[1] - $b[0]; 262 foreach $i (2..$num-1) { 263 $u[$i] = 2*( $h[$i] + $h[$i-1] ) - ($h[$i-1])**2/$u[$i-1]; 264 $v[$i] = $b[$i] - $b[$i-1] - $h[$i-1]*$v[$i-1]/$u[$i-1] 265 } 266 $z[$num] = 0; 267 for ($i = $num-1; $i>0; $i--) { 268 $z[$i] = ( $v[$i] - $h[$i]*$z[$i+1] )/$u[$i]; 269 } 270 $z[0] = 0; 271 my ($a_ref, $b_ref, $c_ref, $d_ref); 272 foreach $i (0..$num-1) { 273 $a_ref->[$i] = $z[$i]/(6*$h[$i]); 274 $b_ref->[$i] = $z[$i+1]/(6*$h[$i]); 275 $c_ref->[$i] = ( ($y_ref->[$i+1])/$h[$i] - $z[$i+1]*$h[$i]/6 ); 276 $d_ref->[$i] = ( ($y_ref->[$i])/$h[$i] - $z[$i]*$h[$i]/6 ); 277 } 278 ($t_ref, $a_ref, $b_ref, $c_ref, $d_ref); 279 } 280 281 282 sub javaScript_cubic_spline { 283 my $x_ref = shift; 284 my $y_ref = shift; 285 my %options = @_; 286 $options{name} = 'func' unless defined($options{name} and $options{name} ); 287 $options{llimit} = 'func' unless defined($options{llimit} and $options{llimit}); 288 $options{rlimit} = 'func' unless defined($options{rlimit} and $options{rlimit}); 289 290 my ($t_ref, $a_ref, $b_ref, $c_ref, $d_ref) = create_cubic_spline ($x_ref, $y_ref); 291 292 293 my $str_t_array = "t = new Array(" . join(",",@$t_ref) . ");"; 294 my $str_a_array = "a = new Array(" . join(",",@$a_ref) . ");"; 295 my $str_b_array = "b = new Array(" . join(",",@$b_ref) . ");"; 296 my $str_c_array = "c = new Array(" . join(",",@$c_ref) . ");"; 297 my $str_d_array = "d = new Array(" . join(",",@$d_ref) . ");"; 298 299 my $output_str = <<END_OF_JAVA_TEXT; 300 <SCRIPT LANGUAGE="JavaScript"> 301 <!-- Begin 302 303 304 305 function $options{name}(x) { 306 $str_t_array 307 $str_a_array 308 $str_b_array 309 $str_c_array 310 $str_d_array 311 312 // Evaluate a cubic spline defined by the vectors above 313 i = 0; 314 while (x > t[i+1] ) { 315 i++ 316 } 317 318 if ( t[i] <= x && x <= t[i+1] && $options{llimit} <= x && x <= $options{rlimit} ) { 319 return ( ( t[i+1] - x )*( d[i] +a[i]*( t[i+1] - x )*( t[i+1] - x ) ) 320 + ( x - t[i] )*( b[i]*( x - t[i])*( x - t[i] ) +c[i] ) 321 ); 322 323 } else { 324 return( "undefined" ) ; 325 } 326 327 } 328 // End 329 --> 330 </SCRIPT> 331 <NOSCRIPT> 332 This problem requires a browser capable of processing javaScript 333 </NOSCRIPT> 334 END_OF_JAVA_TEXT 335 336 $output_str; 337 } 338 339 340 341 342 =head4 Integration by trapezoid rule 343 344 =pod 345 Useage: trapezoid(function_reference, start, end, [steps:30]); 346 347 348 =cut 349 350 sub trapezoid { 351 my $fn_ref = shift; 352 my $x0 = shift; 353 my $x1 = shift; 354 my $steps = shift; 355 $steps = 30 unless defined $steps; 356 my $delta =($x1-$x0)/$steps; 357 my $i; 358 my $sum=0; 359 foreach $i (1..($steps-1) ) { 360 $sum =$sum + &$fn_ref( $x0 + $i*$delta ); 361 } 362 $sum = $sum + 0.5*( &$fn_ref($x0) + &$fn_ref($x1) ); 363 $sum*$delta; 364 } 365 366 =head4 Romberg method of integration 367 368 =pod 369 370 Useage: romberg(function_reference, x0, x1, level); 371 372 =cut 373 374 sub romberg_iter { 375 my $fn_ref = shift; 376 my $x0 = shift; 377 my $x1 = shift; 378 my $j = shift; 379 my $k = shift; 380 my $out; 381 if ($k == 1 ) { 382 $out = trapezoid($fn_ref, $x0,$x1,2**($j-1) ); 383 } else { 384 385 $out = ( 4**($k-1) * romberg_iter($fn_ref, $x0,$x1,$j,$k-1) - 386 romberg_iter($fn_ref, $x0,$x1,$j-1,$k-1) ) / ( 4**($k-1) -1) ; 387 } 388 $out; 389 } 390 391 392 sub romberg { 393 my $fn_ref = shift; 394 my $x0 = shift; 395 my $x1 = shift; 396 my $level = shift; 397 $level = 6 unless defined $level; 398 romberg_iter($fn_ref, $x0,$x1, $level,$level); 399 } 400 401 ######################################################### 402 403 1; 404
| aubreyja at gmail dot com | ViewVC Help |
| Powered by ViewVC 1.0.9 |