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