1/* Filename lc2.mac
2
3   ***************************************************************
4   *							         *
5   *                     <package name>                          *
6   *                <functionality description>                  *
7   *                                                             *
8   *          from: Perturbation Methods, Bifurcation            *
9   *                Theory and Computer Algebra.                 *
10   *           by Rand & Armbruster (Springer 1987)              *
11   *                Programmed by Richard Rand                   *
12   *      These files are released to the public domain          *
13   *            						 *
14   ***************************************************************
15*/
16/* program number 2: lc2(), lindstedt's method for more general
17   differential equations. see page 14 in "perturbation methods,
18   bifurcation theory and computer algebra". */
19
20
21
22
23lc2():=(
24
25/* input the differential equation */
26kill(x,y,xylist,paramlist),
27print("the d.e.'s are of the form x' = -y + e*f, y' = x + e*g"),
28print("where f,g are of the form:
29          quadratic + e*(linear + cubic) + e^2*(quartic) + ..."),
30f:read("enter f:"),
31print("f =",f),
32g:read("enter g:"),
33print("g =",g),
34
35/* set up the series expansions */
36n:read("enter truncation order:"),
37k[0]:1,
38k[1]:0,
39w:sum(k[i]*e^i,i,0,n),
40xy:[x:b[0]*cos(z)+sum(xx[i](z)*e^i,i,1,n),
41    y:b[0]*sin(z)+sum(yy[i](z)*e^i,i,1,n)],
42xylist:[xx[0](z)=b[0]*cos(z),
43        yy[0](z)=b[0]*sin(z)],
44
45/* plug into d.e.'s and collect terms */
46temp1:[-diff(x,z)*w-y+e*ev(f,diff),-diff(y,z)*w+x+e*ev(g,diff)],
47temp2:taylor(temp1,e,0,n),
48for i:1 thru n do eq[i]:coeff(temp2,e,i),
49
50
51/* main loop */
52for i:1 thru n do block(
53
54/* trigonometric simplification */
55temp3:expand(trigreduce(expand(ev(eq[i],xylist,paramlist,diff)))),
56
57/* eliminate secular terms */
58if i=1 then (temp4:temp3, go(skip_to_here_first_time))
59       else newparamlist:
60             solve([coeff(part(temp3,1),sin(z))-coeff(part(temp3,2),cos(z)),
61                    coeff(part(temp3,1),cos(z))+coeff(part(temp3,2),sin(z))],
62                   [b[i-2],k[i]]),
63if i=2 then (paramlist:newparamlist,
64             print("choices for limit cycle amplitude:"),
65             for j:1 thru length(paramlist) do
66                 print(j,")  ",part(paramlist,j,1,2)),
67             r1:read("enter choice number"),
68             paramlist:part(paramlist,r1))
69        else paramlist:append(paramlist,newparamlist),
70temp4:expand(ev(temp3,paramlist)),
71xylist:expand(ev(xylist,paramlist)),
72skip_to_here_first_time,
73
74/* output progress */
75print("done with step",i,"of",n),
76
77/* exit here if last iteration */
78if i=n then go(end),
79
80/* solve the d.e.'s */
81temp4a:subst(dummy(z),yy[i](z),temp4),
82atvalue(dummy(z),z=0,0),
83temp5:desolve(temp4a,[xx[i](z),dummy(z)]),
84temp5a:subst(yy[i](z),dummy(z),temp5),
85temp5b:subst(b[i],xx[i](0),temp5a),
86xylist:append(xylist,[temp5b]),
87
88/* end of main loop */
89end),
90
91/* output results */
92w:ev(w,paramlist),
93soln:taylor(ev([x,y],xylist,paramlist),e,0,n-2),
94print("x =",part(soln,1)),
95print("y =",part(soln,2)),
96print("w =",w))$
97