share/contrib/rand/averm.mac

   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
  26 averm():=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
  63 jump2,
  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
  77 jump1,
  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)$