1/* Filename lc.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 1: lc(), lindstedt's method to calculate limit cycles,
17  see page 7 in "perturbation methods, bifurcation theory  and computer
18  algebra".*/
19
20/* this program applies lindstedt's method to
21the equation:
22               x'' + x + e f(x,x') = 0,
23assuming a limit cycle exists.
24call it with:
25                   lc();
26*/
27
28lc():=(
29
30/* input the differential equation */
31kill(x,xlist,paramlist),
32print("the d.e. is of the form:  x'' + x + e * f(x,x') = 0"),
33f:read("enter f(x,y), representing x' as y"),
34print("the d.e. is: x'' + x + e (",f,") = 0"),
35f:subst('diff(x,z,1)*w,y,f),
36
37/* set up the series expansions */
38n:read("enter truncation order"),
39w:1,
40for i thru n do w:w+k[i]*e^i,
41x:b[0]*cos(z),
42xlist:[xx[0] = x],
43for i thru n do x:x+xx[i]*e^i,
44
45/* plug into the d.e. and collect terms */
46depends(xx,z),
47temp1:diff(x,z,2)+x/w^2+e*ev(f,diff)/w^2,
48temp1:taylor(temp1,e,0,n),
49for i thru n do eq[i]:coeff(temp1,e,i),
50
51/* set up pattern matching rules for use in solving d.e. */
52matchdeclare(n1,true),
53defrule(c,cos(n1*z),cos(n1*z)/(n1*n1-1)),
54defrule(s,sin(n1*z),sin(n1*z)/(n1*n1-1)),
55
56/* load poisson series package and set parameter */
57outofpois(dummy),
58poislim:100,
59
60/* main loop */
61for i:1 thru n do block(
62
63/* trigonometric simplification */
64temp1:outofpois(ev(eq[i],xlist,paramlist,diff)),
65
66/* eliminate secular terms */
67if i = 1
68     then (paramlist:solve(coeff(temp1,sin(z)),b[0]),
69           print("choices for limit cycle amplitude:"),
70	   for j:1 thru length(paramlist) do
71                print(j,")  ",part(paramlist,j,2)),
72           r1:read("enter choice number"),
73	   paramlist:append(solve(coeff(temp1,cos(z)),k[1]),
74		[part(paramlist,r1)]))
75     else  paramlist:append(paramlist,solve([coeff(temp1,cos(z)),
76		coeff(temp1,sin(z))],[k[i],b[i-1]])),
77temp1:expand(ev(temp1,paramlist)),
78xlist:expand(ev(xlist,paramlist)),
79
80/* output progress */
81print("done with step",i,"of",n),
82
83/* exit here if last iteration */
84if i=n then go(end),
85
86/* solve the d.e. */
87temp1:factor(ev(temp1,xx[i] = 0)),
88temp1:applyb1(temp1,c,s),
89temp1:xx[i] = temp1+a[i]*sin(z)+b[i]*cos(z),
90
91/* fit the initial condition */
92temp2:rhs(temp1),
93temp2:diff(temp2,z),
94temp2:solve(ev(temp2,z:0),a[i]),
95xlist:append(xlist,[ev(temp1,temp2)]),
96
97/* end of main loop */
98end),
99
100/* output results */
101w:ev(w,paramlist),
102x:taylor(ev(x,xlist,paramlist),e,0,n-1),
103print("x =",x),
104print("w =",w))$
105