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