share/contrib/rand/takens.mac

   1 /* Filename takens.mac
   2
   3    ***************************************************************
   4    *                                                             *
   5    *                     <package name>                          *
   6    *                <functionality description>                  *
   7    *                                                             *
   8    *          from: "Determinacy of Degenerate Equilibria"       *
   9    *                  by Rand & Keith  Applied Mathematics       *
  10    *                 and Computation 21:1-19 (1987)              *
  11    *                                                             *
  12    *                Programmed by Richard Rand                   *
  13    *      These files are released to the public domain          *
  14    *                                                             *
  15    ***************************************************************
  16 */
  17 /* takens' blow-up calculation */
  18
  19 /* main program */
  20 takens():=block(
  21 /* variable i tags loop */
  22 for i:1 thru 8 do
  23         (setup1(),
  24         if i=1 then inputrhs() ,
  25         if i>1 then blowup(),
  26         setup2(),
  27         deffg(),
  28         print(" takens' test "),
  29         print(" truncate f and g to homogeneous polynomials "),
  30         print(takens_truncate()),
  31         getroots(),
  32         print(test()),
  33         if flag=green then return("done"),
  34         defpq(),
  35         defppqq(),
  36         sroots()) )$
  37
  38 /* subroutines to create variable names at ith loop */
  39 setup1():=(
  40 u:concat('u,i),
  41 v:concat('v,i),
  42 x:concat('x,i),
  43 y:concat('y,i))$
  44
  45 setup2():=(
  46 f:concat('f,i),
  47 g:concat('g,i),
  48 p:concat('p,i),
  49 q:concat('q,i),
  50 r:concat('r,i),
  51 s:concat('s,i),
  52 pp:concat('pp,i),
  53 qq:concat('qq,i))$
  54
  55 /* subroutine to input the rhs's from keyboard */
  56 inputrhs():=(
  57 print(" enter the rhs's to be studied "),
  58 print(" use variables x,y, they will be converted to x1,y1 "),
  59 u::read(u,"="),
  60 print(ev(u)),
  61 v::read(v,"="),
  62 print(ev(v)))$
  63
  64 /* subroutine to truncate f and g to terms of lowest degree */
  65 takens_truncate():=block(
  66 for j from 2 thru 8 do
  67         (temp1:ratexpand([ev(f),ev(g)]),
  68         temp2:taylor(temp1,[ev(x),ev(y)],0,j),
  69         if temp2 # taylor([0,0],dummy,0,1) then return(temp2)))$
  70
  71 /* subroutine to solve gtrunc = 0 */
  72 getroots():=(
  73 print("solving gtrunc = 0"),
  74 ftrunc:part(temp2,1),
  75 gtrunc:part(temp2,2),
  76 gtruncx:diff(gtrunc,ev(x)),
  77 gtruncy:diff(gtrunc,ev(y)),
  78 xroots:solve(gtrunc,ev(x)),
  79 yroots:solve(gtrunc,ev(y)),
  80 rootnum:0,
  81 for k:1 thru length(xroots) do
  82         (rootnum:rootnum+1,
  83         root[rootnum]:part(xroots,k)),
  84 for k:1 thru length(yroots) do
  85         (rootnum:rootnum+1,
  86         root[rootnum]:part(yroots,k)),
  87 print("total no. of roots =",rootnum))$
  88
  89 /* perform takens' test for each root */
  90 test():=(block(
  91 flag:green,
  92 for k:1 thru rootnum do (print(root[k]),
  93         ftest:ev(ftrunc,root[k]),
  94         gxtest:ev(gtruncx,root[k]),
  95         gytest:ev(gtruncy,root[k]),
  96 /*      print("ftrunc =",ftest),
  97         print("gxtrunc =",gxtest),
  98         print("gytrunc =",gytest),             */
  99         if ftest=0 then (print("ftrunc is zero!"), flag:red)  else
 100         if gxtest=0 and gytest=0 then (print("dg trunc is zero!"),flag:red)),
 101         if flag=green then "passed test" else "failed test"))$
 102
 103 /* subroutine to define f and g */
 104 deffg():=(
 105 f::expand(ev(x*u+y*v)),
 106 print(f,"=",ev(f)),
 107 g::expand(ev(x*v-y*u)),
 108 print(g,"=",ev(g)))$
 109
 110 /* subroutine to define p and q */
 111 defpq():=(
 112 trans:[ev(x)=r*cos(s),ev(y)=r*sin(s)],
 113 p::ev(f)/r,p::expand(ev(ev(p),trans)),
 114 print(p,"=",ev(p)),
 115 q::ev(g)/r^2,q::expand(ev(ev(q),trans)),
 116 print(q,"=",ev(q)))$
 117
 118 /* subroutine to define pp and qq */
 119 defppqq():=(
 120 exponent:min(lopow(ev(p),r),lopow(ev(q),r)),
 121 pp::expand(ev(p)/r^exponent),
 122 qq::expand(ev(q)/r^exponent),
 123 print("divide out",r^exponent),
 124 /*       print(pp,"=",ev(pp)),
 125          print(qq,"=",ev(qq)),           */
 126 print("now set",r,"= 0"),
 127 ptemp:ev(ev(pp),r::0),
 128 qtemp:ev(ev(qq),r::0),
 129 print(pp,"=",ptemp),
 130 print("note: previous should be zero!"),
 131 print(qq,"=",qtemp))$
 132
 133 /* subroutine to find roots s of qq=0 when r=0 */
 134 /*     user selects root sstar to be used      */
 135 sroots():=(
 136 stemp:solve(qtemp,ev(s)),
 137 for k:1 thru length(stemp) do
 138         print("root no.",k,",",part(stemp,k)),
 139         print("there are",length(stemp),"roots"),
 140         rootno:read("pick a root no., or 0 to enter one"),
 141         if rootno=0 then sstar:read("enter root") else
 142         sstar:rhs(part(stemp,rootno)),
 143         print(s,"star =",sstar))$
 144
 145 /* subroutine to taylor expand pp and qq about r=0, s=star */
 146 /*        returns new u and v for next iteration           */
 147
 148 blowup():=(
 149 r::ev(x),
 150 s::sstar+ev(y),
 151 pow:read("keep terms of what power?"),
 152 print(u,"="),
 153 u::taylor(ev(ev(pp)),[ev(x),ev(y)],0,pow),
 154 print(ev(u)),
 155 print(v,"="),
 156 v::taylor(ev(ev(qq)),[ev(x),ev(y)],0,pow),
 157 print(ev(v)) )$