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
37setup():=(
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
70g_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
85w_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
103feedw():=(
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
119int(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
123vfun(list,value):=map(lambda([u],u=value),list)$
124