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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : sh002i 1050 #163macros.pl
2 :     #macros for Prills 163 problems
3 :    
4 :     #!/usr/bin/perl -w
5 :     #use strict;
6 :     #use Carp;
7 :     BEGIN {
8 :     be_strict();#all variables must be declared local or global
9 :     }
10 :    
11 :     #my @answer = oldivy(1,2,1,8,4);
12 :     #print ("The old program says:\n");
13 :     #print ($answer[0]);
14 :     #print ("\n");
15 :     #@answer = ivy(1,2,1,8,4);
16 :     #print ("My program says:\n");
17 :     #print ($answer[0]);
18 :     #print ("\n");
19 :    
20 :     #the subroutine is invoked with arguments such as complexmult(2,3,4,5) or
21 :     #complexmult(@data), where @data = (2,3,4,5)
22 :    
23 :     sub complexmult {
24 :     my ($S,$T,$U,$V) = @_;
25 :     #this line defines the input arguments as local variables
26 :     my $R = $S *$U -$T * $V;
27 :     my $I = $S *$V + $T * $U ;
28 :     ($R,$I) ;#this returns ($R,$I) from the subroutine
29 :     }
30 :     ##########
31 :     #sub addtwo adds two strings formally
32 :     #An "indicator" for a string is a
33 :     # number ,e.g. coefficient,which indicates
34 :     #whether the string is to be
35 :     #added or is to be regarded as zero.
36 :     #The non-zero terms are formally added as strings.
37 :     #The input is an array
38 :     #($1staddend, $1stindicator,$2ndaddend,$2ndindicator)
39 :     #The return is an array
40 :     #(formal sum, indicator of formal sum)
41 :    
42 :     sub addtwo {
43 :     my ($A,$a,$B,$b) = @_;
44 :     my $out = "0";
45 :     if ($a != 0) {
46 :     $out= $A;
47 :     }
48 :    
49 :    
50 :     if ($b != 0) {
51 :     if ($a == 0) {
52 :     $out= $B ;
53 :     } else {
54 :     $out = "$A + $B";
55 :     }
56 :     }
57 :     my $ind = abs($a) + abs($b);
58 :     ($out,$ind);
59 :     }
60 :     ########
61 :    
62 :     #sub add generalizes sub addtwo to more addends.
63 :     #It formally adds the nonzero terms.
64 :     # The input is an array of even length
65 :     # consisting of each addend,a string,
66 :     # followed by its indicator.
67 :    
68 :     sub add {
69 :     # this function takes the first two terms, puts them together, and keep repeating until you have emptied the list.
70 :     my @sum = ("0" ,0);
71 :     my @list = @_;
72 :     my $el = @list;
73 :     my $x = "0";
74 :     my $y = 0;
75 :     while ($el > 0) {
76 :     $x = shift(@list);
77 :     $y = shift(@list);
78 :     push(@sum,$x,$y);
79 :     @sum = addtwo(@sum);
80 :     $el = @list;
81 :     }
82 :     @sum ;
83 :     }
84 :     #######
85 :    
86 :     # sub diffop cleans up the typed expression
87 :     #of a diff. operator.
88 :     #input @diffop =($A,$B,$C) is the coefficients.
89 :     #input is given as arguments viz difftop($A,$B,$C);
90 :     #output is the diff. operator as a string $L in TEX
91 :    
92 :     sub diffop
93 :     {
94 :     my ($A,$B,$C) = @_ ;
95 :     my ($LDD, $LD ,$LF) = ($A."y'' ",$B."y' ",$C."y ");
96 :     # re-write 1y'' as y''.
97 :     if ($A==1){
98 :     $LDD = "y''";
99 :     }
100 :     # re-write 1y' as y'
101 :     if ($B==1) {
102 :     $LD = "y'";
103 :     }
104 :     # re-write -1y' as -y'
105 :     elsif ($B==-1) {
106 :     $LD = "-y'";
107 :     }
108 :     # re-write 1y as y
109 :     if ($C==1) {
110 :     $LF = "y";
111 :     }
112 :     # re-write -1y as -y
113 :     elsif ($C==-1) {
114 :     $LF = "-y";
115 :     }
116 :     my ($L,$ind) = add($LDD,$A,$LD,$B,$LF,$C);
117 :     $L;
118 :     }
119 :    
120 :     ########
121 :    
122 :     #sub rad simplifies (a/b)*(sqrt(c))
123 :     #input is given as arguments on rad viz.: rad($a,$b,$c);
124 :     #$a,$b,$c are integers and $c>=0 and $b is not zero.
125 :     #output is an array =(answer as string,new$a,new$b, new$c)
126 :    
127 :     sub rad{
128 :     # initalize primes
129 :     my @p = (2,3,5,7,11,13,17,19,23,29);
130 :     my ($a,$b,$c) = @_;
131 :     my $s = "0" ;
132 :     my $i = 0 ;
133 :     # if a=0 then re-write as (0/1)*(sqrt(1)) = 0.
134 :     if ($c*$a ==0){
135 :     $a = 0;
136 :     $b = 1;
137 :     $c = 1;
138 :     }
139 :    
140 :     # if b<0 then re-write the numerator as negative, denominator as positive
141 :     if ($b < 0){
142 :     $a = - $a ;
143 :     $b = - $b ;
144 :     }
145 :    
146 :     my $j = 1 ;
147 :     while($j == 1){
148 :     # can't reduce sqrt(1).
149 :     if ($c == 1){
150 :     last;
151 :     }
152 :     $j = 0;
153 :     foreach $i (@p){
154 :     # if you can factor a prime out of c, do it.
155 :     # ie, sqrt(75) = 5*sqrt(3)
156 :     if ( $c == $i*$i*int($c/($i*$i))){
157 :     $a = $a *$i;
158 :     $c = $c/($i*$i);
159 :     $j=1;
160 :     }
161 :     }
162 :     }
163 :     $j = 1;
164 :    
165 :     # reduce fraction is lowest terms.
166 :     while($j==1){
167 :     # if the denominator is 1, then we're set.
168 :     if ($b==1){
169 :     last;
170 :     }
171 :     $j = 0;
172 :     foreach $i (@p){
173 :     # if you can factor a prime out of both numerator and denominator, do it.
174 :     if ( abs($a) + $b == $i*int(abs($a) /$i) + $i*int($b/$i) ){
175 :     $a = $a /$i;
176 :     $b = $b /$i;
177 :     $j=1;
178 :     }
179 :     }
180 :     }
181 :    
182 :     # $s = answer string
183 :    
184 :     # if you have ($a/1)*sqrt(1) then the answer is simply "$a".
185 :     if ($c == 1) {
186 :     if ($b == 1){
187 :     $s = "$a";
188 :     }
189 :     # if you have ($a/$b)*sqrt(1) then the answer is "$a/$b".
190 :     else {
191 :     $s = "$a/$b";
192 :     }
193 :     }
194 :    
195 :     if ($c > 1) {
196 :     if ($a != 1){
197 :     if ($b == 1) {
198 :     # if denominator is 1, then answer is "$a*sqrt($c)".
199 :     $s = "$a * sqrt($c)";
200 :     }
201 :     # if denominator is nonzero... answer is all three terms.
202 :     else {
203 :     $s = "$a *sqrt($c) /$b";
204 :     }
205 :     } else {
206 :     # if you have "(1/1)*sqrt($c)" then the answer is "sqrt($c)".
207 :     if ($b == 1) {
208 :     $s = "sqrt($c)";
209 :     }
210 :     # if you have "(1/$b)*sqrt($c)" then answer is "sqrt($c)/$b".
211 :     else {
212 :     $s = "sqrt($c) /$b";
213 :     }
214 :     }
215 :     }
216 :     # return all four variables: answer string, reduced numerator, reduced denominator, squareroot with primes factored out
217 :     my $rh_hash = { displaystr => $s,
218 :     num => $a,
219 :     denom => $b,
220 :     root => $c};
221 :     $rh_hash;
222 :    
223 :     }
224 :    
225 :     #######
226 :     sub frac {
227 :     # use rad subroutine
228 :     my ($a,$b) = @_;
229 :     rad($a,$b,1);
230 :     }
231 :     ##########
232 :    
233 :     ####
234 :     #sub exp simplifies exp($r*t) in form for writing perl
235 :     #or tex.The input is exp($r,$ind); $ind indicates whether
236 :     #we want perl or tex mode.$r is a string that represents
237 :     #a number.
238 :     #If $ind = 0 output is "exp(($r)*t)", simplified if possible.
239 :     #If $ind = 1 output is "exp(($r) t)", simplified if possible.
240 :     sub simpleexp {
241 :     my $r = shift;
242 :     my @rr = @_;
243 :     my %options;
244 :     if ($rr[0] eq 'mode') {
245 :     $options{'mode'} = $rr[1];
246 :     } elsif ( $rr[0] == 1) {
247 :     $options{'mode'} = 'typeset';
248 :     } elsif ($rr[0] == 0 ) {
249 :     $options{'mode'} = 'standard';# means we use * for multiplication
250 :     }
251 :    
252 :     my $y = "0";
253 :     if ($r eq "0"){
254 :     if ($options{'mode'} eq 'standard' ) {
255 :     $y = "1";
256 :     } elsif ($options{'mode'} eq 'typeset' ) {
257 :     $y = ""; # multiplication by 1 is the identity
258 :     } else {
259 :     warn "simpleexp doesn't recognize the mode $options{'mode'}";
260 :     }
261 :    
262 :     # change exp(1t) = exp(t)
263 :     }elsif ($r eq "1"){
264 :     $y = "exp(t)";
265 :     }
266 :     # change exp(-1t) = exp(-t)
267 :     elsif ($r eq "-1"){
268 :     $y = "exp(-t)";
269 :     }
270 :     # in typeset modeset you don't use the *
271 :     # in standard modeset you use *
272 :     else {
273 :     if ($options{'mode'} eq 'typeset') {
274 :     $y = "exp(($r)t)";
275 :     } elsif ($options{'mode'} eq 'standard') {
276 :     $y = "exp(($r)*t)";
277 :     } else {
278 :     warn "simpleexp doesn't recognize the mode $options{'mode'}";
279 :     }
280 :     }
281 :     $y;
282 :     }
283 :    
284 :     sub ivy {
285 :     # $a*y'' + $b*y' + $c*y = 0. y(0)=$m. y'(0)=$n.
286 :     my ($a, $b, $c, $m, $n) = @_;
287 :     my $d = ($b*$b)-(4*$a*$c); # d is the discriminant
288 :     my $c1 = "0";
289 :     my $c2 = "0";
290 :     my $r1 = "0";
291 :     my $r2 = "0";
292 :     my $answer = "";
293 :    
294 :     # c1 = first coefficient, c2 = second coefficient, rr1 = e^r1, rr2 = e^r2.
295 :     # c1*rr1 + c2*rr2 = c1*e^r1 + c2*e^r2
296 :    
297 :     if ($d > 0) {
298 :     # y(t) = [m/2 + sqrt(d)*(2An+Bm)/(2d)]e^[t*(-B+sqrt(d))/(2A)] + [m/2 - sqrt(d)*(2An+Bm)/(2d)]e^[t*(-B-sqrt(d))/(2A)]
299 :     my $piece1 = frac($m,2);
300 :     my $piece2 = rad(2*$a*$n+$b*$m,2*$d,$d);
301 :     $c1 = "$piece1->{displaystr} + $piece2->{displaystr}"; # first coefficient: "m/2 + sqrt(d)*(2An+Bm)/(2d)"
302 :     $c2 = "$piece1->{displaystr} - $piece2->{displaystr}"; # second coefficient: "m/2 - sqrt(d)*(2An+Bm)/(2d)"
303 :     $piece1 = frac(-$b,2*$a); # find (-B/(2A)) in lowest terms
304 :     $piece2 = rad(1,2*$a,$d); # find (sqrt(B*B-4AC)/(2A)) in lowest terms
305 :     $r1 = "$piece1->{displaystr} + $piece2->{displaystr}"; # r1: (-B+sqrt(d))/(2A)
306 :     $r2 = "$piece1->{displaystr} - $piece2->{displaystr}"; # r2: (-B-sqrt(d))/(2A)
307 :     my $rr1 = simpleexp($r1,0); # raise e^r1.
308 :     my $rr2 = simpleexp($r2,0); # raise e^r2.
309 :     $answer = "($c2) *$rr2 + ($c1) *$rr1";
310 :    
311 :     }
312 :    
313 :     if ($d == 0) {
314 :     # y(t) = me^((-B/(2A)t) + [(n+mB)/(2A)]*t*e^((-Bt)/(2*A))
315 :     my $piece1 = frac(-$b,2*$a); # find (-B/(2A)) in lowest terms
316 :     my $piece2 = frac(2*$a*$n+$m*$b,2*$a); # find (2An+Bm)/(2A) in lowest terms
317 :     $c1 = "$m"; # first coefficient: "m"
318 :     $c2 = "$piece2->{displaystr}"; # second coefficient: "(n+mB)/(2A)"
319 :     $r1 = "$piece1->{displaystr}"; # r1: (-B/(2A))
320 :     $r2 = $r1; # r2: (-B/(2A))
321 :     my $rr1 = simpleexp($r1,0); # rr1 = e^r1 = e^(-B/(2A))
322 :     my $rr2 = simpleexp($r2,0); # rr2 = e^r2 = e^(-B/(2A))
323 :     $answer = "($c1) *$rr1 + ($c2)*t *$rr2";
324 :     }
325 :    
326 :     # if the descriminant is negative, then the roots are imaginary.
327 :     # recall, e^x where x=a+ib then e^x = (e^a)*cos(bt) + (e^a)*sin(bt).
328 :     if ($d<0){
329 :     # y(t) = me^(-Bt/(2A))*cos(t*sqrt(4AC-B*B)/(2A))+(2An+Bm)*sqrt(4AC-B*B)/(4AC-B*B)*e^(-Bt/(2A))*sin(t*sqrt(4AC-B*B)/(2A))
330 :     my $piece1 = rad (2*$a*$n+$b*$m,-$d,-$d); # find ((2An+Bm)*sqrt(4AC-B*B))/(4AC-B*B) in lowest terms
331 :     my $piece2 = rad (1,2*$a,-$d); # find (sqrt(4AC-B*B)/(2A)) in lowest terms
332 :     $c1 = "$m"; # first coefficient: "m"
333 :     $c2 = "$piece1->{displaystr}"; # second coefficient: "(2An+Bm)*sqrt(4AC-B*B)/(4AC-B*B)"
334 :     my $cs1 = "cos(($piece2->{displaystr})*t)"; # cos(t*sqrt(4AC-B*B)/(2A))
335 :     my $cs2 = "sin(($piece2->{displaystr})*t)"; # sin(t*sqrt(4AC-B*B)/(2A))
336 :     my $piece3 = frac (-$b,2*$a); # find (-B/(2A)) in lowest terms
337 :     $r1 = "$piece3->{displaystr}"; # r1: (-B/(2A))
338 :     $r2 = $r1; # r2: (-B/(2A))
339 :     my $rr1 = simpleexp($r1,0); # rr1 = e^r1 = e^(-B/(2A))
340 :     my $rr2 = simpleexp($r2,0); # rr2 = e^r2 = e^(-B/(2A))
341 :     $answer = "($c1) *($rr1)*($cs1) + ($c2) *($rr2)*($cs2)";
342 :     }
343 :     $answer;
344 :     }
345 :    
346 :     ############
347 :     #sub ivy solves the initial value problem
348 :     # $a*y'' + $b*y' + $c*y = 0, with y(0) = $m, y'(0) = $n
349 :    
350 :     #The numbers $a,$b,$c,$m,$n should be integers with $a not 0.
351 :     #The inputs are given as arguments viz: ivy ($a,$b,$c,$m,$n).
352 :     #The output is the solution as a string giving a function of t.
353 :    
354 :     #sub oldivy {
355 :     # # the two roots are: -b/(2a) plus/minus sqrt(b*b - 4*a*c)/(2a)
356 :     # my ($a,$b,$c,$m,$n) = @_;
357 :     # my $denom = 2*$a ;
358 :     # my $discriminant = ($b*$b)-(4*$a*$c);
359 :     # my $outa = frac(-$b,2*$a);
360 :     # my $dabs = abs($discriminant);
361 :     # my $ira = rad(1,2*$a,$dabs);
362 :     #
363 :     # # let $r = "-$b/(2*$a)" where the fraction is in lowest form.
364 :     # my $r = $outa->{displaystr};
365 :     #
366 :     # # make sure you take a square root of positive number.
367 :     # # let $w = sqrt(abs(-b/(2a)))
368 :     # my $w = sqrt($dabs);
369 :     #
370 :     # # let $irra = "sqrt($dabs)/(2a)" in lowest form.
371 :     # my $irra = $ira->{displaystr};
372 :     #
373 :     # my @suma = add($r,$b,$irra,$discriminant);
374 :     # my @sumb = add($r,$b,"-$irra",$discriminant);
375 :     #
376 :     # # $r1 is string of root1. $r2 is string of root2.
377 :     #
378 :     # # let $r2 = "$r + $irra" = "(-$b+sqrt($dabs))/(2a)".
379 :     # my $r2 = $suma[0];
380 :     #
381 :     # # let $r1 = "$r - $irra" = "(-$b-sqrt($dabs))/(2a)".
382 :     # my $r1 = $sumb[0];
383 :     #
384 :     # # if the sqrt(discriminant) is an integer (real solutions)...
385 :     # if ($w == int($w)) {
386 :     # my $outb = frac(-$b-$w,$denom);
387 :     # my $outc = frac(-$b+$w,$denom);
388 :     #
389 :     # # then let $r1 = (-$b-sqrt($dabs))/(2a). the "plus" term.
390 :     # $r1 = $outb->{displaystr};
391 :     #
392 :     # # then let $r2 = (-$b+sqrt($dabs))/(2a). the "minus" term.
393 :     # $r2 = $outc->{displaystr};
394 :     # }
395 :     #
396 :     #######
397 :     # # let $exp = "exp( (-$b)/(2a) t)";
398 :     # my $exp = simpleexp ($r ,0 );
399 :     #
400 :     # #first find a fundamental set of solutions.
401 :     # #initialize:
402 :     #
403 :     # # Solution: Ae^(($r1)t) + Be^(($r2)t) = y(t)
404 :     # # $y1 = e^(($r1)t), $y2 = e^(($r2)t)
405 :     # # $y1prime = ($r1)*e^(($r1)t), $y2prime = ($r2)*e^(($r2)t)
406 :     #
407 :     # my $y1 = "0";
408 :     # my $y2 = "0";
409 :     # my $y1prime = "0";
410 :     # my $y2prime = "0";
411 :     #
412 :     # if ($discriminant==0) {
413 :     # $y1 = simpleexp ($r,0);
414 :     # $y2 = "t* $y1 ";
415 :     # $y1prime = "$r *$y1";
416 :     # $y2prime = "(1 + $r*t)*$y1 ";
417 :     # if ($r eq "0") {
418 :     # $y1prime = "0";
419 :     # $y2prime = "1";
420 :     # $y2 = "t";
421 :     # }
422 :     # }
423 :     #
424 :     # # if the discriminant is positive,
425 :     # if ($discriminant>0) {
426 :     # $y1 = simpleexp ($r1,0);
427 :     # $y2 = simpleexp ($r2,0);
428 :     # $y1prime = "($r1)*$y1";
429 :     # $y2prime = "($r2)*$y2";
430 :     # }
431 :     #
432 :     # # if the discriminant is negative, then the roots are imaginary.
433 :     # # recall e^(ix) = cos x + isin x
434 :     # # and e^(ix) + e^(-ix) = 2cos x
435 :     #
436 :     # if ( $discriminant<0) {
437 :     # my $arg = "($irra )*t";
438 :     # if ($w == $denom) {
439 :     # $arg = "t";
440 :     # }
441 :     # my $cos = "cos($arg ) ";
442 :     # my $sin = "sin($arg ) ";
443 :     #
444 :     # if ($b == 0) {
445 :     # $y1 = "$cos";
446 :     # $y2 = "$sin";
447 :     # $y1prime = "(-$irra)*$y2";
448 :     # $y2prime = "($irra)*$y1";
449 :     # } else {
450 :     # $y1 = "($exp)*($cos) ";
451 :     # $y2 = "($exp)*($sin) ";
452 :     # $y1prime = "(-$irra)*($y2) + ($r)*($y1)";
453 :     # $y2prime = "($irra)*($y1) + ($r)*($y2)";
454 :     # }
455 :     # }
456 :     # print ("discriminant=$discriminant\n");
457 :     #
458 :     ######Need to get coef below to solve iv problem
459 :     #### Used Cramer's rule
460 :     ##initialize
461 :     # my $coef1 = "0";
462 :     # my $coef2 = "0";
463 :     # my $t1 = 0;
464 :     # my $t2 = 0;
465 :     # my $out = frac(2,1);
466 :     #
467 :     # my $k = ($m*$b*.5)+($a*$n);
468 :     # my $kk = ($m*$b)+(2*$a*$n);
469 :     #
470 :     # # if the discriminant is zero, then you have twin roots. ($r1=$r2)
471 :     # # so, Ae^(($r1)t)+Bte^(($r2)t) = Ae^(($r1)0)+B*0*e^(($r2)0) = y(0) = A+Bt = A = m.
472 :     #
473 :     # if ($discriminant==0) {
474 :     # $coef1 = "$m";
475 :     # #$coef2 = "$kk /$denom";
476 :     # $out = frac($kk, $denom);
477 :     # $coef2 = $out->{displaystr};
478 :     # $t1 = $m;
479 :     # $t2 = $kk;
480 :     # } elsif ($discriminant>0) {
481 :     # my $coef1numb = $w*$m - $kk;
482 :     # my $coef2numb = $w*$m + $kk;
483 :     # $t1 = $coef1numb;
484 :     # $t2 = $coef2numb;
485 :     # my $teiler = 2*$w;
486 :     # my $test = int($w);
487 :     # if ($w == $test) {
488 :     # #$coef1 = "$coef1numb / $teiler";
489 :     # $out = frac($t1,$teiler);
490 :     # $coef1 = $out->{displaystr};
491 :     # #$coef2 = "$coef2numb / $teiler";
492 :     # $out = frac($t2,$teiler);
493 :     # $coef2 = $out->{displaystr};
494 :     # } else {
495 :     # my $dd = 2*$discriminant;
496 :     # my $irc = rad($kk,$dd,$discriminant);
497 :     # my $irrc = $irc->{displaystr};
498 :     # my $mfirst = frac($m,2);
499 :     # my $mf = $mfirst->{displaystr} ;
500 :     # my @sumc = add($mf,$m,$irrc,$kk);
501 :     # my @sumd = add($mf,$m,"-($irrc)",$kk);
502 :     #
503 :     # # coef1 = coefficient of the first solution, coef2 = coefficient of second solution
504 :     #
505 :     # # let $coef2 = "$mf+$kk/(2*sqrt($discriminant))";
506 :     # $coef2 = $sumc[0];
507 :     #
508 :     # # let $coef1= "$mf-$kk/(2*sqrt($discriminant))";
509 :     # $coef1 = $sumd[0];
510 :     # }
511 :     # } else {
512 :     #
513 :     # # $coef2 = ($kk)/sqrt($dabs) = ($kk/$dabs)*sqrt($dabs)
514 :     # $coef1 = "$m";
515 :     # my $ire = rad($kk,$dabs,$dabs);
516 :     # $coef2 = $ire->{displaystr};
517 :     # $t1 = $m;
518 :     # $t2 = $kk;
519 :     # if ($w ==int($w)) {
520 :     # my $outd = frac($kk,$w);
521 :     # # let $coef2 = "$kk/$w";
522 :     # $coef2 = $outd->{displaystr};
523 :     # }
524 :     # }
525 :     #
526 :     # my $z1 = "($coef1 ) *$y1";
527 :     # if ($y1 eq "1") {
528 :     # $z1 = $coef1;
529 :     # }
530 :     #
531 :     # my $z2 = "($coef2 ) *$y2";
532 :     #
533 :     # # if you have "($coef2)*1" then just write "$coef2".
534 :     # if ($y2 eq "1") {
535 :     # $z2 = $coef2;
536 :     # }
537 :     #
538 :     # # $z1 = Ae^(($r1)t), $z2 = Be^(($r2)t)
539 :     # # $t1 = numerical value of A, $t2 = numerical value of B
540 :     #
541 :     # my @sume = add($z1,$t1,$z2,$t2);
542 :     # my $y = $sume[0] ;
543 :     # $y;
544 :     #}
545 :     ####
546 :    
547 :     sub undeterminedExp {
548 :     my ($A,$B,$C,$r,$q0,$q1,$q2) = @_;
549 :     my $P = "$A*x*x + $B *x + $C "; #characteristic poly.
550 :     my $PP = "$A*2*x + $B ";#derivative of characteristic poly.
551 :     my $exp = simpleexp($r,0);
552 :     #$Pr = $P;
553 :     #$Pr =~ s~~x~~$r~~g;
554 :     #$Pr = PG_answer_eval("$Pr ");
555 :     #$PPr = $PP;
556 :     #$PPr =~ s~~x~~$r~~g;
557 :     #$PPr = PG_answer_eval("$PPr ");
558 :     my $Pr = $A *$r *$r + $B *$r + $C ;
559 :     my $PPr = 2* $A* $r + $B;
560 :     #################
561 :     if ($Pr != 0){
562 :     #$r not a root of $P
563 :     # $v0 = "($q0 /$Pr) ";
564 :     my $n1 = -$q1 * $PPr ;
565 :     my $d1 = $Pr * $Pr;
566 :     # $v1 = " ($q1 /$Pr)* t + ($n1 / $d1) ";
567 :     my $n22 = $Pr *$Pr *$q2;
568 :     my $n21 = (- 2 *$PPr ) *$q2;
569 :     my $n20 = (2*$PPr *$PPr - ($A*2*$Pr)) *$q2 ;
570 :     my $d2 = $Pr **3;
571 :     # $v2 = "($n22/$d2) *t*t + ($n21/$d1) *t + ($n20/$d2) ";
572 :     my $c00n = $q0*$d1 -$q1 *$PPr *$Pr + $n20;
573 :     my $c00d = $d2;
574 :     my $fraca = frac ($c00n ,$c00d );
575 :     my $c00 = $fraca->{displaystr};
576 :     my $c01n = $Pr *$q1 - 2 *$PPr *$q2;
577 :     my $c01d = $d1;
578 :     my $fracb = frac($c01n ,$c01d );
579 :     my $c01 = $fracb->{displaystr};
580 :     $c01 = "($c01) *t";
581 :     my $c02n = $q2;
582 :     my $c02d = $Pr;
583 :     my $fracc = frac($c02n ,$c02d );
584 :     my $c02 = $fracc->{displaystr};
585 :     $c02 = "($c02 ) *t*t";
586 :     my @adda = add ($c00,$c00n,$c01,$c01n,$c02,$c02n);
587 :     my $outa = $adda [0];
588 :     my $y = "( $outa ) * $exp " ;
589 :     #############
590 :     } elsif ($PPr != 0){
591 :     #$r a simple root of $P
592 :     my $q2m = -$q2;
593 :     #$v0 = "( $q0/$PPr)*t ";
594 :     my $q2d = 2*$q2;
595 :     my $d10 = 2*$PPr;
596 :     my $d11 = $PPr **2 ;
597 :     my $q1m = -$q1;
598 :     #$v1 = "($q1 / $d10)*t*t + ($q1m / $d11)*t ";
599 :     my $d23 = 3*$PPr;
600 :     my $d22 = $PPr *$PPr ;
601 :     my $d21 = $PPr *$d22;
602 :     #$v2 = "($q2 / $d23) *t*t*t + ($q2m/$d22) *t*t + ($q2d/ #$d21) *t ";
603 :     ######
604 :     my $c10n = $q0 *$d22 -$q1*$A *$PPr + $A*$A*$q2d;
605 :     my $c10d = $d21;
606 :     my $fracd = frac($c10n ,$c10d );
607 :     my $c10 = $fracd->{displaystr};
608 :     #warn " c10 $c10";
609 :     $c10 = "($c10 )*t";
610 :     my $c11n = $PPr * $q1 - 2 *$A * $q2;
611 :     my $c11d = 2*$PPr*$PPr;
612 :     my $frace = frac($c11n ,$c11d );
613 :     my $c11 = $frace->{displaystr};
614 :     #warn " c11 $c11";
615 :     $c11 = "($c11 )*t*t";
616 :     my $c12n = $q2;
617 :     my $c12d = 3*$PPr;
618 :     my $fracf = frac ($c12n ,$c12d );
619 :     my $c12 = $fracf->{displaystr};
620 :     $c12 = "($c12 ) *t*t*t";
621 :     my @addb = add($c10,$c10n,$c11,$c11n,$c12,$c12n);
622 :     my $outb = $addb[0];
623 :     my $y = "( $outb ) * $exp" ;
624 :     ######
625 :     } else {
626 :     # $v2 = "($q2 /12*$A) *t*t*t*t ";
627 :     #v1 = "($q1 /6*$A) *t*t*t ";
628 :     #$v0 = "($q0 /2*$A) *t*t " ;
629 :     my $c20n = $q0;
630 :     my $c20d = 2*$A;
631 :     my $fracg = frac($q0 ,$c20d );
632 :     my $c20 = $fracg->{displaystr};
633 :     $c20 = "($c20 ) *t*t";
634 :     my $c21n = $q1;
635 :     my $c21d = 6*$A;
636 :     my $frach = frac($c21n ,$c21d );
637 :     my $c21 = $frach->{displaystr};
638 :     $c21 = "($c21) *t*t*t";
639 :     my $c22n = $q2;
640 :     my $c22d = 12*$A;
641 :     my $fraci = frac($c22n ,$c22d );
642 :     my $c22 = $fraci->{displaystr};
643 :     $c22 = "($c22 ) *t*t*t*t";
644 :     my @addc = add($c20,$c20n,$c21,$c21n,$c22,$c22n);
645 :     my $outc = $addc[0];
646 :     my $y = "( $outc ) * $exp" ;
647 :     }
648 :    
649 :     }
650 :     #################
651 :    
652 :     # undeterminedSin is a subroutine to solve
653 :     #undetermined coefficient problems that have
654 :     #sines and cosines.
655 :     #The input is an array ($A,$B,$C,$r,$w,$q1,$q0,$r1,$r0)
656 :     #given as arguments on undeterminedSin
657 :     # $L =$A y'' + $B y' + $C y
658 :     # $rhs = ($q1 t + $q0) cos($w t)exp($r t) +
659 :     # ($r1 t + $r0) sin($w t)exp($r t)
660 :     #The subroutine uses undetermined coefficients
661 :     #to find a solution $y of $L = $rhs .
662 :     #The output \is $y
663 :    
664 :     sub undeterminedSin {
665 :     my ($A,$B,$C,$r,$w,$q1,$q0,$r1,$r0) = @_;
666 :     my $Pr = ($A*$r*$r) + $B *$r + $C;
667 :     my $PPr = (2*$A *$r ) + $B ;
668 :     my $re = $Pr -$A* $w*$w ;
669 :     my $im = $PPr * $w ;
670 :     #If P(x) = A x^2 +Bx +C,
671 :     #P($r+i*$w)=$re+i $im.
672 :     my $D = $re **2 + $im **2 ;
673 :     # $D = |P($r + i*$w)|^2
674 :     my $reprime = $PPr;
675 :     my $imprime = 2*$A*$w;
676 :     #If PP(x) = derivative of P(x),
677 :     #PP($r+i $w)=$reprime+$imprime.
678 :     my $cos = "cos($w *t)";
679 :     my $sin = "sin($w *t)";
680 :     if ($w == 1){
681 :     $cos = "cos(t)";
682 :     $sin = "sin(t)";
683 :     }
684 :    
685 :    
686 :     my $exp = simpleexp($r,0);
687 :    
688 :     ############
689 :     if ($D != 0){
690 :     #We first handle case that$r+i$w not a root of $P(x)
691 :     #This solution is based on:
692 :     #Let L[y] = Ay'' +By'+C,z=r+iw,p=P(z),q=P'(z),S,T complex;p!=0.
693 :     #Then L[((St+T-(Sq/p))/p)*e^{zt}]=(St+T)e^{zt}.
694 :     #Take S=q1-i*r1;T=q0-i*r0.Then take real part.
695 :     my ($renum1 ,$imnum1) = complexmult($q1,-$r1,$re, -$im);
696 :    
697 :     #S*(p conjugate)= $renum1+i imnum1
698 :     my ($renum2 ,$imnum2) = complexmult ($renum1,$imnum1,$reprime,$imprime);
699 :     #The above are real and imag parts of q*S*(p conjugate)
700 :     my $first = ($D *$q0 ) - $renum2 ;
701 :     my $second = (-$r0 *$D ) -$imnum2 ;
702 :     my ($renum3, $imnum3 ) = complexmult( $first,$second,$re,-$im);
703 :     #these are re and im of (DT-qS(p conj))*(p conj.)
704 :     my $n1 = $renum1;
705 :     my $n2 = $renum3;
706 :     my $n3 = -$imnum1;
707 :     my $n4 = -$imnum3;
708 :     my $fraca = frac($n1,$D );
709 :     my $tcosp = $fraca->{displaystr};
710 :     #$tcospart = "($tcosp )*t*$exp *$cos ";
711 :     my $tcospart = "($tcosp)*t*$exp *$cos " ; #####################
712 :     my $DD = $D *$D;
713 :     my $fracb = frac($n2 , $DD );
714 :     my $cospart = $fracb->{displaystr};
715 :     $cospart = "($cospart )*$exp*$cos";
716 :     my $fracc = frac($n3 , $D );
717 :     my $tsinpart = $fracc->{displaystr};
718 :     $tsinpart = "($tsinpart )*t*$exp*$sin";
719 :     my $fracd = frac($n4 , $DD );
720 :     my $sinpart = $fracd->{displaystr};
721 :     $sinpart = "($sinpart )*$exp*$sin";
722 :     my @suma = add($tcospart,$n1,$cospart,$n2,$tsinpart,$n3,$sinpart,$n4 );
723 :     my $out = $suma[0];
724 :     #The solution is $out
725 :     } else{
726 :     #We now handle case that $r+iw is a root of $P
727 :     #In this case $PPr = 0 and $PP($r + i$w) = 2*$A*i*$w
728 :     #Solution based on
729 :     #L[((S/2q)t*t -(AS/q*q)t +(T/q)t)e^{zt}]=
730 :     #(St+T)e^{zt}.Notation as for 1st part.
731 :     my $n3 = $q1 - (2*$w *$r0);
732 :     my $n4 = $r1 + (2*$w *$q0 );
733 :     my $n1 = -$r1 ;
734 :     my $n2 = $q1 ;
735 :     my $T2 = 4*$A *$w ;
736 :     my $T1 = $w * $T2 ;
737 :     my $frace = frac($n1 , $T2 );
738 :     my $t2cospart = $frace->{displaystr};
739 :     $t2cospart = "($t2cospart )*t*t*$exp *$cos ";
740 :     my $fracf = frac($n3 , $T1 );
741 :     my $tcospart = $fracf->{displaystr};
742 :     $tcospart = "($tcospart )*t*$exp *$cos ";
743 :     my $fracg = frac($n2 , $T2 );
744 :     my $t2sinpart = $fracg->{displaystr};
745 :     $t2sinpart = "($t2sinpart )*t*t*$exp*$sin";
746 :     my $frach = frac($n4 , $T1 );
747 :     my $tsinpart = $frach->{displaystr};
748 :     $tsinpart = "($tsinpart )*t*$exp*$sin";
749 :     my @addb = add ($t2cospart,$n1,$tcospart,$n3,$t2sinpart,$n2,$tsinpart,$n4 );
750 :     my $out = $addb[0];
751 :     }
752 :    
753 :     }
754 :    
755 :     sub check_eigenvector {
756 :    
757 :     my $eigenvalue = shift;
758 :     my $matrix = shift;
759 :     my %options = @_;
760 :     assign_option_aliases( \%options, );
761 :    
762 :     set_default_options( \%options,
763 :     'debug' => 0,
764 :     'correct_ans' => undef
765 :     );
766 :    
767 :    
768 :     my @correct_vector = ();
769 :     @correct_vector = @{$options{'correct_ans'}} if defined ($options{'correct_ans'});
770 :    
771 :     my $ans_eval = new AnswerEvaluator;
772 :    
773 :     $ans_eval->{debug} = $options{debug};
774 :     my $corr_ans_points = "( " . join(", ", @correct_vector). " )" ;
775 :     $ans_eval->ans_hash( correct_ans => $corr_ans_points );
776 :     $ans_eval->install_pre_filter(\&is_array);
777 :     $ans_eval->install_pre_filter(\&std_num_array_filter);
778 :    
779 :     $ans_eval->install_evaluator(sub { my $rh_ans = shift;
780 :     my %options = @_;
781 :     my @vector = @{$rh_ans->input()};
782 :     return($rh_ans) unless @correct_vector == @vector;
783 :     # make sure the vectors are the same dimension
784 :    
785 :     my $vec = new Matrix(2,1);
786 :     $vec->assign(1,1, $vector[0]);
787 :     $vec->assign(2,1, $vector[1]);
788 :     my $out_vec = $matrix * $vec;
789 :     my @diff;
790 :     $diff[0] = $out_vec->element(1,1) - $vec->element(1,1)*$eigenvalue;
791 :     $diff[1] = $out_vec->element(2,1) - $vec->element(2,1)*$eigenvalue;
792 :     $rh_ans->{score} = zero_check(\@diff);
793 :     $rh_ans;
794 :    
795 :     });
796 :     $ans_eval->install_post_filter( sub { my $rh_ans= shift;
797 :     if ($rh_ans->error_flag('SYNTAX') ) {
798 :     $rh_ans->{ans_message} = $rh_ans->{error_message};
799 :     $rh_ans->clear_error('SYNTAX');
800 :     $rh_ans;
801 :     }
802 :     });
803 :    
804 :     $ans_eval;
805 :     }
806 :    
807 :    
808 :    
809 :     =pod
810 :    
811 :     rungeKutta4a
812 :    
813 :    
814 :     =cut
815 :    
816 :    
817 :    
818 :     sub rungeKutta4a {
819 :     my $rh_ans = shift;
820 :     my %options = @_;
821 :     my $rf_fun = $rh_ans->{rf_diffeq};
822 :     set_default_options( \%options,
823 :     'initial_t' => 1,
824 :     'initial_y' => 1,
825 :     'dt' => .01,
826 :     'num_of_points' => 10, #number of reported points
827 :     'interior_points' => 5, # number of 'interior' steps between reported points
828 :     'debug' => 1, # remind programmers to always pass the debug parameter
829 :     );
830 :     my $t = $options{initial_t};
831 :     my $y = $options{initial_y};
832 :    
833 :     my $num = $options{'num_of_points'}; # number of points
834 :     my $num2 = $options{'interior_points'}; # number of steps between points.
835 :     my $dt = $options{'dt'};
836 :     my $errors = undef;
837 :     my $rf_rhs = sub { my @in = @_;
838 :     my ( $out, $err) = &$rf_fun(@in);
839 :     $errors .= " $err at ( ".join(" , ", @in) . " )<br>\n" if defined($err);
840 :     $out = 'NaN' if defined($err) and not is_a_number($out);
841 :     $out;
842 :     };
843 :    
844 :     my @output = ([$t, $y]);
845 :     my ($i, $j, $K1,$K2,$K3,$K4);
846 :    
847 :     for ($j=0; $j<$num; $j++) {
848 :     for ($i=0; $i<$num2; $i++) {
849 :     $K1 = $dt*&$rf_rhs($t, $y);
850 :     $K2 = $dt*&$rf_rhs($t+$dt/2,$y+$K1/2);
851 :     $K3 = $dt*&$rf_rhs($t+$dt/2, $y+$K2/2);
852 :     $K4 = $dt*&$rf_rhs($t+$dt, $y+$K3);
853 :     $y = $y + ($K1 + 2*$K2 + 2*$K3 + $K4)/6;
854 :     $t = $t + $dt;
855 :     }
856 :     push(@output, [$t, $y]);
857 :     }
858 :     $rh_ans->{evaluation_points} = \@output;
859 :     $rh_ans->throw_error($errors) if defined($errors);
860 :     $rh_ans;
861 :     }
862 :    
863 :    
864 :     sub level_curve_check {
865 :     my $diffEqRHS = shift; #required differential equation
866 :     my $correctEqn = shift; # required answer in order to check the equation
867 :     my %options = @_;
868 :     assign_option_aliases( \%options,
869 :     'vars' => 'var',
870 :     'numPoints' => 'num_of_points',
871 :     'reltol' => 'relTol',
872 :     );
873 :     set_default_options( \%options,
874 :     'initial_t' => 0,
875 :     'initial_y' => 1,
876 :     'var' => [qw( x y )],
877 :     'num_of_points' => 10,
878 :     'tolType' => (defined($options{tol}) ) ? 'absolute' : 'relative',
879 :     'relTol' => .01,
880 :     'tol' => .01,
881 :     'debug' => 0,
882 :     );
883 :    
884 :     my $initial_t = $options{initial_t};
885 :     my $initial_y = $options{initial_y};
886 :     my $var = $options{var};
887 :     my $numPoints = $options{num_of_points};
888 :     my @VARS = get_var_array( $var );
889 :     my ($tolType, $tol);
890 :     if ($options{tolType} eq 'absolute') {
891 :     $tolType = 'absolute';
892 :     $tol = $options{'tol'};
893 :     delete($options{'relTol'}) if exists( $options{'relTol'} );
894 :     } else {
895 :     $tolType = 'relative';
896 :     $tol = $options{'relTol'};
897 :     delete($options{'tol'}) if exists( $options{'tol'} );
898 :     }
899 :     #prepare the correct answer and check its syntax
900 :     my $rh_correct_ans = new AnswerHash;
901 :     $rh_correct_ans ->{correct_ans} = $correctEqn;
902 :     # check and calculate the function defining the differential equation
903 :     $rh_correct_ans->input( $diffEqRHS );
904 :     $rh_correct_ans = check_syntax($rh_correct_ans);
905 :     warn $rh_correct_ans->{error_message},$rh_correct_ans->pretty_print() if $rh_correct_ans->{error_flag};
906 :    
907 :     $rh_correct_ans->{error_flag} = undef;
908 :    
909 :     $rh_correct_ans = function_from_string2($rh_correct_ans,
910 :     ra_vars => [@VARS],
911 :     store_in =>'rf_diffeq',
912 :     debug=>$options{debug}
913 :     );
914 :     warn "Error in compiling instructor's answer: $diffEqRHS<br> $rh_correct_ans->{error_message}<br>\n$rh_correct_ans->pretty_print()"
915 :     if $rh_correct_ans->{error_flag};
916 :    
917 :    
918 :     # create the test points that should lie on a solution curve of the differential equation
919 :     $rh_correct_ans = rungeKutta4a( $rh_correct_ans,
920 :     initial_t => $initial_t,
921 :     initial_y => $initial_y,
922 :     num_of_points => $numPoints,
923 :     debug=>$options{debug}
924 :     );
925 :     warn "Errors in calculating the solution curve $rh_correct_ans->{student_ans}<BR>\n
926 :     $rh_correct_ans->{error_message}<br>\n",$rh_correct_ans->pretty_print() if $rh_correct_ans->catch_error();
927 :     $rh_correct_ans->clear_error();
928 :    
929 :     # check and compile the correct answer submitted by the instructor.
930 :     my ($check_eval) = fun_cmp('c', vars => [@VARS],
931 :     params => ['c'],
932 :     tolType => $options{tolType},
933 :     relTol => $options{relTol},
934 :     tol => $options{tol},
935 :     debug => $options{debug},
936 :     ); # an evaluator that tests for constants;
937 :     $check_eval->ans_hash(evaluation_points => $rh_correct_ans->{evaluation_points});
938 :     $check_eval->evaluate($rh_correct_ans->{correct_ans});
939 :     if( $check_eval->ans_hash->{score} == 0 or (defined($options{debug}) and $options{debug})) {
940 :     # write error message for professor
941 :     my $out1 = $check_eval->ans_hash->{evaluation_points};
942 :     my $rf_corrEq = $check_eval->ans_hash->{rf_student_ans};
943 :     my $error_string = "This equation $correctEqn is not constant on solution curves of y'(t) = $diffEqRHS\r\n<br>
944 :     starting at ( $initial_t , $initial_y )<br>
945 :     $check_eval->ans_hash->pretty_print()".
946 :     "options<br>\n".pretty_print({ vars => [@VARS],
947 :     params => ['c'],
948 :     tolType => $options{tolType},
949 :     relTol => $options{relTol},
950 :     tol => $options{tol},
951 :     debug => $options{debug},
952 :     });
953 :    
954 :     for (my $i=0; $i<$numPoints;$i++) {
955 :     my ($z, $err) = &$rf_corrEq( $out1->[$i][0], $out1->[$i][1] );
956 :     $z = $err if defined $err;
957 :     $error_string .= "F( ". $out1->[$i][0] . " , ". $out1->[$i][1] . " ) = $z <br>\r\n";
958 :     }
959 :     $error_string .= $rh_correct_ans->error_message();
960 :     warn $error_string, $check_eval->ans_hash->pretty_print;
961 :     }
962 :    
963 :     my ($constant_eval) = fun_cmp('c', vars => [@VARS],
964 :     params => ['c'],
965 :     tolType => $options{tolType},
966 :     relTol => $options{relTol},
967 :     tol => $options{tol},
968 :     debug => $options{debug},
969 :     ); # an evaluator that tests for constants;
970 :     $constant_eval->ans_hash(evaluation_points => $rh_correct_ans->{evaluation_points});
971 :     my $answer_evaluator = new AnswerEvaluator;
972 :     $answer_evaluator->ans_hash( correct_ans => $rh_correct_ans->{correct_ans}, # used for answer only
973 :     rf_correct_ans => sub { my @input = @_; pop(@input); },
974 :     # return the last input which is the constant parameter 'c';
975 :     evaluation_points => $rh_correct_ans->{evaluation_points},
976 :     ra_param_vars => ['c'], # compare with constant function
977 :     ra_vars => [@VARS],
978 :     type => 'level_curve',
979 :     );
980 :     $answer_evaluator->install_evaluator(sub { my $ans_hash = shift;
981 :     my %options = @_;
982 :     $constant_eval->evaluate($ans_hash->{student_ans});
983 :     $constant_eval->ans_hash;
984 :     });
985 :    
986 :     $answer_evaluator->install_post_filter( sub { my $ans_hash = shift; $ans_hash->{correct_ans} = $correctEqn; $ans_hash; } );
987 :     $answer_evaluator->install_post_filter( sub { my $rh_ans= shift;
988 :     my %options = @_;
989 :     if ($rh_ans->catch_error('SYNTAX') ) {
990 :     $rh_ans->{ans_message} = $rh_ans->{error_message};
991 :     $rh_ans->clear_error('SYNTAX');
992 :    
993 :     }
994 :     $rh_ans;
995 :     });
996 :    
997 :     $answer_evaluator;
998 :    
999 :     }
1000 :    
1001 :    
1002 :     1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9