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