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
39reduction1():=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
50setup():=(
51/*the function setup asks for the variables of the d.e., the bifurcation
52point, the critical eigenfunction, the x-interval, the boundary
53conditions and the differential equation. */
54assume_pos:true,
55ls_list:[],
56y:read("enter dependent variable"),
57print("use x as the independent variable and alpha as a parameter to vary"),
58cal:read("enter the critical bifurcation value alpha"),
59print("we define lam = alpha - ",cal),
60cfun: read("enter the critical eigenfunction"),
61depends([zw,y,w],[amp,lam,x]),
62len:read("what is the length of the x-interval"),
63norm:integrate(cfun^2,x,0,len),
64print("specify the boundary conditions"),
65print("your choice for the b.c. on y at x=0 and x=",len),
66print("enter 1 for y=0, 2 for y'=0"),
67bcl:read("b.c. at 0?"),
68bcr:read("b.c. at",len,"?"),
69eq:diff(y,x,2) +
70	 read("the d.e. is of the form y'' + f(y,y',alpha) = 0,enter f"),
71eqlam:ev(eq,alpha=lam+cal),
72print(eqlam),
73database:[diff(w,amp)=0,diff(w,lam)=0],
74sub:y=amp*cfun+w,
75temp1:ev(eqlam,sub,diff),
76nullans:read("do you know apriori that some taylor coefficents are zero, y/n")
77)$
78
79
80g_poly(i,j):=block(
81/* this is a function to determine a particular taylor coefficient of the
82bifurcation equation g(amp,lam) =0. it requires knowledge about the taylor
83coefficients 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
87zero, 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
92critical 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
108w_poly(i,j):=block(
109/* this function allows to determine iteratively any particular taylor
110coefficient 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
121an 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
131the 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
145project(exp):=(
146	  pro1:ev(exp,zw=0),
147	  pro2:integrate(pro1*cfun,x,0,len)/norm,
148	  expand(exp-pro2*cfun))$
149neumannbc(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))$
159g_eq():=  sum(ev(amp^k*lam^l/k!*bifeq[k,l],ls_list[n]),n,1,length(ls_list))$
160