share/contrib/rand/vandpol.mac

   1 /* Filename vandpol.mac
   2
   3    ***************************************************************
   4    *                                                             *
   5    *                     <package name>                          *
   6    *                <functionality description>                  *
   7    *                                                             *
   8    *          from: Computer Algebra in Applied Math.            *
   9    *                   by Rand (Pitman,1984)                     *
  10    *                Programmed by Richard Rand                   *
  11    *      These files are released to the public domain          *
  12    *                                                             *
  13    ***************************************************************
  14 */
  15
  16
  17 /*
  18 (d12) This program computes a perturbation solution for the
  19
  20
  21 limit cycle in Van der Pol's equation.  Call it by typing:
  22
  23
  24                      VANDERPOL(N)
  25
  26
  27 where N is the order of truncation.
  28 */
  29
  30 vanderpol(n):=(setup1(n),setup2(n),
  31           for i thru n do
  32               block(step1(i),step2(i),if i > 1 then output1(i),
  33                     if i = n then go(end),step3(i),step4(i),step5(i),end),
  34           output2(n))$
  35 setup1(n):=(w:1,for i thru n do w:w+k[i]*e^i,x:2*cos(t),
  36        for i thru n do x:x+y[i](t)*e^i)$
  37 setup2(n):=(temp1:diff(x,t,2)+x/w^2-e*(1-x^2)*diff(x,t)/w,
  38        temp1:taylor(temp1,e,0,n),for i thru n do eq[i]:coeff(temp1,e,i))$
  39 step1(i):=temp1
  40       :expand(trigreduce(expand(ev(eq[i],makelist([e[j],f[j]],j,1,i-1),
  41                                    diff))))$
  42 step2(i):=(
  43       if i = 1 then f[i]:solve(coeff(temp1,cos(t)),k[i])
  44           else f[i]:solve([coeff(temp1,cos(t)),coeff(temp1,sin(t))],
  45                           [k[i],b[i-1]]),temp1:ev(temp1,f[i]))$
  46 step3(i):=(temp1:ode2(temp1,y[i](t),t),
  47   temp1:subst(a[i],%k1,temp1),temp1:subst(b[i],%k2,temp1))$
  48 step4(i):=(temp2:rhs(temp1),temp2:diff(temp2,t),
  49       temp2:solve(ev(temp2,t:0),a[i]))$
  50 step5(i):=e[i]:ev(temp1,temp2)$
  51 output1(i):=(print(expand(ev(e[i-1],f[i]))),print(" "))$
  52 output2(n):=(print("w=",ev(w,makelist([f[j]],j,1,n))),print(" "))$