/* PROGRAM NUMBER 11: REDUCTION2(), A LIAPUNOV-SCHMIDT REDUCTION FOR
   STEADY STATE BIFURCATIONS IN SYSTEMS OF ORDINARY DIFFERENTIAL
   EQUATIONS. SEE PAGE 176 AND PAGE 211 IN "PERTURBATION METHODS, 
   BIFURCATION  THEORY AND COMPUTER ALGEBRA".*/





REDUCTION2():=BLOCK(
	   SETUP(),
	   ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"),
	   FOR I:2 THRU ORDER-1 DO  W_POLY(I,0),
	   FOR I:1 THRU ORDER-2 DO  W_POLY(I,1),
	   FOR I:1 THRU ORDER DO G_POLY(I,0),
	   FOR I:1 THRU ORDER-1 DO G_POLY(I,1),
	   G_EQ())$




SETUP():=(/*THE FUNCTION SETUP NEEDS THE DIMENSION N OF THE SYSTEM OF
            O.D.E., THE N-DIMENSIONAL LIST OF VARIABLES Y, THE SYSTEM
	    OF O.D.E'S CALLED EQS, THE CRITICAL EIGENVECTOR CFUN AND
	    THE ADJOINT CRITICAL EIGENVECTOR AFUN.*/
   PRINTFLAG:TRUE,
   LS_LIST:[],
   N:READ("NUMBER OF EQUATIONS"),
   Y:MAKELIST(READ("ENTER VARIABLE NUMBER",I),I,1,N),
   ALPHA:READ("ENTER THE BIFURCATION PARAMETER"),
   CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE",ALPHA),
   PRINT("WE DEFINE LAM =",ALPHA - CAL),
   CFUN:READ("ENTER THE CRITICAL EIGENVECTOR AS A LIST"),
   AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENVECTOR"),
   KILL(W,ZWLIST),
   W:MAKELIST(CONCAT(W,I),I,1,N),
   ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N),
   DEPENDS(APPEND(ZWLIST,W),[AMP,LAM]),
   PRINT("ENTER THE DIFFERENTIAL EQUATION"),
   EQS:MAKELIST(READ("DIFF(",Y[I],",T)="),I,1,N),
   EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL),
   PRINT(EQLAM),
   WNULL:VFUN(W,0),
   ZWNULL:VFUN(ZWLIST,0),
   SUB:MAPLIST("=",Y,AMP*CFUN+W),
   DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W),
   DATABASE:APPEND(DATABASE,MAP(LAMBDA([U],DIFF(U,LAM)=0),W)),
   NORM:AFUN.CFUN,
   TEMP1:EV(EQLAM,SUB,DIFF),
   PRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"),
	  NULLANS:READ(" ARE ZERO, Y/N")
   )$


G_POLY(I,J):=BLOCK(/*PROVIDED ALL  NECESSARY KNOWLEDGE ABOUT THE TAYLOR
		     COEFFICIENTS OF W(AMP,LAM) IS STORED IN THE LIST
		     DATABASE, G_POLY CALCULATES ANY SPECIFIC TAYLOR
		     COEFFICIENT OF THE BIFURCATION EQUATION 
		     G(AMP,LAM)= 0 VIA THE PROJECTION ONTO CFUN.*/
    LS_LIST:CONS([K=I,L=J],LS_LIST),
    IF NULLANS = Y THEN (
	     ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, Y/N"),
	     IF ZEROANS = Y THEN  RETURN(BIFEQ[I,J]:0)),
   WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W),
   TEMP2:DIFF(TEMP1,AMP,I,LAM,J),
   TEMP2:SUBST(WMAX_DIFF,TEMP2),
   TEMP3:SUBST(DATABASE,TEMP2),
   TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)),
   BIFEQ[I,J]:RATSIMP(TEMP4.AFUN))$



W_POLY(I,J):=BLOCK(/*THE FUNCTION W_POLY ALLOWS TO CALCULATE ITERATIVELY,
		     I.E. STARTING WITH THE LOWEST ORDER, ALL TAYLOR
		     COEFFICIENTS OF W(AMP,LAM) BY PROJECTING ONTO THE
		     COMPLEMENT OF CFUN*/
   IF NULLANS = Y THEN (
		  ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," 
                                IDENTICALLY ZERO, Y/N"),
		  IF ZEROANS = Y THEN 
			(ADDBASE:['DIFF(W,AMP,I,LAM,J)=0],
			 DATABASE:APPEND(DATABASE,ADDBASE),
			 RETURN(ADDBASE))),

   WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W),
   TEMP2:DIFF(TEMP1,AMP,I,LAM,J),
   TEMP2:SUBST(WMAX_DIFF,TEMP2),
   TEMP3:SUBST(DATABASE,TEMP2),		
   TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL),
   TEMP5:EV(TEMP4,ZWNULL).AFUN,
   TEMP6:TEMP4-TEMP5/NORM*CFUN,
   TEMP7:SOLVE(TEMP6,ZWLIST),
   TEMP8:EV(ZWLIST,TEMP7),
   TEMP9:RATSIMP(TEMP8.AFUN),
   IF TEMP9 = 0 THEN 
            ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),TEMP8)
		       ELSE 
			    (TEMP10:RATSIMP(TEMP8 - TEMP9/NORM*CFUN),
			     ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),
						U,1,N),TEMP10)),
   DATABASE:APPEND(DATABASE,ADDBASE),
   PRINT(ADDBASE))$



VFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$
 
G_EQ():=  SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$