1/* Filename reduct2.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 11: reduction2(), a liapunov-schmidt reduction for
17   steady state bifurcations in systems of ordinary differential
18   equations. see page 176 and page 211 in "perturbation methods,
19   bifurcation  theory and computer algebra".*/
20
21
22
23
24
25reduction2():=block(
26	   setup(),
27	   order:read("to which order do you want to calculate"),
28	   for i:2 thru order-1 do  w_poly(i,0),
29	   for i:1 thru order-2 do  w_poly(i,1),
30	   for i:1 thru order do g_poly(i,0),
31	   for i:1 thru order-1 do g_poly(i,1),
32	   g_eq())$
33
34
35
36
37setup():=(/*the function setup needs the dimension n of the system of
38            o.d.e., the n-dimensional list of variables y, the system
39	    of o.d.e's called eqs, the critical eigenvector cfun and
40	    the adjoint critical eigenvector afun.*/
41   printflag:true,
42   ls_list:[],
43   n:read("number of equations"),
44   y:makelist(read("enter variable number",i),i,1,n),
45   alpha:read("enter the bifurcation parameter"),
46   cal:read("enter the critical bifurcation value",alpha),
47   print("we define lam =",alpha - cal),
48   cfun:read("enter the critical eigenvector as a list"),
49   afun:read("enter the adjoint critical eigenvector"),
50   kill(w,zwlist),
51   w:makelist(concat(w,i),i,1,n),
52   zwlist:makelist(concat(zw,i),i,1,n),
53   depends(append(zwlist,w),[amp,lam]),
54   print("enter the differential equation"),
55   eqs:makelist(read("diff(",y[i],",t)="),i,1,n),
56   eqlam:ev(eqs,ev(alpha) = lam + cal),
57   print(eqlam),
58   wnull:vfun(w,0),
59   zwnull:vfun(zwlist,0),
60   sub:maplist("=",y,amp*cfun+w),
61   database:map(lambda([u],diff(u,amp)=0),w),
62   database:append(database,map(lambda([u],diff(u,lam)=0),w)),
63   norm:afun.cfun,
64   temp1:ev(eqlam,sub,diff),
65   print("do you know apriori that some taylor coefficents"),
66	  nullans:read(" are zero, y/n")
67   )$
68
69
70g_poly(i,j):=block(/*provided all  necessary knowledge about the taylor
71		     coefficients of w(amp,lam) is stored in the list
72		     database, g_poly calculates any specific taylor
73		     coefficient of the bifurcation equation
74		     g(amp,lam)= 0 via the projection onto cfun.*/
75    ls_list:cons([k=i,l=j],ls_list),
76    if nullans = y then (
77	     zeroans:read("is ls(",i,j,") identically zero, y/n"),
78	     if zeroans = y then  return(bifeq[i,j]:0)),
79   wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = 0),w),
80   temp2:diff(temp1,amp,i,lam,j),
81   temp2:subst(wmax_diff,temp2),
82   temp3:subst(database,temp2),
83   temp4:expand(ev(temp3,amp=0,lam=0,wnull)),
84   bifeq[i,j]:ratsimp(temp4.afun))$
85
86
87
88w_poly(i,j):=block(/*the function w_poly allows to calculate iteratively,
89		     i.e. starting with the lowest order, all taylor
90		     coefficients of w(amp,lam) by projecting onto the
91		     complement of cfun*/
92   if nullans = y then (
93		  zeroans:read("is diff(w,amp,",i,"lam,",j,"
94                                identically zero, y/n"),
95		  if zeroans = y then
96			(addbase:['diff(w,amp,i,lam,j)=0],
97			 database:append(database,addbase),
98			 return(addbase))),
99
100   wmax_diff:map(lambda([u],'diff(u,amp,i,lam,j) = concat(z,u)),w),
101   temp2:diff(temp1,amp,i,lam,j),
102   temp2:subst(wmax_diff,temp2),
103   temp3:subst(database,temp2),
104   temp4:ev(temp3,amp=0,lam=0,wnull),
105   temp5:ev(temp4,zwnull).afun,
106   temp6:temp4-temp5/norm*cfun,
107   temp7:solve(temp6,zwlist),
108   temp8:ev(zwlist,temp7),
109   temp9:ratsimp(temp8.afun),
110   if temp9 = 0 then
111            addbase:map("=",makelist('diff(w[u],amp,i,lam,j),u,1,n),temp8)
112		       else
113			    (temp10:ratsimp(temp8 - temp9/norm*cfun),
114			     addbase:map("=",makelist('diff(w[u],amp,i,lam,j),
115						u,1,n),temp10)),
116   database:append(database,addbase),
117   print(addbase))$
118
119
120
121vfun(list,value):=map(lambda([u],u=value),list)$
122
123g_eq():=  sum(ev(amp^k*lam^l/k!*bifeq[k,l],ls_list[n]),n,1,length(ls_list))$
124