share/contrib/rand/improved.mac

   1 /* Filename improved.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    * Reverse engineered from file newimprv.bk1                   *
  14    * David Billinghurst - Feb 2004                               *
  15    *                                                             *
  16    ***************************************************************
  17 */ /*   This program uses recursive functions to find
  18 the transition curves in Mathieu's equation.  To call it,
  19 type:
  20                         TC()
  21 */
  22
  23 tc():=(input(),sign:1,find(),if n > 0 then (sign:-1,find()))$
  24 input():=(n:read("enter transition curve number n"),
  25       m:read("enter degree of truncation"))$
  26 find():=(remarray(b,e),delta:n^2/4,for i thru m do delta:delta+d(i)*e^i,
  27      print("delta=",delta),print(" "))$
  28 a(j,k):=
  29   if j < 0 or k < 0 then
  30     0
  31   else if j = 0 and k = n then
  32     1
  33   else if j = 0 then
  34     0
  35   else if k = n then
  36     0
  37   else if numberp(b[j,k]) then
  38     b[j,k]
  39   else if k = 0 then
  40     b[j,k]:(-a(j-1,2)/2-sum(d(i)*a(j-i,0),i,1,j))/(n^2/4)
  41   else
  42     b[j,k]:(-(a(j-1,k-2)+a(j-1,k+2)+sign*a(j-1,2-k))/2
  43       -sum(d(i)*a(j-i,k),i,1,j))/((n^2-k^2)/4)
  44 $
  45
  46 d(j):=
  47   if numberp(e[j]) then
  48     e[j]
  49   else if n = 0 then
  50     e[j]:-a(j-1,2)/2
  51   else
  52     e[j]:-(a(j-1,n-2)+a(j-1,n+2)+sign*a(j-1,2-n))/2
  53 $