1/* Filename reduct3.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 12: a liapunov-schmidt reduction for steady state 17 bifurcations in systems of partial differential equations 18 depending on one independent space variable. see page 189 in 19 "perturbation methods, bifurcation theory and computer algebra". 20*/ 21 22 23 24 25/*this file contains reduction3(), a function to perform a liapunov- 26schmidt reduction for steady state bifurcations from n coupled 27partial differential equations on one spatial dimension. 28the following additional functions are supplied: 29 setup() allows the problem to be entered. 30 g_poly(i,j) calculates the coefficient of amp^i lam^j in the bifurcation 31 equation g(amp,lam). 32 w_poly(i,j) calculates the coefficient of amp^i lam^j in w(x;amp,lam). 33 solve_ode(exp) solves certain ordinary differential equations via a fourier 34 mode ansatz. 35 feedw(exp) ensures that <afun,w> = 0 . 36 find_trig(exp) identifies fourier modes. 37 setify(list) transforms a list into a set. 38 g_eq() assembles the bifurcation equation. 39 vfun(list,value) creates the substitution list 40 [list[1] = value, list[2] = value, ...] 41 diffnull(i,j) sets the derivative diff(w,amp,i,lam,j) to zero.*/ 42 43reduction3():=block( 44 setup(), 45 order:read("to which order do you want to calculate"), 46 for i:2 thru order-1 do w_poly(i,0), 47 for i:1 thru order-2 do w_poly(i,1), 48 for i:1 thru order do g_poly(i,0), 49 for i:1 thru order-1 do g_poly(i,1), 50 g_eq())$ 51 52 53setup():=( /* this function performs the input for the liapunov-schmidt 54 reduction*/ 55assume_pos:true, 56ls_list:[], 57n:read("enter the number of differential equations"), 58y:read("enter the dependent variables as a list"), 59xvar:read("enter the spatial coordinate"), 60alpha:read("enter the bifurcation parameter"), 61cal:read("enter the critical bifurcation value"), 62print("we define lam = ",alpha - cal), 63cfun:read("enter the critical eigenfunction as a list"), 64afun:read("enter the adjoint critical eigenfunction as a list"), 65kill(w), 66w:makelist(concat(w,i),i,1,n), 67zwlist:makelist(concat(zw,i),i,1,n), 68nullist:makelist(0,i,1,n), 69depends(append(zwlist,w,y),cons(xvar,[amp,lam])), 70eqs:makelist(read("enter the differential equation number",i),i,1,n), 71eqlam:ev(eqs,ev(alpha) = lam + cal), 72print(eqlam), 73len:read("what is the length of the space interval"), 74wnull:vfun(w,0), 75sub:maplist("=",y,amp*cfun+w), 76database:append(difnull(1,0),difnull(0,1)), 77zwnull:vfun(zwlist,0), 78norm:integrate(afun.cfun,xvar,0,len), 79temp1:ev(eqlam,sub,diff), 80print("do you know apriori that some taylor coefficients are 0"), 81nullans:read("y,n") 82)$ 83 84 85g_poly(i,j):=block( 86/*this is a function to determine a particular taylor coefficient of the 87bifurcation equation g(amp,lam) = 0. it requires kowledge about the taylor 88coefficients of w(amp,lam). this knowledge is stored in the list database*/ 89 ls_list:cons([k=i,l=j],ls_list), 90 if nullans = y then ( 91 zeroans:read("is g_poly(",i,",",j,")identically 92 zero, y/n"), 93 if zeroans = y then return(bifeq[i,j]:0)), 94 temp2:diff(temp1,amp,i,lam,j), 95 derivsubst:true, 96 /*set the derivatives diff(w,amp,i,lam,j) to zero*/ 97 temp3:subst(difnull(i,j),temp2), 98 /*enter all information in database*/ 99 d_base_length:length(database), 100 for ii thru d_base_length do 101 temp3:ev(subst(database[d_base_length+1-ii],temp3),diff), 102 derivsubst:false, 103 temp4:expand(ev(temp3,amp=0,lam=0,wnull,integrate)), 104 /*project onto afun*/ 105 temp5:ratsimp(temp4.afun), 106 bifeq[i,j]:integrate(trigreduce(temp5),xvar,0,len))$ 107 108 109 110 111w_poly(i,j):=block( 112/* this function allows to determine iteratively any particular taylor 113coefficient of the function w(amp,lam). it returns a differential equation 114for the particular coefficient of interest (called zw1,zw2...). this d.e. 115is solved via solve_ode and the result is fed into database from feedw*/ 116 if nullans = y then ( 117 zeroans:read("is diff(w(amp,",i,",lam,",j,") identically 118 zero, y/n"), 119 if zeroans = y then ( 120 addbase:difnull(i,j), 121 database:append(database,addbase), 122 return(addbase))), 123 wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = concat(z,u)),w), 124 temp2:diff(temp1,amp,i,lam,j), 125 derivsubst:true, 126 /*rename the derivatives diff(w,amp,i,lam,j) to zw */ 127 temp3:subst(wmax_diff,temp2), 128 /*enter all information in stored in database*/ 129 d_base_length:length(database), 130 for ii thru d_base_length do 131 temp3:ev(subst(database[d_base_length+1-ii],temp3),diff), 132 derivsubst:false, 133 temp4:ev(temp3,amp=0,lam=0,wnull,integrate), 134 /*this is the projection q onto the range of the 135 linear differential operator in the problem*/ 136 temp5:integrate(ev(temp4,zwnull).afun,xvar,0,len), 137 temp6:temp4-temp5/norm*cfun, 138 temp7:trigreduce(temp6), 139 /*the set of o.d.e.'s to solve */ 140 w_de:expand(temp7), 141 temp8:ev(w_de,vfun(zwlist,0)), 142 /*if the particular solution of w_de is zero then w=0 */ 143 if nullist=temp8 then ( addbase:difnull(i,j), 144 database:append(database,addbase), 145 return(addbase)), 146 temp9:solve_ode(temp8), 147 feedw(temp9) 148 )$ 149 150 151solve_ode(exp):=( 152/*this function solve the d.e. w_de by a fourier mode ansatz*/ 153 trigfun:[], 154 const:false, 155 for i thru n do ( 156 /*determine the fourier modes*/ 157 trig1:exp[i], 158 if trig1 # 0 then ( 159 trig2:apply1(trig1,sinnull,cosnull), 160 if trig2 # 0 then ( 161 const:true, 162 trig1:trig1-trig2), 163 trigfun:append(find_trig(trig1),trigfun))), 164 sol1:delete(dummy,setify(trigfun)), 165 /*make an ansatz*/ 166 ansatz:makelist(sum(am[i,j]*sol1[i],i,1,length(sol1)),j,1,n), 167 sol2:ev(w_de,map("=",zwlist,ansatz),diff), 168 sol3:makelist(ratcoef(sol2,i),i,sol1), 169 eqlist:[], 170 for i thru length(sol3) do eqlist:append(eqlist,sol3[i]), 171 varlist:[], 172 for i thru n do 173 for j thru length(sol1) do varlist:cons(am[j,i],varlist), 174 /*find the amplitudes of the fourier modes*/ 175 sol4:solve(eqlist,varlist), 176 /*solve for the constant fourier mode if necessary*/ 177 cansatz:0, 178 if const = true then (cansatz:makelist(concat(con,i),i,1,n), 179 csol1:ev(w_de,map("=",zwlist,cansatz),diff), 180 csol2:apply1(csol1,sinnull,cosnull), 181 csol3:solve(csol2,cansatz)), 182 ev(ansatz+cansatz,sol4,csol3))$ 183 184 185feedw(exp):=( 186/* this function allows the result of the ode-solver to be entered into 187 the list database. it checks for orthogonality to the critical adjoint 188 eigenfunction and removes solutions of the homogeneous equation (i.e. 189 nonorthogonal terms)*/ 190 f2:integrate(exp.afun,xvar,0,len), 191 if ratsimp(f2)=0 192 then 193 addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n), 194 exp) 195 else (ortho_result:ratsimp(exp- f2/norm*cfun), 196 addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n), 197 ortho_result)), 198 database:append(database,addbase), 199 print(addbase))$ 200/*collect all information stored in bifeq and assemble the bifurcation 201equation*/ 202g_eq():= sum(ev(amp^k*lam^l/k!*bifeq[k,l],ls_list[n]),n,1,length(ls_list))$ 203 204 205setify(list):=(/*transforms a list into a set, i.e. removes duplicates*/ 206 set:[list[1]], 207 for i :2 thru length(list) do 208 if not member(list[i],set) then 209 set:cons(list[i],set), 210 set)$ 211 212find_trig(exp):=(/*finds the fourier modes*/ 213 f_a1:args(exp+dummy), 214 f_a2:apply1(f_a1,sinfind,cosfind) 215 )$ 216 217/*auxiliary functions */ 218matchdeclare([dummy1,dummy2],true)$ 219defrule(cosfind,dummy1*cos(dummy2),cos(dummy2))$ 220defrule(sinfind,dummy1*sin(dummy2),sin(dummy2))$ 221defrule(cosnull,cos(dummy1),0)$ 222defrule(sinnull,sin(dummy1),0)$ 223 224 225vfun(list,value):=map(lambda([u],u=value),list)$ 226difnull(i,j):=map(lambda([u],'diff(u,amp,i,lam,j)=0),w) $ 227