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