/*PROGRAM NUMBER 13: CONTAINS THE BUILDING BLOCKS TO PERFORM STEADY
  STATE BIFURCATIONS FOR MORE PARTIAL DIFFERENTIAL EQUATIONS. SEE PAGES
  198-202 IN "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER
  ALGEBRA" */




/*THE FOLLOWING FUNCTIONS CAN BE USED TO PERFORM A LIAPUNOV-SCHMIDT REDUCTION
  FOR THE 2D BENARD PROBLEM */

/*SETUP() ALLOWS THE PROBLEM TO BE ENTERED,
  G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION
  EQUATION G(AMP,LAM),
  W_POLY(I,J) DETERMINES A DIFFERENTIAL EQUATION, WHOSE SOLUTION IS THE
  COEFFICIENT OF AMP^I LAM^J IN W(AMP,LAM),
  FEEDW() ALLOWS THIS SOLUTION TO BE ENTERED INTO THE LIST DATABASE,
  INT(EXP) FINDS MULTIPLE DEFINITE INTEGRALS EFFECTIVELY,
  VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST 
		[LIST[1] = VALUE, LIST[2] = VALUE ...]
*/

SETUP():=( 
       N:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"),
       Y:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"),
       SPACE:READ("ENTER NUMBER OF SPATIAL COORDINATES"),
       XVAR:READ("ENTER THE SPATIAL COORDINATES AS A LIST"),
       ALPHA:READ("ENTER THE BIFURCATION PARAMETER"),
       CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"),
       PRINT("WE DEFINE LAM = ",ALPHA - CAL),
       CFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"),
       AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS ALIST"),
       KILL(W),
       W:MAKELIST(CONCAT(W,I),I,1,N),
       ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N),
       DEPENDS(APPEND(ZWLIST,W,Y),APPEND([AMP,LAM],XVAR)),
       EQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N),
       EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL),
       PRINT("ENTER THE LOWER LEFT CORNER OF THE",N,"DIMENSIONAL SPACE 
						INTERVAL"),
       LBOUND:READ("[",X[1],"=...,]"),
       PRINT("ENTER THE UPPER RIGHT CORNER OF THE",N,"DIMENSIONAL SPACE 
						INTERVAL"),
       UBOUND:READ("[",X[1],"=...,]"),
       WNULL:VFUN(W,0),
       SUB:MAPLIST("=",Y,AMP*CFUN+W),
       DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W),
       ZWNULL:VFUN(ZWLIST,0),
       NORM:INT(AFUN.CFUN),
       TEMP1:EV(EQLAM,SUB,DIFF),
       EQLAM
       )$



G_POLY(I,J):=BLOCK(
       WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W),
       TEMP2:DIFF(TEMP1,AMP,I,LAM,J),
       DERIVSUBST:TRUE,
       TEMP3:SUBST(WMAX_DIFF,TEMP2),
       M:LENGTH(DATABASE),
       FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF),
       DERIVSUBST:FALSE,
       TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)),
       TEMP5:RATSIMP(TEMP4.AFUN),
       BIFEQ:INT(TEMP5))$




W_POLY(I,J):=BLOCK(
       WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W),
       TEMP2:DIFF(TEMP1,AMP,I,LAM,J),
       DERIVSUBST:TRUE,
       TEMP3:SUBST(WMAX_DIFF,TEMP2),
       M:LENGTH(DATABASE),
       FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF),
       DERIVSUBST:FALSE,
       TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)),
       TEMP5:INT(EV(TEMP4,ZWNULL).AFUN),
       TEMP6:TEMP4-TEMP5/NORM*CFUN,  
       FOR I THRU SPACE DO TEMP7:TRIGREDUCE(TEMP6,XVAR[I]),
       W_DE:RATSIMP(TEMP7 ),
       PRINT("NOW SOLVE THE EQUATIONS",W_DE,"=0! THEY ARE GIVEN IN W_DE"),
       PRINT("CALL FEEDW() TO PROCEED"))$ 



FEEDW():=(
       ADDBASE:[],
       RESULT:MAKELIST((PRINT("ENTER",ZWLIST[K]),READ()),K,1,N),
       F2:INT(RESULT.AFUN),
       IF RATSIMP(F2)=0 
	        THEN 
                     ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),
							RESULT)
	        ELSE (ORTHO_RESULT:RATSIMP(RESULT- F2/NORM*CFUN),
/*PROJECTION Q*/
		     ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),
						ORTHO_RESULT)),
       DATABASE:APPEND(DATABASE,ADDBASE),
       ADDBASE)$

 
INT(EXP):=(%INT:EXP,FOR I THRU SPACE DO %INT:INTEGRATE(TRIGREDUCE(%INT,
					XVAR[I]),XVAR[I]),
			RATSIMP(EV(%INT,UBOUND)-EV(%INT,LBOUND)))$

VFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$