[system] / trunk / webwork / system / courseScripts / PGnumericalmacros.pl Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

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

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
8
9   Numerical methods for the PG language
10
12
13
14
16
17 =cut
18
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
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
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
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
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
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
```