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