/* PROGRAM NUMBER 12: A LIAPUNOV-SCHMIDT REDUCTION FOR STEADY STATE
   BIFURCATIONS IN SYSTEMS OF PARTIAL DIFFERENTIAL EQUATIONS 
   DEPENDING ON ONE INDEPENDENT SPACE VARIABLE. SEE PAGE 189 IN
   "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA".
*/




/*THIS FILE CONTAINS REDUCTION3(), A FUNCTION TO PERFORM A LIAPUNOV-
SCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS FROM N COUPLED
PARTIAL DIFFERENTIAL EQUATIONS ON ONE SPATIAL DIMENSION. 
THE FOLLOWING ADDITIONAL FUNCTIONS ARE SUPPLIED:
  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) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM).
  SOLVE_ODE(EXP) SOLVES CERTAIN ORDINARY DIFFERENTIAL EQUATIONS VIA A FOURIER
  MODE ANSATZ.
  FEEDW(EXP) ENSURES THAT <AFUN,W> = 0 .
  FIND_TRIG(EXP) IDENTIFIES FOURIER MODES.
  SETIFY(LIST) TRANSFORMS A LIST INTO A SET.
  G_EQ() ASSEMBLES THE BIFURCATION EQUATION.
  VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST 
          [LIST[1] = VALUE, LIST[2] = VALUE, ...]
  DIFFNULL(I,J) SETS THE DERIVATIVE DIFF(W,AMP,I,LAM,J) TO ZERO.*/

REDUCTION3():=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():=( /* THIS FUNCTION PERFORMS THE INPUT FOR THE LIAPUNOV-SCHMIDT 
	      REDUCTION*/
ASSUME_POS:TRUE,
LS_LIST:[],
N:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"),
Y:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"),
XVAR:READ("ENTER THE SPATIAL COORDINATE"),
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 A LIST"),
KILL(W),
W:MAKELIST(CONCAT(W,I),I,1,N),
ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N),
NULLIST:MAKELIST(0,I,1,N),
DEPENDS(APPEND(ZWLIST,W,Y),CONS(XVAR,[AMP,LAM])),
EQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N),
EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL),
PRINT(EQLAM),
LEN:READ("WHAT IS THE LENGTH OF THE SPACE INTERVAL"),
WNULL:VFUN(W,0),
SUB:MAPLIST("=",Y,AMP*CFUN+W),
DATABASE:APPEND(DIFNULL(1,0),DIFNULL(0,1)),
ZWNULL:VFUN(ZWLIST,0),
NORM:INTEGRATE(AFUN.CFUN,XVAR,0,LEN),
TEMP1:EV(EQLAM,SUB,DIFF),
PRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICIENTS ARE 0"),
NULLANS:READ("Y,N")
)$


G_POLY(I,J):=BLOCK(
/*THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE
BIFURCATION EQUATION G(AMP,LAM) = 0. IT REQUIRES KOWLEDGE ABOUT THE TAYLOR 
COEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE*/
       LS_LIST:CONS([K=I,L=J],LS_LIST),
       IF NULLANS = Y THEN (
			    ZEROANS:READ("IS G_POLY(",I,",",J,")IDENTICALLY 
						ZERO, Y/N"),
			    IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)),
       TEMP2:DIFF(TEMP1,AMP,I,LAM,J),
       DERIVSUBST:TRUE,
       /*SET THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZERO*/
       TEMP3:SUBST(DIFNULL(I,J),TEMP2),
       /*ENTER ALL INFORMATION IN DATABASE*/
       D_BASE_LENGTH:LENGTH(DATABASE),
       FOR II THRU D_BASE_LENGTH DO
                TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF),
       DERIVSUBST:FALSE,
       TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE)),
       /*PROJECT ONTO AFUN*/
       TEMP5:RATSIMP(TEMP4.AFUN),
       BIFEQ[I,J]:INTEGRATE(TRIGREDUCE(TEMP5),XVAR,0,LEN))$




W_POLY(I,J):=BLOCK(
/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR 
COEFFICIENT OF THE FUNCTION W(AMP,LAM). IT RETURNS A DIFFERENTIAL EQUATION 
FOR THE PARTICULAR COEFFICIENT OF INTEREST (CALLED ZW1,ZW2...). THIS D.E. 
IS SOLVED VIA SOLVE_ODE AND THE RESULT IS FED INTO DATABASE FROM FEEDW*/
       IF NULLANS = Y THEN (
		    ZEROANS:READ("IS DIFF(W(AMP,",I,",LAM,",J,") IDENTICALLY 
					ZERO, Y/N"),
		    IF ZEROANS = Y THEN (
				 ADDBASE:DIFNULL(I,J),
				 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),
       DERIVSUBST:TRUE,
       /*RENAME THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZW */
       TEMP3:SUBST(WMAX_DIFF,TEMP2),
       /*ENTER ALL INFORMATION IN STORED IN DATABASE*/
       D_BASE_LENGTH:LENGTH(DATABASE),
       FOR II THRU D_BASE_LENGTH DO
                TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF),
       DERIVSUBST:FALSE,
       TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE),
       /*THIS IS THE PROJECTION Q ONTO THE RANGE OF THE
         LINEAR DIFFERENTIAL OPERATOR IN THE PROBLEM*/
       TEMP5:INTEGRATE(EV(TEMP4,ZWNULL).AFUN,XVAR,0,LEN),
       TEMP6:TEMP4-TEMP5/NORM*CFUN, 
       TEMP7:TRIGREDUCE(TEMP6),
       /*THE SET OF O.D.E.'S TO SOLVE */
       W_DE:EXPAND(TEMP7),
       TEMP8:EV(W_DE,VFUN(ZWLIST,0)),
       /*IF THE PARTICULAR SOLUTION OF W_DE IS ZERO THEN W=0 */
       IF NULLIST=TEMP8 THEN ( ADDBASE:DIFNULL(I,J),
			       DATABASE:APPEND(DATABASE,ADDBASE),
			       RETURN(ADDBASE)),
       TEMP9:SOLVE_ODE(TEMP8),
       FEEDW(TEMP9)
       )$
 

