share/contrib/rand/benard.mac

   1 /* Filename benard.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 */ /*program number 13: contains the building blocks to perform steady
  16   state bifurcations for more partial differential equations. see pages
  17   198-202 in "perturbation methods, bifurcation theory and computer
  18   algebra" */
  19
  20
  21
  22
  23 /*the following functions can be used to perform a liapunov-schmidt reduction
  24   for the 2d benard problem */
  25
  26 /*setup() allows the problem to be entered,
  27   g_poly(i,j) calculates the coefficient of amp^i lam^j in the bifurcation
  28   equation g(amp,lam),
  29   w_poly(i,j) determines a differential equation, whose solution is the
  30   coefficient of amp^i lam^j in w(amp,lam),
  31   feedw() allows this solution to be entered into the list database,
  32   int(exp) finds multiple definite integrals effectively,
  33   vfun(list,value) creates the substitution list
  34                 [list[1] = value, list[2] = value ...]
  35 */
  36
  37 setup():=(
  38        n:read("enter the number of differential equations"),
  39        y:read("enter the dependent variables as a list"),
  40        space:read("enter number of spatial coordinates"),
  41        xvar:read("enter the spatial coordinates as a list"),
  42        alpha:read("enter the bifurcation parameter"),
  43        cal:read("enter the critical bifurcation value"),
  44        print("we define lam = ",alpha - cal),
  45        cfun:read("enter the critical eigenfunction as a list"),
  46        afun:read("enter the adjoint critical eigenfunction as alist"),
  47        kill(w),
  48        w:makelist(concat(w,i),i,1,n),
  49        zwlist:makelist(concat(zw,i),i,1,n),
  50        depends(append(zwlist,w,y),append([amp,lam],xvar)),
  51        eqs:makelist(read("enter the differential equation number",i),i,1,n),
  52        eqlam:ev(eqs,ev(alpha) = lam + cal),
  53        print("enter the lower left corner of the",n,"dimensional space
  54                                                 interval"),
  55        lbound:read("[",x[1],"=...,]"),
  56        print("enter the upper right corner of the",n,"dimensional space
  57                                                 interval"),
  58        ubound:read("[",x[1],"=...,]"),
  59        wnull:vfun(w,0),
  60        sub:maplist("=",y,amp*cfun+w),
  61        database:map(lambda([u],diff(u,amp)=0),w),
  62        zwnull:vfun(zwlist,0),
  63        norm:int(afun.cfun),
  64        temp1:ev(eqlam,sub,diff),
  65        eqlam
  66        )$
  67
  68
  69
  70 g_poly(i,j):=block(
  71        wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = 0),w),
  72        temp2:diff(temp1,amp,i,lam,j),
  73        derivsubst:true,
  74        temp3:subst(wmax_diff,temp2),
  75        m:length(database),
  76        for i thru m do temp3:ev(subst(database[m+1-i],temp3),diff),
  77        derivsubst:false,
  78        temp4:expand(ev(temp3,amp=0,lam=0,wnull)),
  79        temp5:ratsimp(temp4.afun),
  80        bifeq:int(temp5))$
  81
  82
  83
  84
  85 w_poly(i,j):=block(
  86        wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = concat(z,u)),w),
  87        temp2:diff(temp1,amp,i,lam,j),
  88        derivsubst:true,
  89        temp3:subst(wmax_diff,temp2),
  90        m:length(database),
  91        for i thru m do temp3:ev(subst(database[m+1-i],temp3),diff),
  92        derivsubst:false,
  93        temp4:expand(ev(temp3,amp=0,lam=0,wnull)),
  94        temp5:int(ev(temp4,zwnull).afun),
  95        temp6:temp4-temp5/norm*cfun,
  96        for i thru space do temp7:trigreduce(temp6,xvar[i]),
  97        w_de:ratsimp(temp7 ),
  98        print("now solve the equations",w_de,"=0! they are given in w_de"),
  99        print("call feedw() to proceed"))$
 100
 101
 102
 103 feedw():=(
 104        addbase:[],
 105        result:makelist((print("enter",zwlist[k]),read()),k,1,n),
 106        f2:int(result.afun),
 107        if ratsimp(f2)=0
 108                 then
 109                      addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n),
 110                                                         result)
 111                 else (ortho_result:ratsimp(result- f2/norm*cfun),
 112 /*projection q*/
 113                      addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n),
 114                                                 ortho_result)),
 115        database:append(database,addbase),
 116        addbase)$
 117
 118
 119 int(exp):=(%int:exp,for i thru space do %int:integrate(trigreduce(%int,
 120                                         xvar[i]),xvar[i]),
 121                         ratsimp(ev(%int,ubound)-ev(%int,lbound)))$
 122
 123 vfun(list,value):=map(lambda([u],u=value),list)$