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