share/contrib/rand/reduct2.mac

   1 /* Filename reduct2.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 11: reduction2(), a liapunov-schmidt reduction for
  17    steady state bifurcations in systems of ordinary differential
  18    equations. see page 176 and page 211 in "perturbation methods,
  19    bifurcation  theory and computer algebra".*/
  20
  21
  22
  23
  24
  25 reduction2():=block(
  26            setup(),
  27            order:read("to which order do you want to calculate"),
  28            for i:2 thru order-1 do  w_poly(i,0),
  29            for i:1 thru order-2 do  w_poly(i,1),
  30            for i:1 thru order do g_poly(i,0),
  31            for i:1 thru order-1 do g_poly(i,1),
  32            g_eq())$
  33
  34
  35
  36
  37 setup():=(/*the function setup needs the dimension n of the system of
  38             o.d.e., the n-dimensional list of variables y, the system
  39             of o.d.e's called eqs, the critical eigenvector cfun and
  40             the adjoint critical eigenvector afun.*/
  41    printflag:true,
  42    ls_list:[],
  43    n:read("number of equations"),
  44    y:makelist(read("enter variable number",i),i,1,n),
  45    alpha:read("enter the bifurcation parameter"),
  46    cal:read("enter the critical bifurcation value",alpha),
  47    print("we define lam =",alpha - cal),
  48    cfun:read("enter the critical eigenvector as a list"),
  49    afun:read("enter the adjoint critical eigenvector"),
  50    kill(w,zwlist),
  51    w:makelist(concat(w,i),i,1,n),
  52    zwlist:makelist(concat(zw,i),i,1,n),
  53    depends(append(zwlist,w),[amp,lam]),
  54    print("enter the differential equation"),
  55    eqs:makelist(read("diff(",y[i],",t)="),i,1,n),
  56    eqlam:ev(eqs,ev(alpha) = lam + cal),
  57    print(eqlam),
  58    wnull:vfun(w,0),
  59    zwnull:vfun(zwlist,0),
  60    sub:maplist("=",y,amp*cfun+w),
  61    database:map(lambda([u],diff(u,amp)=0),w),
  62    database:append(database,map(lambda([u],diff(u,lam)=0),w)),
  63    norm:afun.cfun,
  64    temp1:ev(eqlam,sub,diff),
  65    print("do you know apriori that some taylor coefficients"),
  66           nullans:read(" are zero, y/n")
  67    )$
  68
  69
  70 g_poly(i,j):=block(/*provided all  necessary knowledge about the taylor
  71                      coefficients of w(amp,lam) is stored in the list
  72                      database, g_poly calculates any specific taylor
  73                      coefficient of the bifurcation equation
  74                      g(amp,lam)= 0 via the projection onto cfun.*/
  75     ls_list:cons([k=i,l=j],ls_list),
  76     if nullans = y then (
  77              zeroans:read("is ls(",i,j,") identically zero, y/n"),
  78              if zeroans = y then  return(bifeq[i,j]:0)),
  79    wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = 0),w),
  80    temp2:diff(temp1,amp,i,lam,j),
  81    temp2:subst(wmax_diff,temp2),
  82    temp3:subst(database,temp2),
  83    temp4:expand(ev(temp3,amp=0,lam=0,wnull)),
  84    bifeq[i,j]:ratsimp(temp4.afun))$
  85
  86
  87
  88 w_poly(i,j):=block(/*the function w_poly allows to calculate iteratively,
  89                      i.e. starting with the lowest order, all taylor
  90                      coefficients of w(amp,lam) by projecting onto the
  91                      complement of cfun*/
  92    if nullans = y then (
  93                   zeroans:read("is diff(w,amp,",i,"lam,",j,"
  94                                 identically zero, y/n"),
  95                   if zeroans = y then
  96                         (addbase:['diff(w,amp,i,lam,j)=0],
  97                          database:append(database,addbase),
  98                          return(addbase))),
  99
 100    wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = concat(z,u)),w),
 101    temp2:diff(temp1,amp,i,lam,j),
 102    temp2:subst(wmax_diff,temp2),
 103    temp3:subst(database,temp2),
 104    temp4:ev(temp3,amp=0,lam=0,wnull),
 105    temp5:ev(temp4,zwnull).afun,
 106    temp6:temp4-temp5/norm*cfun,
 107    temp7:solve(temp6,zwlist),
 108    temp8:ev(zwlist,temp7),
 109    temp9:ratsimp(temp8.afun),
 110    if temp9 = 0 then
 111             addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n),temp8)
 112                        else
 113                             (temp10:ratsimp(temp8 - temp9/norm*cfun),
 114                              addbase:map("=",makelist('diff(w[u],amp,i,lam,j),
 115                                                 u,1,n),temp10)),
 116    database:append(database,addbase),
 117    print(addbase))$
 118
 119
 120
 121 vfun(list,value):=map(lambda([u],u=value),list)$
 122
 123 g_eq():=  sum(ev(amp^k*lam^l/k!*bifeq[k,l],ls_list[n]),n,1,length(ls_list))$