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