share/contrib/rand/Avg3.mac

   1  /* Mathieu's equation: third order averaging */
   2
   3  /* setup */
   4  depends([a,psi],t);
   5  depends([w1,w2,v1,v2,u1,u2],[abar,psibar,t]);
   6  depends([abar,psibar],t);
   7  x:a*cos(t+psi);
   8  xd:-a*sin(t+psi);
   9  f:-4*x*(del1+e*del2+e^2*del3+cos(2*t));
  10  de1:[diff(a,t)=-e*sin(t+psi)*f,diff(psi,t)=-e/a*cos(t+psi)*f];
  11  transf:[a=abar+e*w1+e^2*v1+e^3*u1,psi=psibar+e*w2+e^2*v2+e^3*u2];
  12  de2:de1,transf,diff$
  13  de3:solve(%,[diff(abar,t),diff(psibar,t)])$
  14  de4:taylor(de3,e,0,3)$
  15  for i:0 thru 3 do eq[i]:coeff(de4,e,i);
  16  neweq[0]:eq[0];
  17  kill(abar,psibar);
  18
  19  /* first order averaging */
  20  ws1:eq[1]$
  21  ws2:expand(trigreduce(expand(ws1)))$
  22  ws3:integrate(ws2,t)$
  23  ws4:ws3,%integconst1:0,%integconst2:0$
  24  ws5:ws4-ratcoef(ws4,t)*t$
  25  ws6:solve(part(ws5,1),[w1,w2]);
  26  ws7:ws1,ws6,diff;
  27  neweq[1]:expand(trigreduce(expand(ws7)));
  28
  29  /* second order averaging */
  30  vs1:eq[2],ws6,diff$
  31  vs2:expand(trigreduce(expand(vs1)))$
  32  vs3:integrate(vs2,t)$
  33  vs4:vs3,%integconst3:0,%integconst4:0$
  34  vs5:vs4-ratcoef(vs4,t)*t$
  35  vs6:solve(part(vs5,1),[v1,v2]);
  36  vs7:vs1,vs6,diff$
  37  neweq[2]:expand(trigreduce(expand(vs7)));
  38
  39  /* third order averaging */
  40  us1:eq[3],ws6,vs6,diff$
  41  us2:expand(trigreduce(expand(us1)))$
  42  us3:integrate(us2,t)$
  43  us4:us3,%integconst5:0,%integconst6:0$
  44  us5:us4-ratcoef(us4,t)*t$
  45  us6:solve(part(us5,1),[u1,u2]);
  46  us7:us1,us6,diff$
  47  neweq[3]:expand(trigreduce(expand(us7)));
  48
  49  neweqs:sum(neweq[i]*e^i,i,0,3);
  50
  51  depends([xbar,ybar],t);
  52  depends([abar,psibar],t);
  53  [xbar = abar*cos(psibar),ybar = -abar*sin(psibar)];
  54  diff(%,t);
  55  ev(%,neweqs);
  56  trigexpand(%);
  57  expand(%);
  58  trigsimp(%);
  59  taylor(%,e,0,3);
  60  ev(%,sin(psibar) = -ybar/abar,cos(psibar) = xbar/abar);
  61  factor(%);
  62