share/contrib/rand/Newimprv.bk1

   1 /* Filename <name>.mac
   2
   3    ***************************************************************
   4    *                                                             *
   5    *                     <package name>                          *
   6    *                <functionality description>                  *
   7    *                                                             *
   8    *          from: Computer Algebra in Applied Math.            *
   9    *                   by Rand (Pitman,1984)                     *
  10    *                Programmed by Richard Rand                   *
  11    *      These files are released to the public domain          *
  12    *                                                             *
  13    ***************************************************************
  14 */ /*   This program uses recursive functions to find
  15 the transition curves in Mathieu's equation.  To call it,
  16 type:
  17                         TC()
  18 */
  19
  20 tc():=(input(),sign:1,find(),IF n > 0 THEN (sign:-1,find()))$
  21 input():=(n:read("ENTER TRANSITION CURVE NUMBER N"),
  22       m:read("ENTER DEGREE OF TRUNCATION"))$
  23 find():=(delta:n^2/4,for i thru m do delta:delta+d[i]*e^i,
  24      print("delta=",delta),print(" "))$
  25 a[j,k]:=IF j < 0 OR k < 0 THEN 0
  26    ELSE (IF j = 0 AND k = n THEN 1
  27              ELSE (IF j = 0 THEN 0
  28                        ELSE (IF k = n THEN 0
  29                                  ELSE (IF k = 0
  30                                            THEN (-a[j-1,2]/2
  31                                            -sum(d[i]*a[j-i,0],i,1,j))
  32                                            /(n^2/4)
  33                                            ELSE (-(a[j-1,k-2]
  34                                            +a[j-1,k+2]+sign*a[j-1,2-k])
  35                                            /2
  36                                            -sum(d[i]*a[j-i,k],i,1,j))
  37                                            /((n^2-k^2)/4)))))$
  38 d[j]:=IF n = 0 THEN -a[j-1,2]/2
  39    ELSE -(a[j-1,n-2]+a[j-1,n+2]+sign*a[j-1,2-n])/2$