share/contrib/rand/mathieu0.mac

   1 /* Filename mathieu0.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 (d4) This program computes the transition curve through
  16
  17
  18 the origin (N = 0) in Mathieu's equation using a perturbation
  19
  20
  21 method.  To call it, type
  22
  23
  24                      MATHIEU0()
  25
  26 */
  27
  28 mathieu0():=(input(),setup1(),setup2(),
  29          for i thru m do (step1(),step2a(),step3a()),output())$
  30 input():=m:read("ENTER DEGREE OF TRUNCATION")$
  31 setup1():=(w:0,for i thru m do w:w+k[i]*e^i,x:1,
  32        for i thru m do x:x+y[i](t)*e^i)$
  33 setup2():=(temp1:diff(x,t,2)+x*(w+e*cos(t)),temp1:taylor(temp1,e,0,m),
  34        for i thru m do eq[i]:coeff(temp1,e,i))$
  35 step1():=temp1
  36       :expand(trigreduce(expand(ev(eq[i],makelist([e[j],f[j]],j,1,i-1),
  37                                    diff))))$
  38 step2a():=(temp1:ode2(temp1,y[i](t),t),temp1:ev(temp1,%k1:0,%k2:0))$
  39 step3a():=(f[i]:solve(coeff(expand(rhs(temp1)),t,2),k[i]),e[i]:ev(temp1,f[i]))$
  40 output():=(print("delta=",ev(w,makelist([f[j]],j,1,m))),print(" "))$