1/* Filename duffing.mac 2 3 *************************************************************** 4 * * 5 * <package name> * 6 * <functionality description> * 7 * * 8 * from: Computer Algebra in Applied Math. * 9 * by Rand (Pitman,1984) * 10 * Programmed by Richard Rand * 11 * These files are released to the public domain * 12 * * 13 *************************************************************** 14*/ /* 15 16(d4) This program computes a perturbation solution for 17 18 19the periodic response of Duffing's equation. 20 21 22Call it by typing: 23 24 25 DUFFING(N) 26 27 28where N is the order of truncation. 29*/ 30 31duffing(n):=(setup1(n),setup2(n), 32 for i thru n do 33 (step1(i),step2(i),step3(i),step4(i),step5(i),output1(i)), 34 output2(n))$ 35setup1(n):=(w:1,for i thru n do w:w+k[i]*e^i,x:a*cos(t), 36 for i thru n do x:x+y[i](t)*e^i)$ 37setup2(n):=(temp1:diff(x,t,2)+x/w^2+e*x^3/w^2-e*f*cos(t)/w^2, 38 temp1:taylor(temp1,e,0,n),for i thru n do eq[i]:coeff(temp1,e,i))$ 39step1(i):=temp1 40 :expand(trigreduce(expand(ev(eq[i],makelist([e[j],f[j]],j,1,i-1), 41 diff))))$ 42step2(i):=(f[i]:solve(coeff(temp1,cos(t)),k[i]),temp1:ev(temp1,f[i]))$ 43step3(i):=(temp1:ode2(temp1,y[i](t),t), 44 temp1:subst(a[i],%k1,temp1),temp1:subst(b[i],%k2,temp1))$ 45step4(i):=(temp2:rhs(temp1),temp2:diff(temp2,t), 46 temp2:solve([ev(rhs(temp1),t:0),ev(temp2,t:0)],[a[i],b[i]]))$ 47step5(i):=e[i]:ev(temp1,temp2)$ 48output1(i):=(print(expand(e[i])),print(" "))$ 49output2(n):=(print("w=",ev(w,makelist([f[j]],j,1,n))),print(" "))$ 50