SOLVE_ODE(EXP):=(
/*THIS FUNCTION SOLVE THE D.E. W_DE BY A FOURIER MODE ANSATZ*/
       TRIGFUN:[],
       CONST:FALSE,
       FOR I THRU N DO (
       /*DETERMINE THE FOURIER MODES*/
		    TRIG1:EXP[I],
		    IF TRIG1 # 0 THEN  (
		                 TRIG2:APPLY1(TRIG1,SINNULL,COSNULL),
				 IF TRIG2 # 0 THEN (
			                      CONST:TRUE,
					      TRIG1:TRIG1-TRIG2),
		    TRIGFUN:APPEND(FIND_TRIG(TRIG1),TRIGFUN))),
       SOL1:DELETE(DUMMY,SETIFY(TRIGFUN)),
       /*MAKE AN ANSATZ*/
       ANSATZ:MAKELIST(SUM(AM[I,J]*SOL1[I],I,1,LENGTH(SOL1)),J,1,N),
       SOL2:EV(W_DE,MAP("=",ZWLIST,ANSATZ),DIFF),
       SOL3:MAKELIST(RATCOEF(SOL2,I),I,SOL1),
       EQLIST:[],
       FOR I THRU LENGTH(SOL3) DO EQLIST:APPEND(EQLIST,SOL3[I]),
       VARLIST:[],
       FOR I THRU N DO 
		FOR J THRU LENGTH(SOL1) DO VARLIST:CONS(AM[J,I],VARLIST),
       /*FIND THE AMPLITUDES OF THE FOURIER MODES*/
       SOL4:SOLVE(EQLIST,VARLIST),
       /*SOLVE FOR THE CONSTANT FOURIER MODE IF NECESSARY*/
       CANSATZ:0,
       IF CONST = TRUE THEN (CANSATZ:MAKELIST(CONCAT(CON,I),I,1,N),
			     CSOL1:EV(W_DE,MAP("=",ZWLIST,CANSATZ),DIFF),
			     CSOL2:APPLY1(CSOL1,SINNULL,COSNULL),
			     CSOL3:SOLVE(CSOL2,CANSATZ)),
       EV(ANSATZ+CANSATZ,SOL4,CSOL3))$


FEEDW(EXP):=(
/* THIS FUNCTION ALLOWS THE RESULT OF THE ODE-SOLVER TO BE ENTERED INTO 
   THE LIST DATABASE. IT CHECKS FOR ORTHOGONALITY TO THE CRITICAL ADJOINT
   EIGENFUNCTION AND REMOVES  SOLUTIONS OF THE HOMOGENEOUS EQUATION (I.E. 
   NONORTHOGONAL TERMS)*/
       F2:INTEGRATE(EXP.AFUN,XVAR,0,LEN),
       IF RATSIMP(F2)=0 
	        THEN 
                     ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),
							EXP)
	        ELSE (ORTHO_RESULT:RATSIMP(EXP- F2/NORM*CFUN),
		     ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),
						ORTHO_RESULT)),
       DATABASE:APPEND(DATABASE,ADDBASE),
       PRINT(ADDBASE))$
/*COLLECT ALL INFORMATION STORED IN BIFEQ AND ASSEMBLE THE BIFURCATION
EQUATION*/
G_EQ():=  SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$


SETIFY(LIST):=(/*TRANSFORMS A LIST INTO A SET, I.E. REMOVES DUPLICATES*/
	SET:[LIST[1]],
        FOR I :2 THRU LENGTH(LIST) DO
	       IF NOT MEMBER(LIST[I],SET) THEN
			SET:CONS(LIST[I],SET),
	SET)$

FIND_TRIG(EXP):=(/*FINDS THE FOURIER MODES*/
	F_A1:ARGS(EXP+DUMMY),
	F_A2:APPLY1(F_A1,SINFIND,COSFIND)
	)$

/*AUXILIARY FUNCTIONS */
MATCHDECLARE([DUMMY1,DUMMY2],TRUE)$
DEFRULE(COSFIND,DUMMY1*COS(DUMMY2),COS(DUMMY2))$
DEFRULE(SINFIND,DUMMY1*SIN(DUMMY2),SIN(DUMMY2))$
DEFRULE(COSNULL,COS(DUMMY1),0)$
DEFRULE(SINNULL,SIN(DUMMY1),0)$


VFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$
DIFNULL(I,J):=MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J)=0),W) $