share/contrib/rand/recursiv.mac

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