1/* Filename twovar.mac
2
3   ***************************************************************
4   *							         *
5   *                     <package name>                          *
6   *                <functionality description>                  *
7   *                                                             *
8   *          from: Perturbation Methods, Bifurcation            *
9   *                Theory and Computer Algebra.                 *
10   *           by Rand & Armbruster (Springer 1987)              *
11   *                Programmed by Richard Rand                   *
12   *      These files are released to the public domain          *
13   *            						 *
14   ***************************************************************
15*/
16/* program number 6: twovar(), a function to perform a two variable
17   expansion to order epsilon for n weakly coupled perturbed nonlinear
18   oscillators. see page 94 in "perturbation methods, bifurcation
19   theory and computer algebra". */
20
21
22
23twovar():=block(
24
25choice:read("do you want to enter new data (y/n)"),
26
27if choice = n then go(jump1),
28
29/* input d.e.'s */
30n:read("number of d.e.'s"),
31
32print("the",n,"d.e.'s will be in the form:"),
33print("x[i]'' + w[i]^2 x[i] = e f[i](x[1],...,x[",n,"],t)"),
34
35for i:1 thru n do
36    x[i]:read("enter symbol for x[",i,"]"),
37
38for i:1 thru n do
39    depends(x[i],t),
40
41for i:1 thru n do
42    w[i]:read("enter w[",i,"]"),
43
44for i:1 thru n do
45    f[i]:read("enter f[",i,"]"),
46
47jump1,
48
49/* update eqs for substitution of resonant values on 2nd time thru */
50for i:1 thru n do(
51   w[i]:ev(w[i]),
52   f[i]:ev(f[i])),
53
54/* echo back the d.e.'s */
55print("the d.e.'s are entered as:"),
56for i:1 thru n do
57    print(x[i],"'' +",w[i]^2*x[i],"=",e*f[i]),
58
59print("the method assumes a solution in the form:"),
60print("x[i] = x0[i] + e x1[i]"),
61print("where x0[i] = a[i](eta) cos w[i] xi + b[i](eta) sin w[i] xi"),
62print("where xi = t and eta = e t"),
63
64/* don't use a or b as parameters in the given d.e.'s */
65depends([a,b],eta),
66
67for i:1 thru n do
68    x0[i]:a[i]*cos(w[i]*xi)+b[i]*sin(w[i]*xi),
69
70for i:1 thru n do
71    g[i]:ev(f[i],t=xi),
72
73for i:1 thru n do(
74    for j:1 thru n do
75        g[i]:ev(g[i],x[j]::x0[j])),
76
77for i:1 thru n do(
78    g[i]:g[i]-2*diff(x0[i],xi,1,eta,1),
79    g[i]:ev(g[i],diff),
80    g[i]:expand(trigreduce(expand(g[i])))),
81
82/* collect secular terms */
83for i:1 thru n do(
84    s[i]:coeff(g[i],sin(w[i]*xi)),
85    c[i]:coeff(g[i],cos(w[i]*xi))),
86
87/* display secular terms */
88print("removal of secular terms in the x1[i] eqs. gives:"),
89
90for i:1 thru n do(
91    print(s[i],"= 0"),
92    print(c[i],"= 0")),
93
94abeqs:append(makelist(s[i],i,1,n),makelist(c[i],i,1,n)),
95
96choice2:read("do you want to transform to polar coordinates (y/n)"),
97
98if choice2=n then go(jump2),
99
100/* transform to polar coordinates */
101depends([r,theta],eta),
102trans:makelist([a[i]=r[i]*cos(theta[i]),b[i]=r[i]*sin(theta[i])],i,1,n),
103inteqs:ev(abeqs,trans,diff),
104rtheqs:solve(inteqs,append(makelist(diff(r[i],eta),i,1,n),
105                     makelist(diff(theta[i],eta),i,1,n))),
106rtheqs:trigsimp(rtheqs),
107rtheqs:expand(trigreduce(rtheqs)),
108print(rtheqs),
109
110jump2,
111
112choice3:read("do you want to search for resonant parameter values (y/n)"),
113
114if choice3=n then go(end),
115
116/* obtain frequencies which appear on rhs's of d.e.'s */
117
118/* define pattern matching rules to isolate freqs */
119matchdeclare([dummy1,dummy2],true),
120defrule(cosarg,dummy1*cos(dummy2),dummy2),
121defrule(sinarg,dummy1*sin(dummy2),dummy2),
122
123for i:1 thru n do(
124/* change sum to a list */
125    temp1:args(g[i]),
126/* remove constant terms */
127    temp2:map(trigidentify,temp1),
128/* isolate arguments of trig terms */
129    temp3:apply1(temp2,cosarg,sinarg),
130    temp4:ev(temp3,xi=1),
131/* remove duplicate arguments */
132    freq[i]:setify(temp4)),
133
134/* display frequencies */
135for i:1 thru n do(
136    print(x[i],"eq's resonant freq =",w[i]),
137    print("freqs on rhs =",freq[i])),
138
139jump3,
140
141par:read("which parameter to search for ?"),
142
143/* compute resonant values of desired parameter */
144resvals:[],
145for i:1 thru n do(
146   for j:1 thru length(freq[i]) do(
147      res:solve(part(freq[i],j)=w[i],par),
148      if (res#[] and res#all) then resvals:append(resvals,res),
149      res:solve(part(freq[i],j)=-w[i],par),
150      if (res#[] and res#all) then resvals:append(resvals,res))),
151
152/* eliminate duplicate parameter values */
153resvalues:setify(resvals),
154
155/* display resonant parameter values */
156print(resvalues),
157
158choice4:read("do you want to search for another parameter (y/n) ?"),
159
160if choice4=y then go(jump3),
161
162end," ")$
163
164/* auxiliary functions */
165
166trigidentify(exp):=if freeof(sin,exp) and freeof(cos,exp) then 0 else exp$
167
168setify(list):=(
169    set:[list[1]],
170    for i:2 thru length(list) do(
171        if not member(list[i],set) then set:cons(list[i],set)),
172    set)$
173