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