share/contrib/rand/twovar.mac

   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
  23 twovar():=block(
  24
  25 choice:read("do you want to enter new data (y/n)"),
  26
  27 if choice = n then go(jump1),
  28
  29 /* input d.e.'s */
  30 n:read("number of d.e.'s"),
  31
  32 print("the",n,"d.e.'s will be in the form:"),
  33 print("x[i]'' + w[i]^2 x[i] = e f[i](x[1],...,x[",n,"],t)"),
  34
  35 for i:1 thru n do
  36     x[i]:read("enter symbol for x[",i,"]"),
  37
  38 for i:1 thru n do
  39     depends(x[i],t),
  40
  41 for i:1 thru n do
  42     w[i]:read("enter w[",i,"]"),
  43
  44 for i:1 thru n do
  45     f[i]:read("enter f[",i,"]"),
  46
  47 jump1,
  48
  49 /* update eqs for substitution of resonant values on 2nd time thru */
  50 for 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 */
  55 print("the d.e.'s are entered as:"),
  56 for i:1 thru n do
  57     print(x[i],"'' +",w[i]^2*x[i],"=",e*f[i]),
  58
  59 print("the method assumes a solution in the form:"),
  60 print("x[i] = x0[i] + e x1[i]"),
  61 print("where x0[i] = a[i](eta) cos w[i] xi + b[i](eta) sin w[i] xi"),
  62 print("where xi = t and eta = e t"),
  63
  64 /* don't use a or b as parameters in the given d.e.'s */
  65 depends([a,b],eta),
  66
  67 for i:1 thru n do
  68     x0[i]:a[i]*cos(w[i]*xi)+b[i]*sin(w[i]*xi),
  69
  70 for i:1 thru n do
  71     g[i]:ev(f[i],t=xi),
  72
  73 for i:1 thru n do(
  74     for j:1 thru n do
  75         g[i]:ev(g[i],x[j]::x0[j])),
  76
  77 for 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 */
  83 for 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 */
  88 print("removal of secular terms in the x1[i] eqs. gives:"),
  89
  90 for i:1 thru n do(
  91     print(s[i],"= 0"),
  92     print(c[i],"= 0")),
  93
  94 abeqs:append(makelist(s[i],i,1,n),makelist(c[i],i,1,n)),
  95
  96 choice2:read("do you want to transform to polar coordinates (y/n)"),
  97
  98 if choice2=n then go(jump2),
  99
 100 /* transform to polar coordinates */
 101 depends([r,theta],eta),
 102 trans:makelist([a[i]=r[i]*cos(theta[i]),b[i]=r[i]*sin(theta[i])],i,1,n),
 103 inteqs:ev(abeqs,trans,diff),
 104 rtheqs:solve(inteqs,append(makelist(diff(r[i],eta),i,1,n),
 105                      makelist(diff(theta[i],eta),i,1,n))),
 106 rtheqs:trigsimp(rtheqs),
 107 rtheqs:expand(trigreduce(rtheqs)),
 108 print(rtheqs),
 109
 110 jump2,
 111
 112 choice3:read("do you want to search for resonant parameter values (y/n)"),
 113
 114 if 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 */
 119 matchdeclare([dummy1,dummy2],true),
 120 defrule(cosarg,dummy1*cos(dummy2),dummy2),
 121 defrule(sinarg,dummy1*sin(dummy2),dummy2),
 122
 123 for 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 */
 135 for i:1 thru n do(
 136     print(x[i],"eq's resonant freq =",w[i]),
 137     print("freqs on rhs =",freq[i])),
 138
 139 jump3,
 140
 141 par:read("which parameter to search for ?"),
 142
 143 /* compute resonant values of desired parameter */
 144 resvals:[],
 145 for 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 */
 153 resvalues:setify(resvals),
 154
 155 /* display resonant parameter values */
 156 print(resvalues),
 157
 158 choice4:read("do you want to search for another parameter (y/n) ?"),
 159
 160 if choice4=y then go(jump3),
 161
 162 end," ")$
 163
 164 /* auxiliary functions */
 165
 166 trigidentify(exp):=if freeof(sin,exp) and freeof(cos,exp) then 0 else exp$
 167
 168 setify(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)$