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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1080 - (view) (download) (as text)

1 : apizer 1080
2 : sh002i 1050 #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 : apizer 1080
28 : sh002i 1050 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 : apizer 1080 die "ERROR in plot_list -- single array of input has odd number of
54 : sh002i 1050 elements";
55 :     }
56 : apizer 1080
57 : sh002i 1050 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 : apizer 1080 }
65 :    
66 : sh002i 1050 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 : apizer 1080 ... > using Horner's method.
99 : sh002i 1050
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 : apizer 1080 Produces a reference to polynomial function
126 : sh002i 1050 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 : apizer 1080 $qvals[2*$i][1] = ( $qvals[2*$i][0] - $qvals[2*$i-1][0] )
150 : sh002i 1050 / ( $zvals[2*$i]- $zvals[2*$i-1] ) unless $i ==0;
151 : apizer 1080
152 : sh002i 1050 }
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 : apizer 1080 horner(\@zvals,\@output);
165 : sh002i 1050 }
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 : apizer 1080
176 : sh002i 1050 &$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 : apizer 1080
207 :    
208 : sh002i 1050 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 : apizer 1080
216 : sh002i 1050 while (@xvals && $x > $xvals[0]) { # find the function for this range of x
217 :     shift(@xvals);
218 :     $fun = shift(@fns);
219 :     }
220 : apizer 1080
221 : sh002i 1050 # 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 : apizer 1080
238 : sh002i 1050 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 : apizer 1080
243 : sh002i 1050 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 : apizer 1080 and can be placed in the header of the HTML output using
261 : sh002i 1050
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 : apizer 1080 my ($x, $t_ref,$a_ref,$b_ref,$c_ref,$d_ref ) = @_;
276 : sh002i 1050 # 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 : apizer 1080
281 : sh002i 1050 $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 : apizer 1080 $out = ( $t_ref->[$i+1] - $x )* ( ($d_ref->[$i]) +($a_ref->[$i])*( $t_ref->[$i+1] - $x )**2 )
288 : sh002i 1050 +
289 : apizer 1080 ( $x - $t_ref->[$i] ) * ( ($b_ref->[$i])*( $x - $t_ref->[$i] )**2 + ($c_ref->[$i]) )
290 :    
291 :     }
292 : sh002i 1050 $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 : apizer 1080 assign_option_aliases(\%options,
336 :    
337 : sh002i 1050 );
338 :     set_default_options(\%options,
339 :     name => 'func',
340 :     llimit => $x_ref->[0],
341 :     rlimit => $x_ref->[$#$x_ref],
342 :     );
343 : apizer 1080
344 : sh002i 1050 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 : apizer 1080
366 : sh002i 1050 // Evaluate a cubic spline defined by the vectors above
367 :     i = 0;
368 :     while (x > t[i+1] ) {
369 :     i++
370 :     }
371 : apizer 1080
372 : sh002i 1050 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 : apizer 1080 + ( x - t[i] )*( b[i]*( x - t[i])*( x - t[i] ) +c[i] )
375 : sh002i 1050 );
376 :    
377 :     } else {
378 :     return( "undefined" ) ;
379 :     }
380 :    
381 :     }
382 : apizer 1080 // End
383 : sh002i 1050 -->
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 : apizer 1080 assign_option_aliases(\%options,
415 : sh002i 1050 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 : apizer 1080
465 : sh002i 1050 $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 : apizer 1080 Uses Newton's method of approximating roots of equations, and Romberg to evaluate definite integrals.
480 : sh002i 1050
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 : apizer 1080 my $i=0;
490 : sh002i 1050 my $funct;
491 :     my $deriv;
492 : apizer 1080 while (abs($delta) > 0.000001 && $i < 5000) {
493 : sh002i 1050 $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 : apizer 1080 }
499 : sh002i 1050 $delta = $funct/$deriv;
500 : apizer 1080 $b = $b - $delta;
501 : sh002i 1050 $i++;
502 : apizer 1080 }
503 : sh002i 1050 if ($i == 5000) {
504 :     warn 'Newtons method does not converge.';
505 : apizer 1080 return;
506 : sh002i 1050 }
507 :     $b;
508 : apizer 1080 }
509 : sh002i 1050
510 :     #########################################################
511 :    
512 :     1;
513 :    

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9