[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 1831 - (download) (as text) (annotate)
Thu Mar 4 13:06:32 2004 UTC (15 years, 11 months ago) by gage
File size: 30161 byte(s)
Changed documentation spacing

    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