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