1/* Filename average.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 7: average(), allows one to perform first and second order
17  averaging. see page 114 in "perturbation methods, bifurcation theory
18  and computer algebra". */
19
20
21
22/* program to perform 1st or 2nd order averaging
23   on an n-dimensional system of nonautonomous ode's */
24
25/* averaging is performed by converting trig terms to
26   complex exponentials, then killing exponentials   */
27
28average():=block(
29      choice1:read("do you want to enter a new problem? (y/n)"),
30      if choice1 = n then go(jump1),
31      kill(x),
32      print("are you considering a weakly nonlinear oscillator of the form"),
33      choice2:read("z'' + omega0^2 z = eps f(z,zdot,t) ? (y/n)"),
34      if choice2 = n then go(jump2),
35
36/* enter data for single oscillator problem */
37      n:2,
38      omega0:read("enter omega0"),
39      f:read("enter f(z,zdot,t)")*eps,
40/* does f(z,zdot,t) depend explicitly on t? */
41      test:diff(f,t),
42      if test=0 then omega1:omega0
43            else(
44      	    omega:read("enter the forcing frequency"),
45            k:read("enter the order of the subharmonic resonance"),
46            omega1:omega/k),
47/* van der pol transformation */
48      zsub:[z=cos(omega1*t)*x1-sin(omega1*t)*x2,
49               zdot=-omega1*sin(omega1*t)*x1-omega1*cos(omega1*t)*x2],
50/* substitute zsub into transformed eqs */
51/* scale time so that averaging occurs over period 2 pi */
52      vf:ev([-1/omega1^2*(eps*kapomega/k^2*z + f)*sin(omega1*t),
53             -1/omega1^2*(eps*kapomega/k^2*z + f)*cos(omega1*t)],
54             zsub,t=tau/omega1,infeval),
55      vf:ev(vf,tau=t),
56      if omega1 # omega0
57      then print("we write eps*kapomega =",omega^2-k^2*omega0^2)
58            else vf:ev(vf,kapomega=0),
59      vf2:expand(exponentialize(vf)),
60      for i:1 thru 2 do vf2[i]:map(factor,vf2[i]),
61      x:[x1,x2],
62      go(jump1),
63
64
65jump2,
66/* enter data for general problem of n ode's */
67      n:read("enter number of differential equations"),
68      x:makelist(concat('x,i),i,1,n),
69      print("the ode's are of the form:"),
70      print("dx/dt = eps f(x,t)"),
71      print("where x =",x),
72      print("scale time t such that averaging occurs over 2 pi"),
73      vf:makelist(read("enter rhs(",i,")=eps*...")*eps,i,1,n),
74      for i:1 thru n do print("d",x[i],"/dt =",vf[i]),
75      vf2:expand(exponentialize(vf)),
76      for i:1 thru n do vf2[i]:map(factor,vf2[i]),
77
78
79jump1,
80/* averaging */
81      m:read("do you want first or second order averaging?  (1/2)"),
82      coefvfeps1:coeff(vf2,eps),
83      coefvfeps2:coeff(vf2,eps,2),
84      g1bar:demoivre(apply1(coefvfeps1,killexp)),
85      if m=1 then result:eps*g1bar
86                else(
87                g1tilde:coefvfeps1-g1bar,
88                w1:integrate(g1tilde,t),
89      remarray(jacob),
90      jacob[i,j] := diff(g1tilde[i],x[j]),
91      jacg1tilde:genmatrix(jacob,n,n),
92      g2bar:demoivre(apply1(expand(jacg1tilde.w1)+coefvfeps2,killexp)),
93      result:eps*g1bar+eps^2*g2bar),
94/* replace x by y */
95      for i:1 thru n do result:subst(concat('y,i),concat('x,i),result),
96      result)$
97
98
99/* auxiliary functions to kill exponential terms */
100      contains_t(exp):= not freeof(t,exp)$
101      matchdeclare(q,contains_t)$
102      defrule(killexp,exp(q),0)$
103