share/contrib/rand/reduct1.mac

   1 /* Filename reduct1.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 10: reduction1(), a liapunov-schmidt reduction for
  17   steady state bifurcations in one differential equation depending
  18   on one independent variable. see page 170 in "perturbation methods,
  19   bifurcation theory and computer algebra". */
  20
  21
  22
  23
  24 /*this file contains reduction1(), a function to perform a liapunov-
  25   schmidt reduction for steady state bifurcations of nonlinear d.e.'s
  26   of the form y'' + f(y,y',alpha) = 0. y = y(x) is defined on a real
  27   interval with dirichlet or neumann boundary conditions and f depends
  28   only linearly on alpha.
  29   it also contains these additional functions:
  30   setup() allows to enter the problem.
  31   g_poly(i,j) calculates the coefficient of amp^i lam^j in the bifurcation
  32   equation.
  33   w_poly(i,j) calculates the coefficient of amp^i lam^j in w(x;amp,lam).
  34   project(exp) ensures that <cfun,exp>=0.
  35   neumannbc(exp) accounts for neumann boundary conditions.
  36   g_eq() assembles the bifurcation equation. */
  37
  38
  39 reduction1():=block(
  40            setup(),
  41            order:read("to which order do you want to calculate"),
  42            for i:2 thru order-1 do  w_poly(i,0),
  43            for i:1 thru order-2 do  w_poly(i,1),
  44            for i:1 thru order do g_poly(i,0),
  45            for i:1 thru order-1 do g_poly(i,1),
  46            g_eq())$
  47
  48
  49
  50 setup():=(
  51 /*the function setup asks for the variables of the d.e., the bifurcation
  52 point, the critical eigenfunction, the x-interval, the boundary
  53 conditions and the differential equation. */
  54 assume_pos:true,
  55 ls_list:[],
  56 y:read("enter dependent variable"),
  57 print("use x as the independent variable and alpha as a parameter to vary"),
  58 cal:read("enter the critical bifurcation value alpha"),
  59 print("we define lam = alpha - ",cal),
  60 cfun: read("enter the critical eigenfunction"),
  61 depends([zw,y,w],[amp,lam,x]),
  62 len:read("what is the length of the x-interval"),
  63 norm:integrate(cfun^2,x,0,len),
  64 print("specify the boundary conditions"),
  65 print("your choice for the b.c. on y at x=0 and x=",len),
  66 print("enter 1 for y=0, 2 for y'=0"),
  67 bcl:read("b.c. at 0?"),
  68 bcr:read("b.c. at",len,"?"),
  69 eq:diff(y,x,2) +
  70          read("the d.e. is of the form y'' + f(y,y',alpha) = 0,enter f"),
  71 eqlam:ev(eq,alpha=lam+cal),
  72 print(eqlam),
  73 database:[diff(w,amp)=0,diff(w,lam)=0],
  74 sub:y=amp*cfun+w,
  75 temp1:ev(eqlam,sub,diff),
  76 nullans:read("do you know apriori that some taylor coefficents are zero, y/n")
  77 )$
  78
  79
  80 g_poly(i,j):=block(
  81 /* this is a function to determine a particular taylor coefficient of the
  82 bifurcation equation g(amp,lam) =0. it requires knowledge about the taylor
  83 coefficients of w(amp,lam). this knowledge is stored in the list database */
  84           ls_list:cons([k=i,l=j],ls_list),
  85           if nullans = y then (
  86                                zeroans:read("is g_poly(",i,",",j,") identically
  87 zero, y/n"),
  88                                if zeroans = y then  return(bifeq[i,j]:0)),
  89           temp2:diff(temp1,amp,i,lam,j),
  90           derivsubst:true,
  91 /*this differential of w will be annihilated by the projection onto the
  92 critical eigenfunction, hence we set it to zero already here */
  93           temp3:subst('diff(w,amp,i,lam,j)=0,temp2),
  94 /* here we insert all knowledge gained through w_poly */
  95           d_base_length:length(database),
  96           for ii thru d_base_length do
  97                      temp3:ev(subst(database[d_base_length+1-ii],temp3),diff),
  98           derivsubst:false,
  99           temp4:ev(temp3,amp=0,lam=0,w=0),
 100 /*projection onto cfun, the critical eigenfunction */
 101           temp5:trigreduce(cfun*temp4),
 102           bifeq[i,j]:ratsimp(integrate(temp5,x,0,len))
 103           )$
 104
 105
 106
 107
 108 w_poly(i,j):=block(
 109 /* this function allows to determine iteratively any particular taylor
 110 coefficient of the function w(amp,lam). the result is fed into database.*/
 111           if nullans = y then (
 112                    zeroans:read("is diff(w,amp,",i,",lam,",j,") identically zero
 113 , y/n"),
 114                            if zeroans = y then
 115                                     (addbase:['diff(w,amp,i,lam,j)=0],
 116                                      database:append(database,addbase),
 117                                      return(addbase))),
 118           temp2:diff(temp1,amp,i,lam,j),
 119           derivsubst:true,
 120 /*we substitute zw for the unknown taylor coefficient and solve
 121 an o.d.e. for zw in temp7 */
 122           temp3:subst(diff(w,amp,i,lam,j)=zw,temp2),
 123 /*now we insert all knowledge gained through previous iterations*/
 124           d_base_length:length(database),
 125           for ii thru d_base_length do
 126                      temp3:ev(subst(database[d_base_length+1-ii],temp3),diff),
 127           derivsubst:false,
 128           temp4:ev(temp3,amp=0,lam=0,w=0),
 129           temp5:trigreduce(temp4),
 130 /*this ensures that the r.h.s. fo the d.e. temp6 is orthogonal to
 131 the solution of the homogeneous equation */
 132           temp6:project(temp5),
 133           temp7:ode2(temp6,zw,x),
 134 /*satisfy boundary conditions*/
 135           if bcl*bcr=1 then temp8:bc2(temp7,x=0,zw=0,x=len,zw=0)
 136                         else
 137                              temp8:neumannbc(temp7),
 138 /*make sure that <cfun,w>= 0*/
 139           temp9:project(temp8),
 140           addbase:['diff(w,amp,i,lam,j)=rhs(temp9)],
 141           database:append(database,addbase),
 142           print(addbase))$
 143
 144
 145 project(exp):=(
 146           pro1:ev(exp,zw=0),
 147           pro2:integrate(pro1*cfun,x,0,len)/norm,
 148           expand(exp-pro2*cfun))$
 149 neumannbc(exp):=(
 150           rexp:rhs(exp),
 151           nbc1:diff(rexp,x),
 152           if bcl=1 and bcr=2 then nbc2:solve([ev(rexp,x=0),ev(nbc1,x=len)],
 153                 [%k1,%k2]),
 154           if bcl=2 and bcr=1 then nbc2:solve([ev(rexp,x=len),ev(nbc1,x=0)],
 155                 [%k1,%k2]),
 156           if bcl=2 and bcr=2 then nbc2:solve([ev(nbc1,x=len),ev(nbc1,x=0)],
 157                 [%k1,%k2]),
 158           ev(exp,nbc2))$
 159 g_eq():=  sum(ev(amp^k*lam^l/k!*bifeq[k,l],ls_list[n]),n,1,length(ls_list))$