1/* Filename average.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 7: average(), allows one to perform first and second order 17 averaging. see page 114 in "perturbation methods, bifurcation theory 18 and computer algebra". */ 19 20 21 22/* program to perform 1st or 2nd order averaging 23 on an n-dimensional system of nonautonomous ode's */ 24 25/* averaging is performed by converting trig terms to 26 complex exponentials, then killing exponentials */ 27 28average():=block( 29 choice1:read("do you want to enter a new problem? (y/n)"), 30 if choice1 = n then go(jump1), 31 kill(x), 32 print("are you considering a weakly nonlinear oscillator of the form"), 33 choice2:read("z'' + omega0^2 z = eps f(z,zdot,t) ? (y/n)"), 34 if choice2 = n then go(jump2), 35 36/* enter data for single oscillator problem */ 37 n:2, 38 omega0:read("enter omega0"), 39 f:read("enter f(z,zdot,t)")*eps, 40/* does f(z,zdot,t) depend explicitly on t? */ 41 test:diff(f,t), 42 if test=0 then omega1:omega0 43 else( 44 omega:read("enter the forcing frequency"), 45 k:read("enter the order of the subharmonic resonance"), 46 omega1:omega/k), 47/* van der pol transformation */ 48 zsub:[z=cos(omega1*t)*x1-sin(omega1*t)*x2, 49 zdot=-omega1*sin(omega1*t)*x1-omega1*cos(omega1*t)*x2], 50/* substitute zsub into transformed eqs */ 51/* scale time so that averaging occurs over period 2 pi */ 52 vf:ev([-1/omega1^2*(eps*kapomega/k^2*z + f)*sin(omega1*t), 53 -1/omega1^2*(eps*kapomega/k^2*z + f)*cos(omega1*t)], 54 zsub,t=tau/omega1,infeval), 55 vf:ev(vf,tau=t), 56 if omega1 # omega0 57 then print("we write eps*kapomega =",omega^2-k^2*omega0^2) 58 else vf:ev(vf,kapomega=0), 59 vf2:expand(exponentialize(vf)), 60 for i:1 thru 2 do vf2[i]:map(factor,vf2[i]), 61 x:[x1,x2], 62 go(jump1), 63 64 65jump2, 66/* enter data for general problem of n ode's */ 67 n:read("enter number of differential equations"), 68 x:makelist(concat('x,i),i,1,n), 69 print("the ode's are of the form:"), 70 print("dx/dt = eps f(x,t)"), 71 print("where x =",x), 72 print("scale time t such that averaging occurs over 2 pi"), 73 vf:makelist(read("enter rhs(",i,")=eps*...")*eps,i,1,n), 74 for i:1 thru n do print("d",x[i],"/dt =",vf[i]), 75 vf2:expand(exponentialize(vf)), 76 for i:1 thru n do vf2[i]:map(factor,vf2[i]), 77 78 79jump1, 80/* averaging */ 81 m:read("do you want first or second order averaging? (1/2)"), 82 coefvfeps1:coeff(vf2,eps), 83 coefvfeps2:coeff(vf2,eps,2), 84 g1bar:demoivre(apply1(coefvfeps1,killexp)), 85 if m=1 then result:eps*g1bar 86 else( 87 g1tilde:coefvfeps1-g1bar, 88 w1:integrate(g1tilde,t), 89 remarray(jacob), 90 jacob[i,j] := diff(g1tilde[i],x[j]), 91 jacg1tilde:genmatrix(jacob,n,n), 92 g2bar:demoivre(apply1(expand(jacg1tilde.w1)+coefvfeps2,killexp)), 93 result:eps*g1bar+eps^2*g2bar), 94/* replace x by y */ 95 for i:1 thru n do result:subst(concat('y,i),concat('x,i),result), 96 result)$ 97 98 99/* auxiliary functions to kill exponential terms */ 100 contains_t(exp):= not freeof(t,exp)$ 101 matchdeclare(q,contains_t)$ 102 defrule(killexp,exp(q),0)$ 103