share/contrib/rand/average.mac

   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
  28 average():=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
  65 jump2,
  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
  79 jump1,
  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)$