/*PROGRAM NUMBER 10: REDUCTION1(), A LIAPUNOV-SCHMIDT REDUCTION FOR
  STEADY STATE BIFURCATIONS IN ONE DIFFERENTIAL EQUATION DEPENDING
  ON ONE INDEPENDENT VARIABLE. SEE PAGE 170 IN "PERTURBATION METHODS,
  BIFURCATION THEORY AND COMPUTER ALGEBRA". */


  

/*THIS FILE CONTAINS REDUCTION1(), A FUNCTION TO PERFORM A LIAPUNOV-
  SCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS OF NONLINEAR D.E.'S
  OF THE FORM Y'' + F(Y,Y',ALPHA) = 0. Y = Y(X) IS DEFINED ON A REAL
  INTERVAL WITH DIRICHLET OR NEUMANN BOUNDARY CONDITIONS AND F DEPENDS 
  ONLY LINEARLY ON ALPHA.
  IT ALSO CONTAINS THESE ADDITIONAL FUNCTIONS:
  SETUP() ALLOWS TO ENTER THE PROBLEM.
  G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION
  EQUATION.
  W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM).
  PROJECT(EXP) ENSURES THAT <CFUN,EXP>=0.
  NEUMANNBC(EXP) ACCOUNTS FOR NEUMANN BOUNDARY CONDITIONS.
  G_EQ() ASSEMBLES THE BIFURCATION EQUATION. */


REDUCTION1():=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 ASKS FOR THE VARIABLES OF THE D.E., THE BIFURCATION
POINT, THE CRITICAL EIGENFUNCTION, THE X-INTERVAL, THE BOUNDARY
CONDITIONS AND THE DIFFERENTIAL EQUATION. */
ASSUME_POS:TRUE,
LS_LIST:[],
Y:READ("ENTER DEPENDENT VARIABLE"),
PRINT("USE X AS THE INDEPENDENT VARIABLE AND ALPHA AS A PARAMETER TO VARY"),
CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE ALPHA"),
PRINT("WE DEFINE LAM = ALPHA - ",CAL),
CFUN: READ("ENTER THE CRITICAL EIGENFUNCTION"),
DEPENDS([ZW,Y,W],[AMP,LAM,X]),
LEN:READ("WHAT IS THE LENGTH OF THE X-INTERVAL"),
NORM:INTEGRATE(CFUN^2,X,0,LEN),
PRINT("SPECIFY THE BOUNDARY CONDITIONS"),
PRINT("YOUR CHOICE FOR THE B.C. ON Y AT X=0 AND X=",LEN),
PRINT("ENTER 1 FOR Y=0, 2 FOR Y'=0"),
BCL:READ("B.C. AT 0?"),
BCR:READ("B.C. AT",LEN,"?"),
EQ:DIFF(Y,X,2) +
	 READ("THE D.E. IS OF THE FORM Y'' + F(Y,Y',ALPHA) = 0,ENTER F"),
EQLAM:EV(EQ,ALPHA=LAM+CAL),
PRINT(EQLAM),
DATABASE:[DIFF(W,AMP)=0,DIFF(W,LAM)=0],
SUB:Y=AMP*CFUN+W,
TEMP1:EV(EQLAM,SUB,DIFF),
NULLANS:READ("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS ARE ZERO, 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 KNOWLEDGE 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,
/*THIS DIFFERENTIAL OF W WILL BE ANNIHILATED BY THE PROJECTION ONTO THE
CRITICAL EIGENFUNCTION, HENCE WE SET IT TO ZERO ALREADY HERE */
	  TEMP3:SUBST('DIFF(W,AMP,I,LAM,J)=0,TEMP2),
/* HERE WE INSERT ALL KNOWLEDGE GAINED THROUGH W_POLY */
	  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,W=0),
/*PROJECTION ONTO CFUN, THE CRITICAL EIGENFUNCTION */
          TEMP5:TRIGREDUCE(CFUN*TEMP4),
	  BIFEQ[I,J]:RATSIMP(INTEGRATE(TEMP5,X,0,LEN))
	  )$


  

W_POLY(I,J):=BLOCK( 
/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR
COEFFICIENT OF THE FUNCTION W(AMP,LAM). THE RESULT IS FED INTO DATABASE.*/
	  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))),
	  TEMP2:DIFF(TEMP1,AMP,I,LAM,J),
	  DERIVSUBST:TRUE,
/*WE SUBSTITUTE ZW FOR THE UNKNOWN TAYLOR COEFFICIENT AND SOLVE
AN O.D.E. FOR ZW IN TEMP7 */
	  TEMP3:SUBST(DIFF(W,AMP,I,LAM,J)=ZW,TEMP2),
/*NOW WE INSERT ALL KNOWLEDGE GAINED THROUGH PREVIOUS ITERATIONS*/
	  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,W=0),
          TEMP5:TRIGREDUCE(TEMP4),
/*THIS ENSURES THAT THE R.H.S. FO THE D.E. TEMP6 IS ORTHOGONAL TO
THE SOLUTION OF THE HOMOGENEOUS EQUATION */
	  TEMP6:PROJECT(TEMP5),
	  TEMP7:ODE2(TEMP6,ZW,X),
/*SATISFY BOUNDARY CONDITIONS*/
	  IF BCL*BCR=1 THEN TEMP8:BC2(TEMP7,X=0,ZW=0,X=LEN,ZW=0)
	                ELSE
			     TEMP8:NEUMANNBC(TEMP7),
/*MAKE SURE THAT <CFUN,W>= 0*/
	  TEMP9:PROJECT(TEMP8),
	  ADDBASE:['DIFF(W,AMP,I,LAM,J)=RHS(TEMP9)],
	  DATABASE:APPEND(DATABASE,ADDBASE),
	  PRINT(ADDBASE))$


PROJECT(EXP):=(
	  PRO1:EV(EXP,ZW=0),
	  PRO2:INTEGRATE(PRO1*CFUN,X,0,LEN)/NORM,
	  EXPAND(EXP-PRO2*CFUN))$
NEUMANNBC(EXP):=(
	  REXP:RHS(EXP),
	  NBC1:DIFF(REXP,X),
	  IF BCL=1 AND BCR=2 THEN NBC2:SOLVE([EV(REXP,X=0),EV(NBC1,X=LEN)],
		[%K1,%K2]),
	  IF BCL=2 AND BCR=1 THEN NBC2:SOLVE([EV(REXP,X=LEN),EV(NBC1,X=0)],
		[%K1,%K2]),
	  IF BCL=2 AND BCR=2 THEN NBC2:SOLVE([EV(NBC1,X=LEN),EV(NBC1,X=0)],
		[%K1,%K2]),
	  EV(EXP,NBC2))$
G_EQ():=  SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$