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