1/* Filename averm.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 8: averm(), allows mth order averaging. see page 128 in
17   "perturbation methods, bifurcation theory and computer algebra". */
18
19
20/* program to perform mth order averaging
21   on an n-dimensional system of nonautonomous ode's */
22
23/* averaging is performed by converting trig terms to
24   complex exponentials, then killing exponentials   */
25
26averm():=block(
27      choice1:read("do you want to enter a new problem? (y/n)"),
28      if choice1 = n then go(jump1),
29      kill(x),
30      print("are you considering a weakly nonlinear oscillator of the form"),
31      choice2:read("z'' + omega0^2 z = eps f(z,zdot,t) ? (y/n)"),
32      if choice2 = n then go(jump2),
33
34/* enter data for single oscillator problem */
35      n:2,
36      omega0:read("enter omega0"),
37      f:read("enter f(z,zdot,t)")*eps,
38/* does f(z,zdot,t) depend explicitly on t? */
39      test:diff(f,t),
40      if test=0 then omega1:omega0
41            else(
42      	    omega:read("enter the forcing frequency"),
43            k:read("enter the order of the subharmonic resonance"),
44            omega1:omega/k),
45/* van der pol transformation */
46      zsub:[z=cos(omega1*t)*x1-sin(omega1*t)*x2,
47               zdot=-omega1*sin(omega1*t)*x1-omega1*cos(omega1*t)*x2],
48/* substitute zsub into transformed eqs */
49/* scale time so that averaging occurs over period 2 pi */
50      vf:ev([-1/omega1^2*(eps*kapomega/k^2*z + f)*sin(omega1*t),
51             -1/omega1^2*(eps*kapomega/k^2*z + f)*cos(omega1*t)],
52             zsub,t=tau/omega1,infeval),
53      vf:ev(vf,tau=t),
54      if omega1 # omega0
55      then print("we write eps*kapomega =",omega^2-k^2*omega0^2)
56            else vf:ev(vf,kapomega=0),
57      vf2:expand(exponentialize(vf)),
58      for i:1 thru 2 do vf2[i]:map(factor,vf2[i]),
59      x:[x1,x2],
60      go(jump1),
61
62
63jump2,
64/* enter data for general problem of n ode's */
65      n:read("enter number of differential equations"),
66      x:makelist(concat('x,i),i,1,n),
67      print("the ode's are of the form:"),
68      print("dx/dt = eps f(x,t)"),
69      print("where x =",x),
70      print("scale time t such that averaging occurs over 2 pi"),
71      vf:makelist(read("enter rhs(",i,")=eps*...")*eps,i,1,n),
72      for i:1 thru n do print("d",x[i],"/dt =",vf[i]),
73      vf2:expand(exponentialize(vf)),
74      for i:1 thru n do vf2[i]:map(factor,vf2[i]),
75
76
77jump1,
78/* averaging */
79      m:read("enter order of averaging"),
80      if m=1
81/* first order averaging */
82      then (temp0:demoivre(apply1(f,killexp)),
83           result:taylor(temp0,eps,0,1))
84/* averaging of order m > 1 */
85      else
86      y:makelist(concat('y,i),i,1,n),
87      wlist:makelist(concat('w,i),i,1,n),
88      depends(wlist,cons(t,y)),
89      trafo:y,
90      xsub:maplist("=",y,x),
91/* wnull sets wlist to zero */
92      wnull:vfun(wlist,0),
93      jacob[i,j] := diff(wlist[i],y[j]),
94      jac:genmatrix(jacob,n,n),
95/* main loop */
96      for i :1 thru m-1 do (
97	  temp1:maplist("=",x,y+eps^i*wlist),
98/* here x is the list [x1,x2,...,xn], y is the list [y1,y2,...,yn],
99        wlist is the transformation vector [w1,w2,...,wn] */
100	  temp2:taylor(sum((-eps)^(i*j)*jac^^j,j,0,m-1).
101		                (ev(vf2,temp1) - diff(wlist,t)*eps^i),eps,0,m),
102/* jac is the jacobian d wlist/dy of the transformation wlist */
103          temp3:mattolist(temp2,n),
104          temp4:map(expand,taylor(ev(temp3,wnull,diff),eps,0,i)),
105/* the ith order mean */
106	  mean:apply1(temp4,killexp),
107	  temp6:expand(temp4-mean),
108	  temp7:ev(integrate(temp6,t),eps=1),
109/* the ith order transformation */
110          temp8:maplist("=",wlist,temp7),
111	  temp9:ratsimp(temp8),
112/* the transformed de */
113	  vf2:expand(ev(temp3,temp9,diff,xsub,infeval)),
114/* the sum of all transformations */
115	  trafo:expand(trafo+ev(eps^i*wlist,temp9))),
116/* end of main loop */
117      print("the transformation: ",x,"="),
118      print(ratsimp(demoivre(trafo))),
119/* the final averaging */
120      result:apply1(vf2,killexp),
121/* replace x by y */
122      for i:1 thru n do result:subst(concat('y,i),concat('x,i),result),
123      result)$
124
125
126/* auxiliary functions */
127      vfun(list,value):=map(lambda([u],u=value),list)$
128      mattolist(mat,dim):=if dim>1 then makelist(mat[i,1],i,1,dim) else [mat]$
129      contains_t(exp):= not freeof(t,exp)$
130      matchdeclare(q,contains_t)$
131      defrule(killexp,exp(q),0)$
132