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 */ 20takens():=block( 21/* variable i tags loop */ 22for 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(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 */ 39setup1():=( 40u:concat('u,i), 41v:concat('v,i), 42x:concat('x,i), 43y:concat('y,i))$ 44 45setup2():=( 46f:concat('f,i), 47g:concat('g,i), 48p:concat('p,i), 49q:concat('q,i), 50r:concat('r,i), 51s:concat('s,i), 52pp:concat('pp,i), 53qq:concat('qq,i))$ 54 55/* subroutine to input the rhs's from keyboard */ 56inputrhs():=( 57print(" enter the rhs's to be studied "), 58print(" use variables x,y, they will be converted to x1,y1 "), 59u::read(u,"="), 60print(ev(u)), 61v::read(v,"="), 62print(ev(v)))$ 63 64/* subroutine to truncate f and g to terms of lowest degree */ 65truncate():=block( 66for 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 */ 72getroots():=( 73print("solving gtrunc = 0"), 74ftrunc:part(temp2,1), 75gtrunc:part(temp2,2), 76gtruncx:diff(gtrunc,ev(x)), 77gtruncy:diff(gtrunc,ev(y)), 78xroots:solve(gtrunc,ev(x)), 79yroots:solve(gtrunc,ev(y)), 80rootnum:0, 81for k:1 thru length(xroots) do 82 (rootnum:rootnum+1, 83 root[rootnum]:part(xroots,k)), 84for k:1 thru length(yroots) do 85 (rootnum:rootnum+1, 86 root[rootnum]:part(yroots,k)), 87print("total no. of roots =",rootnum))$ 88 89/* perform takens' test for each root */ 90test():=(block( 91flag:green, 92for 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 */ 104deffg():=( 105f::expand(ev(x*u+y*v)), 106print(f,"=",ev(f)), 107g::expand(ev(x*v-y*u)), 108print(g,"=",ev(g)))$ 109 110/* subroutine to define p and q */ 111defpq():=( 112trans:[ev(x)=r*cos(s),ev(y)=r*sin(s)], 113p::ev(f)/r,p::expand(ev(ev(p),trans)), 114print(p,"=",ev(p)), 115q::ev(g)/r^2,q::expand(ev(ev(q),trans)), 116print(q,"=",ev(q)))$ 117 118/* subroutine to define pp and qq */ 119defppqq():=( 120exponent:min(lopow(ev(p),r),lopow(ev(q),r)), 121pp::expand(ev(p)/r^exponent), 122qq::expand(ev(q)/r^exponent), 123print("divide out",r^exponent), 124/* print(pp,"=",ev(pp)), 125 print(qq,"=",ev(qq)), */ 126print("now set",r,"= 0"), 127ptemp:ev(ev(pp),r::0), 128qtemp:ev(ev(qq),r::0), 129print(pp,"=",ptemp), 130print("note: previous should be zero!"), 131print(qq,"=",qtemp))$ 132 133/* subroutine to find roots s of qq=0 when r=0 */ 134/* user selects root sstar to be used */ 135sroots():=( 136stemp:solve(qtemp,ev(s)), 137for 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 148blowup():=( 149r::ev(x), 150s::sstar+ev(y), 151pow:read("keep terms of what power?"), 152print(u,"="), 153u::taylor(ev(ev(pp)),[ev(x),ev(y)],0,pow), 154print(ev(u)), 155print(v,"="), 156v::taylor(ev(ev(qq)),[ev(x),ev(y)],0,pow), 157print(ev(v)) )$ 158