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
63