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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4997 - (download) (as text) (annotate)
Mon Jun 11 18:16:40 2007 UTC (12 years, 7 months ago) by gage
File size: 25640 byte(s)
Fixing docementation so that it can be read from the web.

    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   my $error_string = "This equation $correctEqn is not constant on solution curves of  y'(t) = $diffEqRHS\r\n<br>
  787                             starting at ( $initial_t , $initial_y )<br>
  788                             $check_eval->ans_hash->pretty_print()".
  789           "options<br>\n".pretty_print({  vars  =>  [@VARS],
  790                   params  =>  ['c'],
  791                   tolType =>  $options{tolType},
  792                   relTol  =>  $options{relTol},
  793                   tol   =>  $options{tol},
  794                   debug =>  $options{debug},
  795                     });
  796 
  797   for (my $i=0; $i<$numPoints;$i++) {
  798       my ($z, $err) = &$rf_corrEq( $out1->[$i][0], $out1->[$i][1] );
  799       $z = $err if defined $err;
  800       $error_string .= "F( ". $out1->[$i][0] . " , ". $out1->[$i][1] . " ) = $z <br>\r\n";
  801   }
  802   $error_string .= $rh_correct_ans->error_message();
  803   warn $error_string, $check_eval->ans_hash->pretty_print;
  804     }
  805 
  806     my ($constant_eval) = fun_cmp('c', vars => [@VARS],
  807           params  =>  ['c'],
  808           tolType =>  $options{tolType},
  809           relTol  =>  $options{relTol},
  810           tol   =>  $options{tol},
  811           debug =>  $options{debug},
  812           );  # an evaluator that tests for constants;
  813     $constant_eval->ans_hash(evaluation_points => $rh_correct_ans->{evaluation_points});
  814     my $answer_evaluator = new AnswerEvaluator;
  815     $answer_evaluator->ans_hash(  correct_ans     =>  $rh_correct_ans->{correct_ans},       # used for answer only
  816         rf_correct_ans    =>  sub { my @input = @_; pop(@input); },
  817         # return the last input which is the constant parameter 'c';
  818         evaluation_points =>  $rh_correct_ans->{evaluation_points},
  819         ra_param_vars     =>  ['c'],                                # compare with constant function
  820         ra_vars       =>  [@VARS],
  821         type        =>  'level_curve',
  822         );
  823     $answer_evaluator->install_evaluator(sub { my $ans_hash = shift;
  824                  my %options = @_;
  825                  $constant_eval->evaluate($ans_hash->{student_ans});
  826                  $constant_eval->ans_hash;
  827              });
  828 
  829     $answer_evaluator->install_post_filter( sub { my $ans_hash = shift; $ans_hash->{correct_ans} = $correctEqn; $ans_hash; } );
  830     $answer_evaluator->install_post_filter( sub {     my $rh_ans= shift;
  831                 my %options = @_;
  832                 if ($rh_ans->catch_error('SYNTAX') ) {
  833                     $rh_ans->{ans_message} = $rh_ans->{error_message};
  834                     $rh_ans->clear_error('SYNTAX');
  835 
  836                 }
  837                 $rh_ans;
  838                   });
  839 
  840   main::PG_restricted_eval('$main::useOldAnswerMacros = '.$saveUseOldAnswerMacros);
  841     $answer_evaluator;
  842 
  843 }
  844 
  845 
  846 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9