1/* Filename averm.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 8: averm(), allows mth order averaging. see page 128 in 17 "perturbation methods, bifurcation theory and computer algebra". */ 18 19 20/* program to perform mth order averaging 21 on an n-dimensional system of nonautonomous ode's */ 22 23/* averaging is performed by converting trig terms to 24 complex exponentials, then killing exponentials */ 25 26averm():=block( 27 choice1:read("do you want to enter a new problem? (y/n)"), 28 if choice1 = n then go(jump1), 29 kill(x), 30 print("are you considering a weakly nonlinear oscillator of the form"), 31 choice2:read("z'' + omega0^2 z = eps f(z,zdot,t) ? (y/n)"), 32 if choice2 = n then go(jump2), 33 34/* enter data for single oscillator problem */ 35 n:2, 36 omega0:read("enter omega0"), 37 f:read("enter f(z,zdot,t)")*eps, 38/* does f(z,zdot,t) depend explicitly on t? */ 39 test:diff(f,t), 40 if test=0 then omega1:omega0 41 else( 42 omega:read("enter the forcing frequency"), 43 k:read("enter the order of the subharmonic resonance"), 44 omega1:omega/k), 45/* van der pol transformation */ 46 zsub:[z=cos(omega1*t)*x1-sin(omega1*t)*x2, 47 zdot=-omega1*sin(omega1*t)*x1-omega1*cos(omega1*t)*x2], 48/* substitute zsub into transformed eqs */ 49/* scale time so that averaging occurs over period 2 pi */ 50 vf:ev([-1/omega1^2*(eps*kapomega/k^2*z + f)*sin(omega1*t), 51 -1/omega1^2*(eps*kapomega/k^2*z + f)*cos(omega1*t)], 52 zsub,t=tau/omega1,infeval), 53 vf:ev(vf,tau=t), 54 if omega1 # omega0 55 then print("we write eps*kapomega =",omega^2-k^2*omega0^2) 56 else vf:ev(vf,kapomega=0), 57 vf2:expand(exponentialize(vf)), 58 for i:1 thru 2 do vf2[i]:map(factor,vf2[i]), 59 x:[x1,x2], 60 go(jump1), 61 62 63jump2, 64/* enter data for general problem of n ode's */ 65 n:read("enter number of differential equations"), 66 x:makelist(concat('x,i),i,1,n), 67 print("the ode's are of the form:"), 68 print("dx/dt = eps f(x,t)"), 69 print("where x =",x), 70 print("scale time t such that averaging occurs over 2 pi"), 71 vf:makelist(read("enter rhs(",i,")=eps*...")*eps,i,1,n), 72 for i:1 thru n do print("d",x[i],"/dt =",vf[i]), 73 vf2:expand(exponentialize(vf)), 74 for i:1 thru n do vf2[i]:map(factor,vf2[i]), 75 76 77jump1, 78/* averaging */ 79 m:read("enter order of averaging"), 80 if m=1 81/* first order averaging */ 82 then (temp0:demoivre(apply1(f,killexp)), 83 result:taylor(temp0,eps,0,1)) 84/* averaging of order m > 1 */ 85 else 86 y:makelist(concat('y,i),i,1,n), 87 wlist:makelist(concat('w,i),i,1,n), 88 depends(wlist,cons(t,y)), 89 trafo:y, 90 xsub:maplist("=",y,x), 91/* wnull sets wlist to zero */ 92 wnull:vfun(wlist,0), 93 jacob[i,j] := diff(wlist[i],y[j]), 94 jac:genmatrix(jacob,n,n), 95/* main loop */ 96 for i :1 thru m-1 do ( 97 temp1:maplist("=",x,y+eps^i*wlist), 98/* here x is the list [x1,x2,...,xn], y is the list [y1,y2,...,yn], 99 wlist is the transformation vector [w1,w2,...,wn] */ 100 temp2:taylor(sum((-eps)^(i*j)*jac^^j,j,0,m-1). 101 (ev(vf2,temp1) - diff(wlist,t)*eps^i),eps,0,m), 102/* jac is the jacobian d wlist/dy of the transformation wlist */ 103 temp3:mattolist(temp2,n), 104 temp4:map(expand,taylor(ev(temp3,wnull,diff),eps,0,i)), 105/* the ith order mean */ 106 mean:apply1(temp4,killexp), 107 temp6:expand(temp4-mean), 108 temp7:ev(integrate(temp6,t),eps=1), 109/* the ith order transformation */ 110 temp8:maplist("=",wlist,temp7), 111 temp9:ratsimp(temp8), 112/* the transformed de */ 113 vf2:expand(ev(temp3,temp9,diff,xsub,infeval)), 114/* the sum of all transformations */ 115 trafo:expand(trafo+ev(eps^i*wlist,temp9))), 116/* end of main loop */ 117 print("the transformation: ",x,"="), 118 print(ratsimp(demoivre(trafo))), 119/* the final averaging */ 120 result:apply1(vf2,killexp), 121/* replace x by y */ 122 for i:1 thru n do result:subst(concat('y,i),concat('x,i),result), 123 result)$ 124 125 126/* auxiliary functions */ 127 vfun(list,value):=map(lambda([u],u=value),list)$ 128 mattolist(mat,dim):=if dim>1 then makelist(mat[i,1],i,1,dim) else [mat]$ 129 contains_t(exp):= not freeof(t,exp)$ 130 matchdeclare(q,contains_t)$ 131 defrule(killexp,exp(q),0)$ 132