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