Parent Directory
| 
Revision Log
moved PG modules and macro files from webwork-modperl to pg -sam
1 #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 |