#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'program1' <<'END_OF_FILE' X X/*PROGRAM NUMBER 1: LC(), LINDSTEDT'S METHOD TO CALCULATE LIMIT CYCLES, X SEE PAGE 7 IN "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER X ALGEBRA".*/ X X/* THIS PROGRAM APPLIES LINDSTEDT'S METHOD TO XTHE EQUATION: X X'' + X + E F(X,X') = 0, XASSUMING A LIMIT CYCLE EXISTS. XCALL IT WITH: X LC(); X*/ X XLC():=( X X/* INPUT THE DIFFERENTIAL EQUATION */ XKILL(X,XLIST,PARAMLIST), XPRINT("THE D.E. IS OF THE FORM: X'' + X + E * F(X,X') = 0"), XF:READ("ENTER F(X,Y), REPRESENTING X' AS Y"), XPRINT("THE D.E. IS: X'' + X + E (",F,") = 0"), XF:SUBST('DIFF(X,Z,1)*W,Y,F), X X/* SET UP THE SERIES EXPANSIONS */ XN:READ("ENTER TRUNCATION ORDER"), XW:1, XFOR I THRU N DO W:W+K[I]*E^I, XX:B[0]*COS(Z), XXLIST:[XX[0] = X], XFOR I THRU N DO X:X+XX[I]*E^I, X X/* PLUG INTO THE D.E. AND COLLECT TERMS */ XDEPENDS(XX,Z), XTEMP1:DIFF(X,Z,2)+X/W^2+E*EV(F,DIFF)/W^2, XTEMP1:TAYLOR(TEMP1,E,0,N), XFOR I THRU N DO EQ[I]:COEFF(TEMP1,E,I), X X/* SET UP PATTERN MATCHING RULES FOR USE IN SOLVING D.E. */ XMATCHDECLARE(N1,TRUE), XDEFRULE(C,COS(N1*Z),COS(N1*Z)/(N1*N1-1)), XDEFRULE(S,SIN(N1*Z),SIN(N1*Z)/(N1*N1-1)), X X/* LOAD POISSON SERIES PACKAGE AND SET PARAMETER */ XOUTOFPOIS(DUMMY), XPOISLIM:100, X X/* MAIN LOOP */ XFOR I:1 THRU N DO BLOCK( X X/* TRIGONOMETRIC SIMPLIFICATION */ XTEMP1:OUTOFPOIS(EV(EQ[I],XLIST,PARAMLIST,DIFF)), X X/* ELIMINATE SECULAR TERMS */ XIF I = 1 X THEN (PARAMLIST:SOLVE(COEFF(TEMP1,SIN(Z)),B[0]), X PRINT("CHOICES FOR LIMIT CYCLE AMPLITUDE:"), X FOR J:1 THRU LENGTH(PARAMLIST) DO X PRINT(J,") ",PART(PARAMLIST,J,2)), X R1:READ("ENTER CHOICE NUMBER"), X PARAMLIST:APPEND(SOLVE(COEFF(TEMP1,COS(Z)),K[1]), X [PART(PARAMLIST,R1)])) X ELSE PARAMLIST:APPEND(PARAMLIST,SOLVE([COEFF(TEMP1,COS(Z)), X COEFF(TEMP1,SIN(Z))],[K[I],B[I-1]])), XTEMP1:EXPAND(EV(TEMP1,PARAMLIST)), XXLIST:EXPAND(EV(XLIST,PARAMLIST)), X X/* OUTPUT PROGRESS */ XPRINT("DONE WITH STEP",I,"OF",N), X X/* EXIT HERE IF LAST ITERATION */ XIF I=N THEN GO(END), X X/* SOLVE THE D.E. */ XTEMP1:FACTOR(EV(TEMP1,XX[I] = 0)), XTEMP1:APPLYB1(TEMP1,C,S), XTEMP1:XX[I] = TEMP1+A[I]*SIN(Z)+B[I]*COS(Z), X X/* FIT THE INITIAL CONDITION */ XTEMP2:RHS(TEMP1), XTEMP2:DIFF(TEMP2,Z), XTEMP2:SOLVE(EV(TEMP2,Z:0),A[I]), XXLIST:APPEND(XLIST,[EV(TEMP1,TEMP2)]), X X/* END OF MAIN LOOP */ XEND), X X/* OUTPUT RESULTS */ XW:EV(W,PARAMLIST), XX:TAYLOR(EV(X,XLIST,PARAMLIST),E,0,N-1), XPRINT("X =",X), XPRINT("W =",W))$ X X X X X X END_OF_FILE if test 2446 -ne `wc -c <'program1'`; then echo shar: \"'program1'\" unpacked with wrong size! fi # end of 'program1' fi if test -f 'program2' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program2'\" else echo shar: Extracting \"'program2'\" \(1880 characters\) sed "s/^X//" >'program2' <<'END_OF_FILE' X X X X X X XLC2():=( X XKILL(X,Y,XYLIST,PARAMLIST), XPRINT("THE D.E.'S ARE OF THE FORM X' = -Y + E*F, Y' = X + E*G"), XPRINT("WHERE F,G ARE OF THE FORM: X QUADRATIC + E*(LINEAR + CUBIC) + E^2*(QUARTIC) + ..."), XF:READ("ENTER F:"), XPRINT("F =",F), XG:READ("ENTER G:"), XPRINT("G =",G), X XN:READ("ENTER TRUNCATION ORDER:"), XK[0]:1, XK[1]:0, XW:SUM(K[I]*E^I,I,0,N), XXY:[X:B[0]*COS(Z)+SUM(XX[I](Z)*E^I,I,1,N), X Y:B[0]*SIN(Z)+SUM(YY[I](Z)*E^I,I,1,N)], XXYLIST:[XX[0](Z)=B[0]*COS(Z), X YY[0](Z)=B[0]*SIN(Z)], X XTEMP1:[-DIFF(X,Z)*W-Y+E*EV(F,DIFF),-DIFF(Y,Z)*W+X+E*EV(G,DIFF)], XTEMP2:TAYLOR(TEMP1,E,0,N), XFOR I:1 THRU N DO EQ[I]:COEFF(TEMP2,E,I), X X XFOR I:1 THRU N DO BLOCK( X XTEMP3:EXPAND(TRIGREDUCE(EXPAND(EV(EQ[I],XYLIST,PARAMLIST,DIFF)))), X XIF I=1 THEN (TEMP4:TEMP3, GO(SKIP_TO_HERE_FIRST_TIME)) X ELSE NEWPARAMLIST: X SOLVE([COEFF(PART(TEMP3,1),SIN(Z))-COEFF(PART(TEMP3,2),COS(Z)), X COEFF(PART(TEMP3,1),COS(Z))+COEFF(PART(TEMP3,2),SIN(Z))], X [B[I-2],K[I]]), XIF I=2 THEN (PARAMLIST:NEWPARAMLIST, X PRINT("CHOICES FOR LIMIT CYCLE AMPLITUDE:"), X FOR J:1 THRU LENGTH(PARAMLIST) DO X PRINT(J,") ",PART(PARAMLIST,J,1,2)), X R1:READ("ENTER CHOICE NUMBER"), X PARAMLIST:PART(PARAMLIST,R1)) X ELSE PARAMLIST:APPEND(PARAMLIST,NEWPARAMLIST), XTEMP4:EXPAND(EV(TEMP3,PARAMLIST)), XXYLIST:EXPAND(EV(XYLIST,PARAMLIST)), XSKIP_TO_HERE_FIRST_TIME, X XPRINT("DONE WITH STEP",I,"OF",N), X XIF I=N THEN GO(END), X XTEMP4A:SUBST(DUMMY(Z),YY[I](Z),TEMP4), XATVALUE(DUMMY(Z),Z=0,0), XTEMP5:DESOLVE(TEMP4A,[XX[I](Z),DUMMY(Z)]), XTEMP5A:SUBST(YY[I](Z),DUMMY(Z),TEMP5), XTEMP5B:SUBST(B[I],XX[I](0),TEMP5A), XXYLIST:APPEND(XYLIST,[TEMP5B]), X XEND), X XW:EV(W,PARAMLIST), XSOLN:TAYLOR(EV([X,Y],XYLIST,PARAMLIST),E,0,N-2), XPRINT("X =",PART(SOLN,1)), XPRINT("Y =",PART(SOLN,2)), XPRINT("W =",W))$ END_OF_FILE if test 1880 -ne `wc -c <'program2'`; then echo shar: \"'program2'\" unpacked with wrong size! fi # end of 'program2' fi if test -f 'program3' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program3'\" else echo shar: Extracting \"'program3'\" \(2208 characters\) sed "s/^X//" >'program3' <<'END_OF_FILE' X X/* PROGRAM NUMBER 3: CM(), CENTER MANIFOLD REDUCTION FOR ORDINARY X DIFFERENTIAL EQUATIONS. SEE PAGE 32 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X XCM():=( X X/* INPUT PROBLEM */ XN:READ("ENTER NO. OF EQS."), XM:READ("ENTER DIMENSION OF CENTER MANIFOLD"), XPRINT("THE D.E.'S MUST BE ARRANGED SO THAT THE FIRST",M,"EQS."), XPRINT("REPRESENT THE CENTER MANIFOLD. I.E. ALL ASSOCIATED"), XPRINT("EIGENVALUES ARE ZERO OR HAVE ZERO REAL PARTS."), XFOR I:1 THRU N DO X X[I]:READ("ENTER SYMBOL FOR VARIABLE NO.",I), XL:READ("ENTER ORDER OF TRUNCATION"), XFOR I:1 THRU N DO ( X PRINT("ENTER RHS OF EQ.",I), X PRINT("D",X[I],"/DT ="), X G[I]:READ()), X/* SET UP D.E.'S */ XFOR I:1 THRU N DO X DEPENDS(X[I],T), XFOR I:1 THRU N DO X (EQ[I]:DIFF(X[I],T)=G[I], X PRINT(EQ[I])), X X/* FORM POWER SERIES */ XSUB:MAKELIST(K[I],I,1,M), XVAR:PRODUCT(X[I]^K[I],I,1,M), XUNK:[], XFOR P:M+1 THRU N DO( X TEMP:A[P,SUB]*VAR, X FOR I:1 THRU M DO X TEMP:SUM(EV(TEMP,K[I]=J),J,0,L), X TEMP2:TAYLOR(TEMP,MAKELIST(X[I],I,1,M),0,L), X /* REMOVE CONSTANT AND LINEAR TERMS */ X TEMP3:TEMP2-PART(TEMP2,1)-PART(TEMP2,2), X SOLN[P]:EXPAND(TEMP3), X /* PREPARE LIST OF UNKNOWNS */ X SETXTO1:MAKELIST(X[I]=1,I,1,M), X /* TURN SUM INTO A LIST */ X UNKN[P]:SUBST("[","+",EV(TEMP3,SETXTO1)), X UNK:APPEND(UNK,UNKN[P])), XSOL:MAKELIST(X[P]=SOLN[P],P,M+1,N), X X/* SUBSTITUTE INTO D.E.'S */ XCMDE:MAKELIST(EQ[I],I,1,M), XREST:MAKELIST(LHS(EQ[I])-RHS(EQ[I]),I,M+1,N), XTEMP4:EV(REST,SOL,DIFF), XTEMP5:EV(TEMP4,CMDE,DIFF), XTEMP6:EV(TEMP5,SOL), XTEMP7:TAYLOR(TEMP6,MAKELIST(X[I],I,1,M),0,L), X X/* COLLECT TERMS */ XCOUNTER:1, X/* MAKE LIST OF TERMS */ XTERMS:SUBST("[","+",SOLN[N]), XTERMS:EV(TERMS,A[DUMMY,SUB]:=1), XFOR I:1 THRU N-M DO ( X EXP[I]:EXPAND(PART(TEMP7,I)), X FOR J:1 THRU LENGTH(TERMS) DO( X COND[COUNTER]:RATCOEF(EXP[I],PART(TERMS,J)), X COUNTER:COUNTER+1)), XCONDS:MAKELIST(COND[I],I,1,COUNTER-1), XCONDS:EV(CONDS,MAKELIST(X[I]=0,I,1,M)), X X/* SOLVE FOR CENTER MANIFOLD */ XACOEFFS:SOLVE(CONDS,UNK), XCENTERMANIFOLD:EV(SOL,ACOEFFS), XPRINT("CENTER MANIFOLD:"), XPRINT(CENTERMANIFOLD), X X/* GET FLOW ON CM */ XCMDE2:EV(CMDE,CENTERMANIFOLD), XPRINT("FLOW ON THE C.M.:"), XPRINT(CMDE2))$ X X X X X X X X END_OF_FILE if test 2208 -ne `wc -c <'program3'`; then echo shar: \"'program3'\" unpacked with wrong size! fi # end of 'program3' fi if test -f 'program4' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program4'\" else echo shar: Extracting \"'program4'\" \(1104 characters\) sed "s/^X//" >'program4' <<'END_OF_FILE' X/*PROGRAM NUMBER 4: TRANSFORM(), THIS PROGRAM ALLOWS TO PERFORM ARBITRARY X COORDINATE TRANSFORMATIONS. SEE PAGE 44 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X XTRANSFORM():=( X X/* INPUT DATA */ XN:READ("ENTER NUMBER OF EQUATIONS"), XFOR I:1 THRU N DO X X[I]:READ("ENTER SYMBOL FOR ORIGINAL VARIABLE",I), XFOR I:1 THRU N DO X Y[I]:READ("ENTER SYMBOL FOR TRANSFORMED VARIABLE",I), XPRINT("THE RHS'S OF THE D.E.'S ARE FUNCTIONS OF THE ORIGINAL VARIABLES:"), XFOR I:1 THRU N DO ( X PRINT("ENTER RHS OF",X[I],"D.E."), X PRINT("D",X[I],"/DT ="), X F[I]:READ(), X PRINT("D",X[I],"/DT =",F[I])), XPRINT("THE TRANSFORMATION IS ENTERED NEXT:"), XFOR I:1 THRU N DO ( X PRINT("ENTER",X[I],"AS A FUNCTION OF THE NEW VARIABLES"), X PRINT(X[I],"="), X G[I]:READ(), X PRINT(X[I],"=",G[I])), X X/* DO IT */ XFOR I:1 THRU N DO DEPENDS([X[I],Y[I]],T), XFOR I:1 THRU N DO EQ[I]:DIFF(X[I],T)=F[I], XTRANS:MAKELIST(X[I]=G[I],I,1,N), XFOR I:1 THRU N DO TREQ[I]:EV(EQ[I],TRANS,DIFF), XTREQS:MAKELIST(TREQ[I],I,1,N), XDERIVS:MAKELIST(DIFF(Y[I],T),I,1,N), XNEWEQS:SOLVE(TREQS,DERIVS))$ X X X X X X X X END_OF_FILE if test 1104 -ne `wc -c <'program4'`; then echo shar: \"'program4'\" unpacked with wrong size! fi # end of 'program4' fi if test -f 'program5' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program5'\" else echo shar: Extracting \"'program5'\" \(4943 characters\) sed "s/^X//" >'program5' <<'END_OF_FILE' X X/*PROGRAM NUMBER 5: NF(), NORMAL FORM TRANSFORMATIONS. SEE PAGE 62 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X/* THIS FILE CONTAINS NF(), A NORMAL FORM UTILITY FUNCTION. IT ALSO CONTAINS X THESE ADDITIONAL FUNCTIONS: X GEN(N) WILL GENERATE A HOMOGENEOUS ORDER N TRANSFORMATION. X DECOMPOSE() ISOLATES THE COEFFICIENTS OF THE NEW EQUATIONS. X VARS(N) GENERATES A LIST OF UNKNOWN COEFFICIENTS OF DEGREE N. X HOPFK(), FOR K=2,3,4,5,6,7 SOLVES FOR THE COEFFICIENTS OF A SYSTEM OF X 2 DE'S SO AS TO PUT THE EQS IN HOPF NORMAL FORM */ X XNF():= BLOCK( X X/* NEW VARIABLE NAMES? */ XTEST : READ ("DO YOU WANT TO ENTER NEW VARIABLE NAMES (Y/N)?"), XIF TEST = N THEN GO(JUMP), XN : READ ("HOW MANY EQS"), XFOR I THRU N DO (X[I] : READ ("SYMBOL FOR OLD X[",I,"]")), XFOR I THRU N DO( Y[I] : READ ("SYMBOL FOR NEW X[",I,"]")), XFOR I THRU N DO DEPENDS([X[I],Y[I]],T), XKILL(FLAG), /* FLAG USED IN GEN */ X XJUMP, X X/* NEW D.E.'S? */ XPRINT ("DO YOU WANT TO ENTER NEW D.E.'S (Y/N)?"), XTEST:READ(), XIF TEST = N THEN GO(LOOP), XFOR I THRU N DO X (RHS[I]:READ("ENTER RHS OF EQ. NO.",I,",D",X[I],"/DT ="), X PRINT("D",X[I],"/DT =",RHS[I])), XKILL(RHS2), XRHS2[I,J] := RHS[I], XRHS3:GENMATRIX(RHS2,N,1), X XLOOP, X X/* NEAR-IDENTITY TRANSFORMATION */ XPRINT("INPUT NEAR-IDENTITY TRANSFORMATION X(USE PREV FOR PREVIOUS TRANSFORMATION)"), XFOR I THRU N DO X (ROW:I, X PREV :TR[I], X TR[I] :READ (X[I],"=",Y[I],"+ ?"), X PRINT (X[I],"=",Y[I]+TR[I])), X X/* INPUT TRUNCATION ORDER */ XTRANS : MAKELIST (X[I]=Y[I]+TR[I],I,1,N), XM : READ("ENTER TRUNCATION ORDER (HIGHEST ORDER TERMS TO BE KEPT)"), X X X/* TRANSFORM THE D.E.'S */ XTEMP2 :EV(RHS3,TRANS), X X/* SOLVE FOR THE TRANSFORMED DERIVATIVES */ XKILL(JACOB), XJACOB[I,J]:=DIFF(TR[I],Y[J]), XJACOB2:GENMATRIX(JACOB,N,N), XTEMP3:SUM((-1)^I*JACOB2^^I,I,0,M-1).TEMP2, X X/* TAYLOR EXPAND THE RESULTING EQS */ XNEWRHS : TAYLOR(TEMP3,MAKELIST(Y[I],I,1,N),0,M), XNEWDES:MAKELIST(DIFF(Y[I],T)=NEWRHS[I,1],I,1,N), XFOR I:1 THRU N DO X PRINT(PART(NEWDES,I)), X X/* ENTER ANOTHER TRANSFORMATION? */ XBRANCH:READ("DO YOU WANT TO ENTER ANOTHER TRANSFORMATION (Y/N)"), XIF BRANCH = Y THEN GO(LOOP), XNEWDES)$ X X XGEN(NN):=( XIF NOT NUMBERP(FLAG[NN]) THEN ( X SUB:MAKELIST(K[I],I,1,N), X VAR:PRODUCT(Y[I]^K[I],I,1,N), X TEMPGEN:A[ROWDUMMY,SUB]*VAR, X FOR I:1 THRU N DO X TEMPGEN:SUM(EV(TEMPGEN,K[I]=J),J,0,NN), X TEMPGEN2:LAST(TAYLOR(TEMPGEN,MAKELIST(Y[I],I,1,N),0,NN)), X TEMPGEN3[NN]:EXPAND(TEMPGEN2), X FLAG[NN] : 1), XEV(TEMPGEN3[NN],ROWDUMMY=ROW)) $ X XDECOMPOSE():=( XKILL(C), XIF NOT NUMBERP(FLAG[M]) THEN GEN(M), XTEMP8:SUBST("[","+",TEMPGEN3[M]), XTERMS:EV(TEMP8,A[DUMMY,SUB]:=1), XCOEFFS:EV(TEMP8,A[DUMMY,SUB]:=C[DUMMY,SUB],MAKELIST(Y[I]=1,I,1,N)), XFOR I:1 THRU N DO( X FOR J:1 THRU LENGTH(TERMS) DO( X EV(PART(COEFFS,J),ROWDUMMY=I)::RATCOEF(EXPAND(NEWRHS[I,1]),PART(TERMS, X J)))))$ X XVARS(NN):=( XTEMP5:SUM(EV(TEMPGEN3[NN]),ROWDUMMY,1,N), XTEMP6:SUBST("[","+",TEMP5), XTEMP7:EV(TEMP6,MAKELIST(Y[I]=1,I,1,N)))$ X X X HOPF2():=(DECOMPOSE(), X SOLVE([C[1,[2,0]],C[1,[1,1]],C[1,[0,2]],C[2,[2,0]], X C[2,[1,1]],C[2,[0,2]]], X VARS(2)))$ X X HOPF3():=(DECOMPOSE(), X SOLVE([C[1,[3,0]]=C[1,[1,2]],C[1,[3,0]]=C[2,[2,1]], X C[1,[3,0]]=C[2,[0,3]],C[1,[0,3]]=C[1,[2,1]], X C[1,[0,3]]=-C[2,[3,0]],C[1,[0,3]]=-C[2,[1,2]]], X VARS(3)))$ X X HOPF4():=(DECOMPOSE(), X SOLVE([C[1,[4,0]],C[1,[3,1]],C[1,[2,2]],C[1,[1,3]], X C[1,[0,4]],C[2,[4,0]],C[2,[3,1]],C[2,[2,2]], X C[2,[1,3]],C[2,[0,4]]], X VARS(4)))$ X X HOPF5():=(DECOMPOSE(), X SOLVE([C[1,[5,0]]=C[1,[3,2]]/2,C[1,[5,0]]=C[1,[1,4]], X C[1,[5,0]]=C[2,[4,1]],C[1,[5,0]]=C[2,[2,3]]/2, X C[1,[5,0]]=C[2,[0,5]],C[2,[5,0]]=C[2,[3,2]]/2, X C[2,[5,0]]=C[2,[1,4]],C[2,[5,0]]=-C[1,[4,1]], X C[2,[5,0]]=-C[1,[2,3]]/2,C[2,[5,0]]=-C[1,[0,5]]], X VARS(5)))$ X X HOPF6():=(DECOMPOSE(), X SOLVE([C[1,[6,0]],C[1,[5,1]],C[1,[4,2]],C[1,[3,3]], X C[1,[2,4]],C[1,[1,5]],C[1,[0,6]],C[2,[6,0]], X C[2,[5,1]],C[2,[4,2]],C[2,[3,3]],C[2,[2,4]], X C[2,[1,5]],C[2,[0,6]]], X VARS(6)))$ X X HOPF7():=(DECOMPOSE(), X SOLVE([C[1,[7,0]]=C[1,[5,2]]/3,C[1,[7,0]]=C[1,[3,4]]/3, X C[1,[7,0]]=C[1,[1,6]],C[1,[7,0]]=C[2,[6,1]], X C[1,[7,0]]=C[2,[4,3]]/3,C[1,[7,0]]=C[2,[2,5]]/3, X C[1,[7,0]]=C[2,[0,7]],C[2,[7,0]]=C[2,[5,2]]/3, X C[2,[7,0]]=C[2,[3,4]]/3,C[2,[7,0]]=C[2,[1,6]], X C[2,[7,0]]=-C[1,[6,1]],C[2,[7,0]]=-C[1,[4,3]]/3, X C[2,[7,0]]=-C[1,[2,5]]/3,C[2,[7,0]]=-C[1,[0,7]]], X VARS(7)))$ X X X X X X X END_OF_FILE if test 4943 -ne `wc -c <'program5'`; then echo shar: \"'program5'\" unpacked with wrong size! fi # end of 'program5' fi if test -f 'program6' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program6'\" else echo shar: Extracting \"'program6'\" \(4046 characters\) sed "s/^X//" >'program6' <<'END_OF_FILE' X/* PROGRAM NUMBER 6: TWOVAR(), A FUNCTION TO PERFORM A TWO VARIABLE X EXPANSION TO ORDER EPSILON FOR N WEAKLY COUPLED PERTURBED NONLINEAR X OSCILLATORS. SEE PAGE 94 IN "PERTURBATION METHODS, BIFURCATION X THEORY AND COMPUTER ALGEBRA". */ X X X XTWOVAR():=BLOCK( X XCHOICE:READ("DO YOU WANT TO ENTER NEW DATA (Y/N)"), X XIF CHOICE = N THEN GO(JUMP1), X X/* INPUT D.E.'S */ XN:READ("NUMBER OF D.E.'S"), X XPRINT("THE",N,"D.E.'S WILL BE IN THE FORM:"), XPRINT("X[I]'' + W[I]^2 X[I] = E F[I](X[1],...,X[",N,"],T)"), X XFOR I:1 THRU N DO X X[I]:READ("ENTER SYMBOL FOR X[",I,"]"), X XFOR I:1 THRU N DO X DEPENDS(X[I],T), X XFOR I:1 THRU N DO X W[I]:READ("ENTER W[",I,"]"), X XFOR I:1 THRU N DO X F[I]:READ("ENTER F[",I,"]"), X XJUMP1, X X/* UPDATE EQS FOR SUBSTITUTION OF RESONANT VALUES ON 2ND TIME THRU */ XFOR I:1 THRU N DO( X W[I]:EV(W[I]), X F[I]:EV(F[I])), X X/* ECHO BACK THE D.E.'S */ XPRINT("THE D.E.'S ARE ENTERED AS:"), XFOR I:1 THRU N DO X PRINT(X[I],"'' +",W[I]^2*X[I],"=",E*F[I]), X XPRINT("THE METHOD ASSUMES A SOLUTION IN THE FORM:"), XPRINT("X[I] = X0[I] + E X1[I]"), XPRINT("WHERE X0[I] = A[I](ETA) COS W[I] XI + B[I](ETA) SIN W[I] XI"), XPRINT("WHERE XI = T AND ETA = E T"), X X/* DON'T USE A OR B AS PARAMETERS IN THE GIVEN D.E.'S */ XDEPENDS([A,B],ETA), X XFOR I:1 THRU N DO X X0[I]:A[I]*COS(W[I]*XI)+B[I]*SIN(W[I]*XI), X XFOR I:1 THRU N DO X G[I]:EV(F[I],T=XI), X XFOR I:1 THRU N DO( X FOR J:1 THRU N DO X G[I]:EV(G[I],X[J]::X0[J])), X XFOR I:1 THRU N DO( X G[I]:G[I]-2*DIFF(X0[I],XI,1,ETA,1), X G[I]:EV(G[I],DIFF), X G[I]:EXPAND(TRIGREDUCE(EXPAND(G[I])))), X X/* COLLECT SECULAR TERMS */ XFOR I:1 THRU N DO( X S[I]:COEFF(G[I],SIN(W[I]*XI)), X C[I]:COEFF(G[I],COS(W[I]*XI))), X X/* DISPLAY SECULAR TERMS */ XPRINT("REMOVAL OF SECULAR TERMS IN THE X1[I] EQS. GIVES:"), X XFOR I:1 THRU N DO( X PRINT(S[I],"= 0"), X PRINT(C[I],"= 0")), X XABEQS:APPEND(MAKELIST(S[I],I,1,N),MAKELIST(C[I],I,1,N)), X XCHOICE2:READ("DO YOU WANT TO TRANSFORM TO POLAR COORDINATES (Y/N)"), X XIF CHOICE2=N THEN GO(JUMP2), X X/* TRANSFORM TO POLAR COORDINATES */ XDEPENDS([R,THETA],ETA), XTRANS:MAKELIST([A[I]=R[I]*COS(THETA[I]),B[I]=R[I]*SIN(THETA[I])],I,1,N), XINTEQS:EV(ABEQS,TRANS,DIFF), XRTHEQS:SOLVE(INTEQS,APPEND(MAKELIST(DIFF(R[I],ETA),I,1,N), X MAKELIST(DIFF(THETA[I],ETA),I,1,N))), XRTHEQS:TRIGSIMP(RTHEQS), XRTHEQS:EXPAND(TRIGREDUCE(RTHEQS)), XPRINT(RTHEQS), X XJUMP2, X XCHOICE3:READ("DO YOU WANT TO SEARCH FOR RESONANT PARAMETER VALUES (Y/N)"), X XIF CHOICE3=N THEN GO(END), X X/* OBTAIN FREQUENCIES WHICH APPEAR ON RHS'S OF D.E.'S */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE FREQS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR I:1 THRU N DO( X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(G[I]), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X TEMP3:APPLY1(TEMP2,COSARG,SINARG), X TEMP4:EV(TEMP3,XI=1), X/* REMOVE DUPLICATE ARGUMENTS */ X FREQ[I]:SETIFY(TEMP4)), X X/* DISPLAY FREQUENCIES */ XFOR I:1 THRU N DO( X PRINT(X[I],"EQ'S RESONANT FREQ =",W[I]), X PRINT("FREQS ON RHS =",FREQ[I])), X XJUMP3, X XPAR:READ("WHICH PARAMETER TO SEARCH FOR ?"), X X/* COMPUTE RESONANT VALUES OF DESIRED PARAMETER */ XRESVALS:[], XFOR I:1 THRU N DO( X FOR J:1 THRU LENGTH(FREQ[I]) DO( X RES:SOLVE(PART(FREQ[I],J)=W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES), X RES:SOLVE(PART(FREQ[I],J)=-W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES))), X X/* ELIMINATE DUPLICATE PARAMETER VALUES */ XRESVALUES:SETIFY(RESVALS), X X/* DISPLAY RESONANT PARAMETER VALUES */ XPRINT(RESVALUES), X XCHOICE4:READ("DO YOU WANT TO SEARCH FOR ANOTHER PARAMETER (Y/N) ?"), X XIF CHOICE4=Y THEN GO(JUMP3), X XEND," ")$ X X/* AUXILIARY FUNCTIONS */ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X END_OF_FILE if test 4046 -ne `wc -c <'program6'`; then echo shar: \"'program6'\" unpacked with wrong size! fi # end of 'program6' fi if test -f 'program7' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program7'\" else echo shar: Extracting \"'program7'\" \(3198 characters\) sed "s/^X//" >'program7' <<'END_OF_FILE' X/*PROGRAM NUMBER 7: AVERAGE(), ALLOWS TO PERFORM FIRST AND SECOND ORDER X AVERAGING. SEE PAGE 114 IN "PERTURBATION METHODS, BIFURCATION THEORY X AND COMPUTER ALGEBRA". */ X X X X/* PROGRAM TO PERFORM 1ST OR 2ND ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERAGE():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("DO YOU WANT FIRST OR SECOND ORDER AVERAGING? (1/2)"), X COEFVFEPS1:COEFF(VF2,EPS), X COEFVFEPS2:COEFF(VF2,EPS,2), X G1BAR:DEMOIVRE(APPLY1(COEFVFEPS1,KILLEXP)), X IF M=1 THEN RESULT:EPS*G1BAR X ELSE( X G1TILDE:COEFVFEPS1-G1BAR, X W1:INTEGRATE(G1TILDE,T), X REMARRAY(JACOB), X JACOB[I,J] := DIFF(G1TILDE[I],X[J]), X JACG1TILDE:GENMATRIX(JACOB,N,N), X G2BAR:DEMOIVRE(APPLY1(EXPAND(JACG1TILDE.W1)+COEFVFEPS2,KILLEXP)), X RESULT:EPS*G1BAR+EPS^2*G2BAR), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS TO KILL EXPONENTIAL TERMS */ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X END_OF_FILE if test 3198 -ne `wc -c <'program7'`; then echo shar: \"'program7'\" unpacked with wrong size! fi # end of 'program7' fi if test -f 'program8' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program8'\" else echo shar: Extracting \"'program8'\" \(4358 characters\) sed "s/^X//" >'program8' <<'END_OF_FILE' X X/* PROGRAM NUMBER 8: AVERM(), ALLOWS MTH ORDER AVERAGING. SEE PAGE 128 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X/* PROGRAM TO PERFORM MTH ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERM():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("ENTER ORDER OF AVERAGING"), X IF M=1 X/* FIRST ORDER AVERAGING */ X THEN (TEMP0:DEMOIVRE(APPLY1(F,KILLEXP)), X RESULT:TAYLOR(TEMP0,EPS,0,1)) X/* AVERAGING OF ORDER M > 1 */ X ELSE X Y:MAKELIST(CONCAT('Y,I),I,1,N), X WLIST:MAKELIST(CONCAT('W,I),I,1,N), X DEPENDS(WLIST,CONS(T,Y)), X TRAFO:Y, X XSUB:MAPLIST("=",Y,X), X/* WNULL SETS WLIST TO ZERO */ X WNULL:VFUN(WLIST,0), X JACOB[I,J] := DIFF(WLIST[I],Y[J]), X JAC:GENMATRIX(JACOB,N,N), X/* MAIN LOOP */ X FOR I :1 THRU M-1 DO ( X TEMP1:MAPLIST("=",X,Y+EPS^I*WLIST), X/* HERE X IS THE LIST [X1,X2,...,XN], Y IS THE LIST [Y1,Y2,...,YN], X WLIST IS THE TRANSFORMATION VECTOR [W1,W2,...,WN] */ X TEMP2:TAYLOR(SUM((-EPS)^(I*J)*JAC^^J,J,0,M-1). X (EV(VF2,TEMP1) - DIFF(WLIST,T)*EPS^I),EPS,0,M), X/* JAC IS THE JACOBIAN D WLIST/DY OF THE TRANSFORMATION WLIST */ X TEMP3:MATTOLIST(TEMP2,N), X TEMP4:MAP(EXPAND,TAYLOR(EV(TEMP3,WNULL,DIFF),EPS,0,I)), X/* THE ITH ORDER MEAN */ X MEAN:APPLY1(TEMP4,KILLEXP), X TEMP6:EXPAND(TEMP4-MEAN), X TEMP7:EV(INTEGRATE(TEMP6,T),EPS=1), X/* THE ITH ORDER TRANSFORMATION */ X TEMP8:MAPLIST("=",WLIST,TEMP7), X TEMP9:RATSIMP(TEMP8), X/* THE TRANSFORMED DE */ X VF2:EXPAND(EV(TEMP3,TEMP9,DIFF,XSUB,INFEVAL)), X/* THE SUM OF ALL TRANSFORMATIONS */ X TRAFO:EXPAND(TRAFO+EV(EPS^I*WLIST,TEMP9))), X/* END OF MAIN LOOP */ X PRINT("THE TRANSFORMATION: ",X,"="), X PRINT(RATSIMP(DEMOIVRE(TRAFO))), X/* THE FINAL AVERAGING */ X RESULT:APPLY1(VF2,KILLEXP), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS */ X VFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X MATTOLIST(MAT,DIM):=IF DIM>1 THEN MAKELIST(MAT[I,1],I,1,DIM) ELSE [MAT]$ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X X X X X X X END_OF_FILE if test 4358 -ne `wc -c <'program8'`; then echo shar: \"'program8'\" unpacked with wrong size! fi # end of 'program8' fi if test -f 'program9' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program9'\" else echo shar: Extracting \"'program9'\" \(4328 characters\) sed "s/^X//" >'program9' <<'END_OF_FILE' X X/* PROGRAM NUMBER 9: LIE(), LIE TRANSFORMATIONS FOR HAMILTONIAN SYSTEMS. X SEE PAGE 144 IN "PERTURBATION METHODS, BIFURCATION THEORY AND X COMPUTER ALGEBRA". */ X X XLIE():=BLOCK( X X/* INPUT PROBLEM ? */ XCHOICE1:READ("DO YOU WANT TO INPUT A NEW PROBLEM (Y/N) ?"), X XIF CHOICE1=N THEN GO(JUMP1), X X/* INPUT PROBLEM */ XN:READ("ENTER NUMBER OF DEGREES OF FREEDOM"), XFOR II:1 THRU N DO ( X Q[II]:READ("ENTER SYMBOL FOR Q[",II,"]"), X P[II]:READ("ENTER SYMBOL FOR P[",II,"]")), XKILL(W), XPRINT("THE HAMILTONIAN DEPENDS ON THE Q'S, P'S, T AND E (SMALL PARAMETER)"), XPRINT("THE E=0 PROBLEM MUST BE OF THE FORM:"), XPRINT("H =",SUM((P[II]^2+W[II]^2*Q[II]^2)/2,II,1,N)), XHORIGINAL:READ("ENTER THE HAMILTONIAN"), XPRINT("H =",HORIGINAL), X X/* TRANSFORM TO ACTION-ANGLE VARIABLES */ X/* FIND THE W[II]'S */ XH0:EV(HORIGINAL,E=0), XFOR II:1 THRU N DO X W[II]:SQRT(DIFF(H0,Q[II],2)), XPRINT("THE ACTION-ANGLE VARIABLES ARE J'S FOR ACTION, PHI'S FOR ANGLE"), XFOR II:1 THRU N DO X TR[II]:[Q[II]=SQRT(2*J[II]/W[II])*SIN(PHI[II]), X P[II]=SQRT(2*J[II]*W[II])*COS(PHI[II])], XTRAN:MAKELIST(TR[II],II,1,N), XH:EV(HORIGINAL,TRAN,ASSUME_POS:TRUE,INFEVAL), XH:TRIGSIMP(H), XH:EXPAND(TRIGREDUCE(EXPAND(H))), XPRINT("H =",H), X XJUMP1, X X/* INPUT TRUNCATION ORDER */ XNTRUNC:READ("ENTER HIGHEST ORDER TERMS IN E TO BE KEPT"), XFOR II:0 THRU NTRUNC DO X H[II]:RATCOEF(H,E,II), X X/* LIE TRANSFORMS */ X/* NEAR IDENTITY TRANSFORMATION FROM (J,PHI)'S TO (I,PSI)'S */ X X/* UPDATE VARIABLES */ XFOR II:1 THRU N DO( X P[II]:I[II], X Q[II]:PSI[II]), X X/* REPLACE J AND PHI BY I AND PSI IN H'S */ XFOR II:0 THRU NTRUNC DO X H[II]:EV(H[II],MAKELIST(J[III]=I[III],III,1,N), X MAKELIST(PHI[III]=PSI[III],III,1,N)), X XK[0]:H[0], X X/* DECLARE WGEN[I] TO BE A FN OF T, Q'S AND P'S */ XKILL(WGEN), XDEPENDS(WGEN,[T]), XFOR II:1 THRU N DO X DEPENDS(WGEN,[Q[II],P[II]]), X X/* E=0 PROBLEM IS OF FORM SUM(W[II]*I[II]) */ X/* CHOOSE WGEN[II] TO KILL AS MUCH AS POSSIBLE IN EQ(II) */ X/* EQUATE K[II] TO UNREMOVABLE TERMS */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE ARGS OF TRIG TERMS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR II:1 THRU NTRUNC DO( X X EQN[II]:EXPAND(TRIGREDUCE(EXPAND(EQ(II)))), X TEMP:EXPAND(EV(EQN[II],WGEN[II]=0)), X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(TEMP), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X ARG1:APPLY1(TEMP2,COSARG,SINARG), X/* REMOVE DUPLICATE ARGUMENTS */ X ARG2:SETIFY(ARG1), X/* REMOVE RESONANT ARGUMENTS */ X ARG3:SUBLIST(ARG2,NOTRESP), X X/* CHOOSE WGEN TO ELIMINATE NONRESONANT TERMS */ X LENG:LENGTH(ARG3), X WGENTEMP1:0, X FOR JJ:1 THRU LENG DO( X WGENTEMP2:AAA*COS(PART(ARG3,JJ))+BBB*SIN(PART(ARG3,JJ)), X TEMP4:EV(EQN[II],WGEN[II]=WGENTEMP2,DIFF), X TEMP5:SOLVE([RATCOEF(TEMP4,COS(PART(ARG3,JJ))), X RATCOEF(TEMP4,SIN(PART(ARG3,JJ)))],[AAA,BBB]), X WGENTEMP1:WGENTEMP1+EV(WGENTEMP2,TEMP5)), X WGEN[II]:WGENTEMP1, X PRINT("WGEN[",II,"] ="), X PRINT(WGEN[II]), X K[II]:EXPAND(EV(EQN[II],DIFF)), X K[II]:EXPAND(RATSIMP(K[II])), X PRINT("THE TRANSFORMED HAMILTONIAN K[",II,"] ="), X PRINT(K[II])), X XKAMILTONIAN:SUM(K[II]*E^II,II,0,NTRUNC), XPRINT ("THE TRANSFORMED HAMILTONIAN IS "), XPRINT ("K =",KAMILTONIAN), X XCHOICE2:READ("DO YOU WANT TO SEE THE NEAR IDENTITY TRANSFORMATION (Y/N) ?"), XIF CHOICE2=N THEN GO(END), X X/* THE NEAR IDENTITY TRANSFORMATION */ XFOR II:1 THRU N DO( X JTRANS[II]:SUM(S(III,P[II])*E^III,III,0,NTRUNC), X PHITRANS[II]:SUM(S(III,Q[II])*E^III,III,0,NTRUNC)), XFOR II:1 THRU N DO( X PRINT(J[II],"=",JTRANS[II]), X PRINT(PHI[II],"=",PHITRANS[II])), X XEND, XKAMILTONIAN)$ X X XPOISSON(F,G):= X SUM(DIFF(F,Q[II])*DIFF(G,P[II])-DIFF(F,P[II])*DIFF(G,Q[II]),II,1,N)$ X XL(II,F):=POISSON(WGEN[II],F)$ X XS(II,F):=(IF II=0 THEN F ELSE SUM(L(II-M,S(M,F)),M,0,II-1)/II)$ X XEQ(II):=(H[II]+(DIFF(WGEN[II],T)+POISSON(WGEN[II],H[0]))/II X +SUM(L(II-M,K[M])+M*S(II-M,H[M]),M,1,II-1)/II)$ X XLZAP(ANY):=DIFF(ANY,T)+POISSON(ANY,H[0])$ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XNOTRESP(EXP):=NOT IS(LZAP(EXP) = 0)$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X X X X X X X END_OF_FILE if test 4328 -ne `wc -c <'program9'`; then echo shar: \"'program9'\" unpacked with wrong size! fi # end of 'program9' fi if test -f 'program10' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program10'\" else echo shar: Extracting \"'program10'\" \(5930 characters\) sed "s/^X//" >'program10' <<'END_OF_FILE' X X/*PROGRAM NUMBER 10: REDUCTION1(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN ONE DIFFERENTIAL EQUATION DEPENDING X ON ONE INDEPENDENT VARIABLE. SEE PAGE 170 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X X X/*THIS FILE CONTAINS REDUCTION1(), A FUNCTION TO PERFORM A LIAPUNOV- X SCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS OF NONLINEAR D.E.'S X OF THE FORM Y'' + F(Y,Y',ALPHA) = 0. Y = Y(X) IS DEFINED ON A REAL X INTERVAL WITH DIRICHLET OR NEUMANN BOUNDARY CONDITIONS AND F DEPENDS X ONLY LINEARLY ON ALPHA. X IT ALSO CONTAINS THESE ADDITIONAL FUNCTIONS: X SETUP() ALLOWS TO ENTER THE PROBLEM. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION. X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X PROJECT(EXP) ENSURES THAT =0. X NEUMANNBC(EXP) ACCOUNTS FOR NEUMANN BOUNDARY CONDITIONS. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. */ X X XREDUCTION1():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X XSETUP():=( X/*THE FUNCTION SETUP ASKS FOR THE VARIABLES OF THE D.E., THE BIFURCATION XPOINT, THE CRITICAL EIGENFUNCTION, THE X-INTERVAL, THE BOUNDARY XCONDITIONS AND THE DIFFERENTIAL EQUATION. */ XASSUME_POS:TRUE, XLS_LIST:[], XY:READ("ENTER DEPENDENT VARIABLE"), XPRINT("USE X AS THE INDEPENDENT VARIABLE AND ALPHA AS A PARAMETER TO VARY"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE ALPHA"), XPRINT("WE DEFINE LAM = ALPHA - ",CAL), XCFUN: READ("ENTER THE CRITICAL EIGENFUNCTION"), XDEPENDS([ZW,Y,W],[AMP,LAM,X]), XLEN:READ("WHAT IS THE LENGTH OF THE X-INTERVAL"), XNORM:INTEGRATE(CFUN^2,X,0,LEN), XPRINT("SPECIFY THE BOUNDARY CONDITIONS"), XPRINT("YOUR CHOICE FOR THE B.C. ON Y AT X=0 AND X=",LEN), XPRINT("ENTER 1 FOR Y=0, 2 FOR Y'=0"), XBCL:READ("B.C. AT 0?"), XBCR:READ("B.C. AT",LEN,"?"), XEQ:DIFF(Y,X,2) + X READ("THE D.E. IS OF THE FORM Y'' + F(Y,Y',ALPHA) = 0,ENTER F"), XEQLAM:EV(EQ,ALPHA=LAM+CAL), XPRINT(EQLAM), XDATABASE:[DIFF(W,AMP)=0,DIFF(W,LAM)=0], XSUB:Y=AMP*CFUN+W, XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), XNULLANS:READ(" ARE ZERO, Y/N") X)$ X X XG_POLY(I,J):=BLOCK( X/* THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) =0. IT REQUIRES KNOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE */ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, X Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*THIS DIFFERENTIAL OF W WILL BE ANNIHILATED BY THE PROJECTION ONTO THE XCRITICAL EIGENFUNCTION, HENCE WE SET IT TO ZERO ALREADY HERE */ X TEMP3:SUBST('DIFF(W,AMP,I,LAM,J)=0,TEMP2), X/* HERE WE INSERT ALL KNOWLEDGE GAINED THROUGH W_POLY */ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X/*PROJECTION ONTO CFUN, THE CRITICAL EIGENFUNCTION */ X TEMP5:TRIGREDUCE(CFUN*TEMP4), X BIFEQ[I,J]:RATSIMP(INTEGRATE(TEMP5,X,0,LEN)) X )$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). THE RESULT IS FED INTO DATABASE.*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*WE SUBSTITUTE ZW FOR THE UNKNOWN TAYLOR COEFFICIENT AND SOLVE XAN O.D.E. FOR ZW IN TEMP7 */ X TEMP3:SUBST(DIFF(W,AMP,I,LAM,J)=ZW,TEMP2), X/*NOW WE INSERT ALL KNOWLEDGE GAINED THROUGH PREVIOUS ITERATIONS*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X TEMP5:TRIGREDUCE(TEMP4), X/*THIS ENSURES THAT THE R.H.S. FO THE D.E. TEMP6 IS ORTHOGONAL TO XTHE SOLUTION OF THE HOMOGENEOUS EQUATION */ X TEMP6:PROJECT(TEMP5), X TEMP7:ODE2(TEMP6,ZW,X), X/*SATISFY BOUNDARY CONDITIONS*/ X IF BCL*BCR=1 THEN TEMP8:BC2(TEMP7,X=0,ZW=0,X=LEN,ZW=0) X ELSE X TEMP8:NEUMANNBC(TEMP7), X/*MAKE SURE THAT = 0*/ X TEMP9:PROJECT(TEMP8), X ADDBASE:['DIFF(W,AMP,I,LAM,J)=RHS(TEMP9)], X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X XPROJECT(EXP):=( X PRO1:EV(EXP,ZW=0), X PRO2:INTEGRATE(PRO1*CFUN,X,0,LEN)/NORM, X EXPAND(EXP-PRO2*CFUN))$ XNEUMANNBC(EXP):=( X REXP:RHS(EXP), X NBC1:DIFF(REXP,X), X IF BCL=1 AND BCR=2 THEN NBC2:SOLVE([EV(REXP,X=0),EV(NBC1,X=LEN)], X [%K1,%K2]), X IF BCL=2 AND BCR=1 THEN NBC2:SOLVE([EV(REXP,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X IF BCL=2 AND BCR=2 THEN NBC2:SOLVE([EV(NBC1,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X EV(EXP,NBC2))$ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X END_OF_FILE if test 5930 -ne `wc -c <'program10'`; then echo shar: \"'program10'\" unpacked with wrong size! fi # end of 'program10' fi if test -f 'program11' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program11'\" else echo shar: Extracting \"'program11'\" \(4141 characters\) sed "s/^X//" >'program11' <<'END_OF_FILE' X/* PROGRAM NUMBER 11: REDUCTION2(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN SYSTEMS OF ORDINARY DIFFERENTIAL X EQUATIONS. SEE PAGE 176 AND PAGE 211 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA".*/ X X X X X XREDUCTION2():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X X XSETUP():=(/*THE FUNCTION SETUP NEEDS THE DIMENSION N OF THE SYSTEM OF X O.D.E., THE N-DIMENSIONAL LIST OF VARIABLES Y, THE SYSTEM X OF O.D.E'S CALLED EQS, THE CRITICAL EIGENVECTOR CFUN AND X THE ADJOINT CRITICAL EIGENVECTOR AFUN.*/ X PRINTFLAG:TRUE, X LS_LIST:[], X N:READ("NUMBER OF EQUATIONS"), X Y:MAKELIST(READ("ENTER VARIABLE NUMBER",I),I,1,N), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE",ALPHA), X PRINT("WE DEFINE LAM =",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENVECTOR AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENVECTOR"), X KILL(W,ZWLIST), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W),[AMP,LAM]), X PRINT("ENTER THE DIFFERENTIAL EQUATION"), X EQS:MAKELIST(READ("DIFF(",Y[I],",T)="),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT(EQLAM), X WNULL:VFUN(W,0), X ZWNULL:VFUN(ZWLIST,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X DATABASE:APPEND(DATABASE,MAP(LAMBDA([U],DIFF(U,LAM)=0),W)), X NORM:AFUN.CFUN, X TEMP1:EV(EQLAM,SUB,DIFF), X PRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), X NULLANS:READ(" ARE ZERO, Y/N") X )$ X X XG_POLY(I,J):=BLOCK(/*PROVIDED ALL NECESSARY KNOWLEDGE ABOUT THE TAYLOR X COEFFICIENTS OF W(AMP,LAM) IS STORED IN THE LIST X DATABASE, G_POLY CALCULATES ANY SPECIFIC TAYLOR X COEFFICIENT OF THE BIFURCATION EQUATION X G(AMP,LAM)= 0 VIA THE PROJECTION ONTO CFUN.*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X BIFEQ[I,J]:RATSIMP(TEMP4.AFUN))$ X X X XW_POLY(I,J):=BLOCK(/*THE FUNCTION W_POLY ALLOWS TO CALCULATE ITERATIVELY, X I.E. STARTING WITH THE LOWEST ORDER, ALL TAYLOR X COEFFICIENTS OF W(AMP,LAM) BY PROJECTING ONTO THE X COMPLEMENT OF CFUN*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL), X TEMP5:EV(TEMP4,ZWNULL).AFUN, X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:SOLVE(TEMP6,ZWLIST), X TEMP8:EV(ZWLIST,TEMP7), X TEMP9:RATSIMP(TEMP8.AFUN), X IF TEMP9 = 0 THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),TEMP8) X ELSE X (TEMP10:RATSIMP(TEMP8 - TEMP9/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J), X U,1,N),TEMP10)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X X X X X X END_OF_FILE if test 4141 -ne `wc -c <'program11'`; then echo shar: \"'program11'\" unpacked with wrong size! fi # end of 'program11' fi if test -f 'program12' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program12'\" else echo shar: Extracting \"'program12'\" \(8632 characters\) sed "s/^X//" >'program12' <<'END_OF_FILE' X X/* PROGRAM NUMBER 12: A LIAPUNOV-SCHMIDT REDUCTION FOR STEADY STATE X BIFURCATIONS IN SYSTEMS OF PARTIAL DIFFERENTIAL EQUATIONS X DEPENDING ON ONE INDEPENDENT SPACE VARIABLE. SEE PAGE 189 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". X*/ X X X X X/*THIS FILE CONTAINS REDUCTION3(), A FUNCTION TO PERFORM A LIAPUNOV- XSCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS FROM N COUPLED XPARTIAL DIFFERENTIAL EQUATIONS ON ONE SPATIAL DIMENSION. XTHE FOLLOWING ADDITIONAL FUNCTIONS ARE SUPPLIED: X SETUP() ALLOWS THE PROBLEM TO BE ENTERED. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM). X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X SOLVE_ODE(EXP) SOLVES CERTAIN ORDINARY DIFFERENTIAL EQUATIONS VIA A FOURIER X MODE ANSATZ. X FEEDW(EXP) ENSURES THAT = 0 . X FIND_TRIG(EXP) IDENTIFIES FOURIER MODES. X SETIFY(LIST) TRANSFORMS A LIST INTO A SET. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE, ...] X DIFFNULL(I,J) SETS THE DERIVATIVE DIFF(W,AMP,I,LAM,J) TO ZERO.*/ X XREDUCTION3():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X XSETUP():=( /* THIS FUNCTION PERFORMS THE INPUT FOR THE LIAPUNOV-SCHMIDT X REDUCTION*/ XASSUME_POS:TRUE, XLS_LIST:[], XN:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), XY:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), XXVAR:READ("ENTER THE SPATIAL COORDINATE"), XALPHA:READ("ENTER THE BIFURCATION PARAMETER"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), XPRINT("WE DEFINE LAM = ",ALPHA - CAL), XCFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), XAFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS A LIST"), XKILL(W), XW:MAKELIST(CONCAT(W,I),I,1,N), XZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), XNULLIST:MAKELIST(0,I,1,N), XDEPENDS(APPEND(ZWLIST,W,Y),CONS(XVAR,[AMP,LAM])), XEQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), XEQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), XPRINT(EQLAM), XLEN:READ("WHAT IS THE LENGTH OF THE SPACE INTERVAL"), XWNULL:VFUN(W,0), XSUB:MAPLIST("=",Y,AMP*CFUN+W), XDATABASE:APPEND(DIFNULL(1,0),DIFNULL(0,1)), XZWNULL:VFUN(ZWLIST,0), XNORM:INTEGRATE(AFUN.CFUN,XVAR,0,LEN), XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICIENTS ARE 0"), XNULLANS:READ("Y,N") X)$ X X XG_POLY(I,J):=BLOCK( X/*THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) = 0. IT REQUIRES KOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS G_POLY(",I,",",J,")IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*SET THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZERO*/ X TEMP3:SUBST(DIFNULL(I,J),TEMP2), X /*ENTER ALL INFORMATION IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE)), X /*PROJECT ONTO AFUN*/ X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ[I,J]:INTEGRATE(TRIGREDUCE(TEMP5),XVAR,0,LEN))$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). IT RETURNS A DIFFERENTIAL EQUATION XFOR THE PARTICULAR COEFFICIENT OF INTEREST (CALLED ZW1,ZW2...). THIS D.E. XIS SOLVED VIA SOLVE_ODE AND THE RESULT IS FED INTO DATABASE FROM FEEDW*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W(AMP,",I,",LAM,",J,") IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN ( X ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*RENAME THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZW */ X TEMP3:SUBST(WMAX_DIFF,TEMP2), X /*ENTER ALL INFORMATION IN STORED IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE), X /*THIS IS THE PROJECTION Q ONTO THE RANGE OF THE X LINEAR DIFFERENTIAL OPERATOR IN THE PROBLEM*/ X TEMP5:INTEGRATE(EV(TEMP4,ZWNULL).AFUN,XVAR,0,LEN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:TRIGREDUCE(TEMP6), X /*THE SET OF O.D.E.'S TO SOLVE */ X W_DE:EXPAND(TEMP7), X TEMP8:EV(W_DE,VFUN(ZWLIST,0)), X /*IF THE PARTICULAR SOLUTION OF W_DE IS ZERO THEN W=0 */ X IF NULLIST=TEMP8 THEN ( ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE)), X TEMP9:SOLVE_ODE(TEMP8), X FEEDW(TEMP9) X )$ X X XSOLVE_ODE(EXP):=( X/*THIS FUNCTION SOLVE THE D.E. W_DE BY A FOURIER MODE ANSATZ*/ X TRIGFUN:[], X CONST:FALSE, X FOR I THRU N DO ( X /*DETERMINE THE FOURIER MODES*/ X TRIG1:EXP[I], X IF TRIG1 # 0 THEN ( X TRIG2:APPLY1(TRIG1,SINNULL,COSNULL), X IF TRIG2 # 0 THEN ( X CONST:TRUE, X TRIG1:TRIG1-TRIG2), X TRIGFUN:APPEND(FIND_TRIG(TRIG1),TRIGFUN))), X SOL1:DELETE(DUMMY,SETIFY(TRIGFUN)), X /*MAKE AN ANSATZ*/ X ANSATZ:MAKELIST(SUM(AM[I,J]*SOL1[I],I,1,LENGTH(SOL1)),J,1,N), X SOL2:EV(W_DE,MAP("=",ZWLIST,ANSATZ),DIFF), X SOL3:MAKELIST(RATCOEF(SOL2,I),I,SOL1), X EQLIST:[], X FOR I THRU LENGTH(SOL3) DO EQLIST:APPEND(EQLIST,SOL3[I]), X VARLIST:[], X FOR I THRU N DO X FOR J THRU LENGTH(SOL1) DO VARLIST:CONS(AM[J,I],VARLIST), X /*FIND THE AMPLITUDES OF THE FOURIER MODES*/ X SOL4:SOLVE(EQLIST,VARLIST), X /*SOLVE FOR THE CONSTANT FOURIER MODE IF NECESSARY*/ X CANSATZ:0, X IF CONST = TRUE THEN (CANSATZ:MAKELIST(CONCAT(CON,I),I,1,N), X CSOL1:EV(W_DE,MAP("=",ZWLIST,CANSATZ),DIFF), X CSOL2:APPLY1(CSOL1,SINNULL,COSNULL), X CSOL3:SOLVE(CSOL2,CANSATZ)), X EV(ANSATZ+CANSATZ,SOL4,CSOL3))$ X X XFEEDW(EXP):=( X/* THIS FUNCTION ALLOWS THE RESULT OF THE ODE-SOLVER TO BE ENTERED INTO X THE LIST DATABASE. IT CHECKS FOR ORTHOGONALITY TO THE CRITICAL ADJOINT X EIGENFUNCTION AND REMOVES SOLUTIONS OF THE HOMOGENEOUS EQUATION (I.E. X NONORTHOGONAL TERMS)*/ X F2:INTEGRATE(EXP.AFUN,XVAR,0,LEN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X EXP) X ELSE (ORTHO_RESULT:RATSIMP(EXP- F2/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X/*COLLECT ALL INFORMATION STORED IN BIFEQ AND ASSEMBLE THE BIFURCATION XEQUATION*/ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X XSETIFY(LIST):=(/*TRANSFORMS A LIST INTO A SET, I.E. REMOVES DUPLICATES*/ X SET:[LIST[1]], X FOR I :2 THRU LENGTH(LIST) DO X IF NOT MEMBER(LIST[I],SET) THEN X SET:CONS(LIST[I],SET), X SET)$ X XFIND_TRIG(EXP):=(/*FINDS THE FOURIER MODES*/ X F_A1:ARGS(EXP+DUMMY), X F_A2:APPLY1(F_A1,SINFIND,COSFIND) X )$ X X/*AUXILIARY FUNCTIONS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE)$ XDEFRULE(COSFIND,DUMMY1*COS(DUMMY2),COS(DUMMY2))$ XDEFRULE(SINFIND,DUMMY1*SIN(DUMMY2),SIN(DUMMY2))$ XDEFRULE(COSNULL,COS(DUMMY1),0)$ XDEFRULE(SINNULL,SIN(DUMMY1),0)$ X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ XDIFNULL(I,J):=MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J)=0),W) $ X END_OF_FILE if test 8632 -ne `wc -c <'program12'`; then echo shar: \"'program12'\" unpacked with wrong size! fi # end of 'program12' fi if test -f 'program13' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program13'\" else echo shar: Extracting \"'program13'\" \(4145 characters\) sed "s/^X//" >'program13' <<'END_OF_FILE' X X/*PROGRAM NUMBER 13: CONTAINS THE BUILDING BLOCKS TO PERFORM STEADY X STATE BIFURCATIONS FOR MORE PARTIAL DIFFERENTIAL EQUATIONS. SEE PAGES X 198-202 IN "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER X ALGEBRA" */ X X X X X/*THE FOLLOWING FUNCTIONS CAN BE USED TO PERFORM A LIAPUNOV-SCHMIDT REDUCTION X FOR THE 2D BENARD PROBLEM */ X X/*SETUP() ALLOWS THE PROBLEM TO BE ENTERED, X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM), X W_POLY(I,J) DETERMINES A DIFFERENTIAL EQUATION, WHOSE SOLUTION IS THE X COEFFICIENT OF AMP^I LAM^J IN W(AMP,LAM), X FEEDW() ALLOWS THIS SOLUTION TO BE ENTERED INTO THE LIST DATABASE, X INT(EXP) FINDS MULTIPLE DEFINITE INTEGRALS EFFECTIVELY, X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE ...] X*/ X XSETUP():=( X N:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), X Y:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), X SPACE:READ("ENTER NUMBER OF SPATIAL COORDINATES"), X XVAR:READ("ENTER THE SPATIAL COORDINATES AS A LIST"), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), X PRINT("WE DEFINE LAM = ",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS ALIST"), X KILL(W), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W,Y),APPEND([AMP,LAM],XVAR)), X EQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT("ENTER THE LOWER LEFT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X LBOUND:READ("[",X[1],"=...,]"), X PRINT("ENTER THE UPPER RIGHT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X UBOUND:READ("[",X[1],"=...,]"), X WNULL:VFUN(W,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X ZWNULL:VFUN(ZWLIST,0), X NORM:INT(AFUN.CFUN), X TEMP1:EV(EQLAM,SUB,DIFF), X EQLAM X )$ X X X XG_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ:INT(TEMP5))$ X X X X XW_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:INT(EV(TEMP4,ZWNULL).AFUN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X FOR I THRU SPACE DO TEMP7:TRIGREDUCE(TEMP6,XVAR[I]), X W_DE:RATSIMP(TEMP7 ), X PRINT("NOW SOLVE THE EQUATIONS",W_DE,"=0! THEY ARE GIVEN IN W_DE"), X PRINT("CALL FEEDW() TO PROCEED"))$ X X X XFEEDW():=( X ADDBASE:[], X RESULT:MAKELIST((PRINT("ENTER",ZWLIST[K]),READ()),K,1,N), X F2:INT(RESULT.AFUN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X RESULT) X ELSE (ORTHO_RESULT:RATSIMP(RESULT- F2/NORM*CFUN), X/*PROJECTION Q*/ X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X ADDBASE)$ X X XINT(EXP):=(%INT:EXP,FOR I THRU SPACE DO %INT:INTEGRATE(TRIGREDUCE(%INT, X XVAR[I]),XVAR[I]), X RATSIMP(EV(%INT,UBOUND)-EV(%INT,LBOUND)))$ X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ END_OF_FILE if test 4145 -ne `wc -c <'program13'`; then echo shar: \"'program13'\" unpacked with wrong size! fi # end of 'program13' fi if test -f 'program14' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program14'\" else echo shar: Extracting \"'program14'\" \(39932 characters\) sed "s/^X//" >'program14' <<'END_OF_FILE' X/* PROGRAM NUMBER 6: TWOVAR(), A FUNCTION TO PERFORM A TWO VARIABLE X EXPANSION TO ORDER EPSILON FOR N WEAKLY COUPLED PERTURBED NONLINEAR X OSCILLATORS. SEE PAGE 94 IN "PERTURBATION METHODS, BIFURCATION X THEORY AND COMPUTER ALGEBRA". */ X X X XTWOVAR():=BLOCK( X XCHOICE:READ("DO YOU WANT TO ENTER NEW DATA (Y/N)"), X XIF CHOICE = N THEN GO(JUMP1), X X/* INPUT D.E.'S */ XN:READ("NUMBER OF D.E.'S"), X XPRINT("THE",N,"D.E.'S WILL BE IN THE FORM:"), XPRINT("X[I]'' + W[I]^2 X[I] = E F[I](X[1],...,X[",N,"],T)"), X XFOR I:1 THRU N DO X X[I]:READ("ENTER SYMBOL FOR X[",I,"]"), X XFOR I:1 THRU N DO X DEPENDS(X[I],T), X XFOR I:1 THRU N DO X W[I]:READ("ENTER W[",I,"]"), X XFOR I:1 THRU N DO X F[I]:READ("ENTER F[",I,"]"), X XJUMP1, X X/* UPDATE EQS FOR SUBSTITUTION OF RESONANT VALUES ON 2ND TIME THRU */ XFOR I:1 THRU N DO( X W[I]:EV(W[I]), X F[I]:EV(F[I])), X X/* ECHO BACK THE D.E.'S */ XPRINT("THE D.E.'S ARE ENTERED AS:"), XFOR I:1 THRU N DO X PRINT(X[I],"'' +",W[I]^2*X[I],"=",E*F[I]), X XPRINT("THE METHOD ASSUMES A SOLUTION IN THE FORM:"), XPRINT("X[I] = X0[I] + E X1[I]"), XPRINT("WHERE X0[I] = A[I](ETA) COS W[I] XI + B[I](ETA) SIN W[I] XI"), XPRINT("WHERE XI = T AND ETA = E T"), X X/* DON'T USE A OR B AS PARAMETERS IN THE GIVEN D.E.'S */ XDEPENDS([A,B],ETA), X XFOR I:1 THRU N DO X X0[I]:A[I]*COS(W[I]*XI)+B[I]*SIN(W[I]*XI), X XFOR I:1 THRU N DO X G[I]:EV(F[I],T=XI), X XFOR I:1 THRU N DO( X FOR J:1 THRU N DO X G[I]:EV(G[I],X[J]::X0[J])), X XFOR I:1 THRU N DO( X G[I]:G[I]-2*DIFF(X0[I],XI,1,ETA,1), X G[I]:EV(G[I],DIFF), X G[I]:EXPAND(TRIGREDUCE(EXPAND(G[I])))), X X/* COLLECT SECULAR TERMS */ XFOR I:1 THRU N DO( X S[I]:COEFF(G[I],SIN(W[I]*XI)), X C[I]:COEFF(G[I],COS(W[I]*XI))), X X/* DISPLAY SECULAR TERMS */ XPRINT("REMOVAL OF SECULAR TERMS IN THE X1[I] EQS. GIVES:"), X XFOR I:1 THRU N DO( X PRINT(S[I],"= 0"), X PRINT(C[I],"= 0")), X XABEQS:APPEND(MAKELIST(S[I],I,1,N),MAKELIST(C[I],I,1,N)), X XCHOICE2:READ("DO YOU WANT TO TRANSFORM TO POLAR COORDINATES (Y/N)"), X XIF CHOICE2=N THEN GO(JUMP2), X X/* TRANSFORM TO POLAR COORDINATES */ XDEPENDS([R,THETA],ETA), XTRANS:MAKELIST([A[I]=R[I]*COS(THETA[I]),B[I]=R[I]*SIN(THETA[I])],I,1,N), XINTEQS:EV(ABEQS,TRANS,DIFF), XRTHEQS:SOLVE(INTEQS,APPEND(MAKELIST(DIFF(R[I],ETA),I,1,N), X MAKELIST(DIFF(THETA[I],ETA),I,1,N))), XRTHEQS:TRIGSIMP(RTHEQS), XRTHEQS:EXPAND(TRIGREDUCE(RTHEQS)), XPRINT(RTHEQS), X XJUMP2, X XCHOICE3:READ("DO YOU WANT TO SEARCH FOR RESONANT PARAMETER VALUES (Y/N)"), X XIF CHOICE3=N THEN GO(END), X X/* OBTAIN FREQUENCIES WHICH APPEAR ON RHS'S OF D.E.'S */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE FREQS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR I:1 THRU N DO( X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(G[I]), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X TEMP3:APPLY1(TEMP2,COSARG,SINARG), X TEMP4:EV(TEMP3,XI=1), X/* REMOVE DUPLICATE ARGUMENTS */ X FREQ[I]:SETIFY(TEMP4)), X X/* DISPLAY FREQUENCIES */ XFOR I:1 THRU N DO( X PRINT(X[I],"EQ'S RESONANT FREQ =",W[I]), X PRINT("FREQS ON RHS =",FREQ[I])), X XJUMP3, X XPAR:READ("WHICH PARAMETER TO SEARCH FOR ?"), X X/* COMPUTE RESONANT VALUES OF DESIRED PARAMETER */ XRESVALS:[], XFOR I:1 THRU N DO( X FOR J:1 THRU LENGTH(FREQ[I]) DO( X RES:SOLVE(PART(FREQ[I],J)=W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES), X RES:SOLVE(PART(FREQ[I],J)=-W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES))), X X/* ELIMINATE DUPLICATE PARAMETER VALUES */ XRESVALUES:SETIFY(RESVALS), X X/* DISPLAY RESONANT PARAMETER VALUES */ XPRINT(RESVALUES), X XCHOICE4:READ("DO YOU WANT TO SEARCH FOR ANOTHER PARAMETER (Y/N) ?"), X XIF CHOICE4=Y THEN GO(JUMP3), X XEND," ")$ X X/* AUXILIARY FUNCTIONS */ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X X X X X X X X X/************************************************************************/ X/************************************************************************/ X X X X/*PROGRAM NUMBER 7: AVERAGE(), ALLOWS TO PERFORM FIRST AND SECOND ORDER X AVERAGING. SEE PAGE 114 IN "PERTURBATION METHODS, BIFURCATION THEORY X AND COMPUTER ALGEBRA". */ X X X X/* PROGRAM TO PERFORM 1ST OR 2ND ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERAGE():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("DO YOU WANT FIRST OR SECOND ORDER AVERAGING? (1/2)"), X COEFVFEPS1:COEFF(VF2,EPS), X COEFVFEPS2:COEFF(VF2,EPS,2), X G1BAR:DEMOIVRE(APPLY1(COEFVFEPS1,KILLEXP)), X IF M=1 THEN RESULT:EPS*G1BAR X ELSE( X G1TILDE:COEFVFEPS1-G1BAR, X W1:INTEGRATE(G1TILDE,T), X REMARRAY(JACOB), X JACOB[I,J] := DIFF(G1TILDE[I],X[J]), X JACG1TILDE:GENMATRIX(JACOB,N,N), X G2BAR:DEMOIVRE(APPLY1(EXPAND(JACG1TILDE.W1)+COEFVFEPS2,KILLEXP)), X RESULT:EPS*G1BAR+EPS^2*G2BAR), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS TO KILL EXPONENTIAL TERMS */ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X X X X X X X X X/*************************************************************************/ X/*************************************************************************/ X X X X/* PROGRAM NUMBER 8: AVERM(), ALLOWS MTH ORDER AVERAGING. SEE PAGE 128 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X/* PROGRAM TO PERFORM MTH ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERM():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("ENTER ORDER OF AVERAGING"), X IF M=1 X/* FIRST ORDER AVERAGING */ X THEN (TEMP0:DEMOIVRE(APPLY1(F,KILLEXP)), X RESULT:TAYLOR(TEMP0,EPS,0,1)) X/* AVERAGING OF ORDER M > 1 */ X ELSE X Y:MAKELIST(CONCAT('Y,I),I,1,N), X WLIST:MAKELIST(CONCAT('W,I),I,1,N), X DEPENDS(WLIST,CONS(T,Y)), X TRAFO:Y, X XSUB:MAPLIST("=",Y,X), X/* WNULL SETS WLIST TO ZERO */ X WNULL:VFUN(WLIST,0), X JACOB[I,J] := DIFF(WLIST[I],Y[J]), X JAC:GENMATRIX(JACOB,N,N), X/* MAIN LOOP */ X FOR I :1 THRU M-1 DO ( X TEMP1:MAPLIST("=",X,Y+EPS^I*WLIST), X/* HERE X IS THE LIST [X1,X2,...,XN], Y IS THE LIST [Y1,Y2,...,YN], X WLIST IS THE TRANSFORMATION VECTOR [W1,W2,...,WN] */ X TEMP2:TAYLOR(SUM((-EPS)^(I*J)*JAC^^J,J,0,M-1). X (EV(VF2,TEMP1) - DIFF(WLIST,T)*EPS^I),EPS,0,M), X/* JAC IS THE JACOBIAN D WLIST/DY OF THE TRANSFORMATION WLIST */ X TEMP3:MATTOLIST(TEMP2,N), X TEMP4:MAP(EXPAND,TAYLOR(EV(TEMP3,WNULL,DIFF),EPS,0,I)), X/* THE ITH ORDER MEAN */ X MEAN:APPLY1(TEMP4,KILLEXP), X TEMP6:EXPAND(TEMP4-MEAN), X TEMP7:EV(INTEGRATE(TEMP6,T),EPS=1), X/* THE ITH ORDER TRANSFORMATION */ X TEMP8:MAPLIST("=",WLIST,TEMP7), X TEMP9:RATSIMP(TEMP8), X/* THE TRANSFORMED DE */ X VF2:EXPAND(EV(TEMP3,TEMP9,DIFF,XSUB,INFEVAL)), X/* THE SUM OF ALL TRANSFORMATIONS */ X TRAFO:EXPAND(TRAFO+EV(EPS^I*WLIST,TEMP9))), X/* END OF MAIN LOOP */ X PRINT("THE TRANSFORMATION: ",X,"="), X PRINT(RATSIMP(DEMOIVRE(TRAFO))), X/* THE FINAL AVERAGING */ X RESULT:APPLY1(VF2,KILLEXP), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS */ X VFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X MATTOLIST(MAT,DIM):=IF DIM>1 THEN MAKELIST(MAT[I,1],I,1,DIM) ELSE [MAT]$ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X X X X X X X X X/***********************************************************************/ X/***********************************************************************/ X X X X/* PROGRAM NUMBER 9: LIE(), LIE TRANSFORMATIONS FOR HAMILTONIAN SYSTEMS. X SEE PAGE 144 IN "PERTURBATION METHODS, BIFURCATION THEORY AND X COMPUTER ALGEBRA". */ X X XLIE():=BLOCK( X X/* INPUT PROBLEM ? */ XCHOICE1:READ("DO YOU WANT TO INPUT A NEW PROBLEM (Y/N) ?"), X XIF CHOICE1=N THEN GO(JUMP1), X X/* INPUT PROBLEM */ XN:READ("ENTER NUMBER OF DEGREES OF FREEDOM"), XFOR II:1 THRU N DO ( X Q[II]:READ("ENTER SYMBOL FOR Q[",II,"]"), X P[II]:READ("ENTER SYMBOL FOR P[",II,"]")), XKILL(W), XPRINT("THE HAMILTONIAN DEPENDS ON THE Q'S, P'S, T AND E (SMALL PARAMETER)"), XPRINT("THE E=0 PROBLEM MUST BE OF THE FORM:"), XPRINT("H =",SUM((P[II]^2+W[II]^2*Q[II]^2)/2,II,1,N)), XHORIGINAL:READ("ENTER THE HAMILTONIAN"), XPRINT("H =",HORIGINAL), X X/* TRANSFORM TO ACTION-ANGLE VARIABLES */ X/* FIND THE W[II]'S */ XH0:EV(HORIGINAL,E=0), XFOR II:1 THRU N DO X W[II]:SQRT(DIFF(H0,Q[II],2)), XPRINT("THE ACTION-ANGLE VARIABLES ARE J'S FOR ACTION, PHI'S FOR ANGLE"), XFOR II:1 THRU N DO X TR[II]:[Q[II]=SQRT(2*J[II]/W[II])*SIN(PHI[II]), X P[II]=SQRT(2*J[II]*W[II])*COS(PHI[II])], XTRAN:MAKELIST(TR[II],II,1,N), XH:EV(HORIGINAL,TRAN,ASSUME_POS:TRUE,INFEVAL), XH:TRIGSIMP(H), XH:EXPAND(TRIGREDUCE(EXPAND(H))), XPRINT("H =",H), X XJUMP1, X X/* INPUT TRUNCATION ORDER */ XNTRUNC:READ("ENTER HIGHEST ORDER TERMS IN E TO BE KEPT"), XFOR II:0 THRU NTRUNC DO X H[II]:RATCOEF(H,E,II), X X/* LIE TRANSFORMS */ X/* NEAR IDENTITY TRANSFORMATION FROM (J,PHI)'S TO (I,PSI)'S */ X X/* UPDATE VARIABLES */ XFOR II:1 THRU N DO( X P[II]:I[II], X Q[II]:PSI[II]), X X/* REPLACE J AND PHI BY I AND PSI IN H'S */ XFOR II:0 THRU NTRUNC DO X H[II]:EV(H[II],MAKELIST(J[III]=I[III],III,1,N), X MAKELIST(PHI[III]=PSI[III],III,1,N)), X XK[0]:H[0], X X/* DECLARE WGEN[I] TO BE A FN OF T, Q'S AND P'S */ XKILL(WGEN), XDEPENDS(WGEN,[T]), XFOR II:1 THRU N DO X DEPENDS(WGEN,[Q[II],P[II]]), X X/* E=0 PROBLEM IS OF FORM SUM(W[II]*I[II]) */ X/* CHOOSE WGEN[II] TO KILL AS MUCH AS POSSIBLE IN EQ(II) */ X/* EQUATE K[II] TO UNREMOVABLE TERMS */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE ARGS OF TRIG TERMS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR II:1 THRU NTRUNC DO( X X EQN[II]:EXPAND(TRIGREDUCE(EXPAND(EQ(II)))), X TEMP:EXPAND(EV(EQN[II],WGEN[II]=0)), X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(TEMP), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X ARG1:APPLY1(TEMP2,COSARG,SINARG), X/* REMOVE DUPLICATE ARGUMENTS */ X ARG2:SETIFY(ARG1), X/* REMOVE RESONANT ARGUMENTS */ X ARG3:SUBLIST(ARG2,NOTRESP), X X/* CHOOSE WGEN TO ELIMINATE NONRESONANT TERMS */ X LENG:LENGTH(ARG3), X WGENTEMP1:0, X FOR JJ:1 THRU LENG DO( X WGENTEMP2:AAA*COS(PART(ARG3,JJ))+BBB*SIN(PART(ARG3,JJ)), X TEMP4:EV(EQN[II],WGEN[II]=WGENTEMP2,DIFF), X TEMP5:SOLVE([RATCOEF(TEMP4,COS(PART(ARG3,JJ))), X RATCOEF(TEMP4,SIN(PART(ARG3,JJ)))],[AAA,BBB]), X WGENTEMP1:WGENTEMP1+EV(WGENTEMP2,TEMP5)), X WGEN[II]:WGENTEMP1, X PRINT("WGEN[",II,"] ="), X PRINT(WGEN[II]), X K[II]:EXPAND(EV(EQN[II],DIFF)), X K[II]:EXPAND(RATSIMP(K[II])), X PRINT("THE TRANSFORMED HAMILTONIAN K[",II,"] ="), X PRINT(K[II])), X XKAMILTONIAN:SUM(K[II]*E^II,II,0,NTRUNC), XPRINT ("THE TRANSFORMED HAMILTONIAN IS "), XPRINT ("K =",KAMILTONIAN), X XCHOICE2:READ("DO YOU WANT TO SEE THE NEAR IDENTITY TRANSFORMATION (Y/N) ?"), XIF CHOICE2=N THEN GO(END), X X/* THE NEAR IDENTITY TRANSFORMATION */ XFOR II:1 THRU N DO( X JTRANS[II]:SUM(S(III,P[II])*E^III,III,0,NTRUNC), X PHITRANS[II]:SUM(S(III,Q[II])*E^III,III,0,NTRUNC)), XFOR II:1 THRU N DO( X PRINT(J[II],"=",JTRANS[II]), X PRINT(PHI[II],"=",PHITRANS[II])), X XEND, XKAMILTONIAN)$ X X XPOISSON(F,G):= X SUM(DIFF(F,Q[II])*DIFF(G,P[II])-DIFF(F,P[II])*DIFF(G,Q[II]),II,1,N)$ X XL(II,F):=POISSON(WGEN[II],F)$ X XS(II,F):=(IF II=0 THEN F ELSE SUM(L(II-M,S(M,F)),M,0,II-1)/II)$ X XEQ(II):=(H[II]+(DIFF(WGEN[II],T)+POISSON(WGEN[II],H[0]))/II X +SUM(L(II-M,K[M])+M*S(II-M,H[M]),M,1,II-1)/II)$ X XLZAP(ANY):=DIFF(ANY,T)+POISSON(ANY,H[0])$ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XNOTRESP(EXP):=NOT IS(LZAP(EXP) = 0)$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X X X X X X X X X X/**************************************************************************/ X/**************************************************************************/ X X X/*PROGRAM NUMBER 10: REDUCTION1(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN ONE DIFFERENTIAL EQUATION DEPENDING X ON ONE INDEPENDENT VARIABLE. SEE PAGE 170 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X X X/*THIS FILE CONTAINS REDUCTION1(), A FUNCTION TO PERFORM A LIAPUNOV- X SCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS OF NONLINEAR D.E.'S X OF THE FORM Y'' + F(Y,Y',ALPHA) = 0. Y = Y(X) IS DEFINED ON A REAL X INTERVAL WITH DIRICHLET OR NEUMANN BOUNDARY CONDITIONS AND F DEPENDS X ONLY LINEARLY ON ALPHA. X IT ALSO CONTAINS THESE ADDITIONAL FUNCTIONS: X SETUP() ALLOWS TO ENTER THE PROBLEM. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION. X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X PROJECT(EXP) ENSURES THAT =0. X NEUMANNBC(EXP) ACCOUNTS FOR NEUMANN BOUNDARY CONDITIONS. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. */ X X XREDUCTION1():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X XSETUP():=( X/*THE FUNCTION SETUP ASKS FOR THE VARIABLES OF THE D.E., THE BIFURCATION XPOINT, THE CRITICAL EIGENFUNCTION, THE X-INTERVAL, THE BOUNDARY XCONDITIONS AND THE DIFFERENTIAL EQUATION. */ XASSUME_POS:TRUE, XLS_LIST:[], XY:READ("ENTER DEPENDENT VARIABLE"), XPRINT("USE X AS THE INDEPENDENT VARIABLE AND ALPHA AS A PARAMETER TO VARY"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE ALPHA"), XPRINT("WE DEFINE LAM = ALPHA - ",CAL), XCFUN: READ("ENTER THE CRITICAL EIGENFUNCTION"), XDEPENDS([ZW,Y,W],[AMP,LAM,X]), XLEN:READ("WHAT IS THE LENGTH OF THE X-INTERVAL"), XNORM:INTEGRATE(CFUN^2,X,0,LEN), XPRINT("SPECIFY THE BOUNDARY CONDITIONS"), XPRINT("YOUR CHOICE FOR THE B.C. ON Y AT X=0 AND X=",LEN), XPRINT("ENTER 1 FOR Y=0, 2 FOR Y'=0"), XBCL:READ("B.C. AT 0?"), XBCR:READ("B.C. AT",LEN,"?"), XEQ:DIFF(Y,X,2) + X READ("THE D.E. IS OF THE FORM Y'' + F(Y,Y',ALPHA) = 0,ENTER F"), XEQLAM:EV(EQ,ALPHA=LAM+CAL), XPRINT(EQLAM), XDATABASE:[DIFF(W,AMP)=0,DIFF(W,LAM)=0], XSUB:Y=AMP*CFUN+W, XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), XNULLANS:READ(" ARE ZERO, Y/N") X)$ X X XG_POLY(I,J):=BLOCK( X/* THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) =0. IT REQUIRES KNOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE */ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, X Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*THIS DIFFERENTIAL OF W WILL BE ANNIHILATED BY THE PROJECTION ONTO THE XCRITICAL EIGENFUNCTION, HENCE WE SET IT TO ZERO ALREADY HERE */ X TEMP3:SUBST('DIFF(W,AMP,I,LAM,J)=0,TEMP2), X/* HERE WE INSERT ALL KNOWLEDGE GAINED THROUGH W_POLY */ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X/*PROJECTION ONTO CFUN, THE CRITICAL EIGENFUNCTION */ X TEMP5:TRIGREDUCE(CFUN*TEMP4), X BIFEQ[I,J]:RATSIMP(INTEGRATE(TEMP5,X,0,LEN)) X )$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). THE RESULT IS FED INTO DATABASE.*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*WE SUBSTITUTE ZW FOR THE UNKNOWN TAYLOR COEFFICIENT AND SOLVE XAN O.D.E. FOR ZW IN TEMP7 */ X TEMP3:SUBST(DIFF(W,AMP,I,LAM,J)=ZW,TEMP2), X/*NOW WE INSERT ALL KNOWLEDGE GAINED THROUGH PREVIOUS ITERATIONS*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X TEMP5:TRIGREDUCE(TEMP4), X/*THIS ENSURES THAT THE R.H.S. FO THE D.E. TEMP6 IS ORTHOGONAL TO XTHE SOLUTION OF THE HOMOGENEOUS EQUATION */ X TEMP6:PROJECT(TEMP5), X TEMP7:ODE2(TEMP6,ZW,X), X/*SATISFY BOUNDARY CONDITIONS*/ X IF BCL*BCR=1 THEN TEMP8:BC2(TEMP7,X=0,ZW=0,X=LEN,ZW=0) X ELSE X TEMP8:NEUMANNBC(TEMP7), X/*MAKE SURE THAT = 0*/ X TEMP9:PROJECT(TEMP8), X ADDBASE:['DIFF(W,AMP,I,LAM,J)=RHS(TEMP9)], X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X XPROJECT(EXP):=( X PRO1:EV(EXP,ZW=0), X PRO2:INTEGRATE(PRO1*CFUN,X,0,LEN)/NORM, X EXPAND(EXP-PRO2*CFUN))$ XNEUMANNBC(EXP):=( X REXP:RHS(EXP), X NBC1:DIFF(REXP,X), X IF BCL=1 AND BCR=2 THEN NBC2:SOLVE([EV(REXP,X=0),EV(NBC1,X=LEN)], X [%K1,%K2]), X IF BCL=2 AND BCR=1 THEN NBC2:SOLVE([EV(REXP,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X IF BCL=2 AND BCR=2 THEN NBC2:SOLVE([EV(NBC1,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X EV(EXP,NBC2))$ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X X X X X X X X/**********************************************************************/ X/**********************************************************************/ X X X X X/* PROGRAM NUMBER 11: REDUCTION2(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN SYSTEMS OF ORDINARY DIFFERENTIAL X EQUATIONS. SEE PAGE 176 AND PAGE 211 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA".*/ X X X X X XREDUCTION2():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X X XSETUP():=(/*THE FUNCTION SETUP NEEDS THE DIMENSION N OF THE SYSTEM OF X O.D.E., THE N-DIMENSIONAL LIST OF VARIABLES Y, THE SYSTEM X OF O.D.E'S CALLED EQS, THE CRITICAL EIGENVECTOR CFUN AND X THE ADJOINT CRITICAL EIGENVECTOR AFUN.*/ X PRINTFLAG:TRUE, X LS_LIST:[], X N:READ("NUMBER OF EQUATIONS"), X Y:MAKELIST(READ("ENTER VARIABLE NUMBER",I),I,1,N), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE",ALPHA), X PRINT("WE DEFINE LAM =",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENVECTOR AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENVECTOR"), X KILL(W,ZWLIST), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W),[AMP,LAM]), X PRINT("ENTER THE DIFFERENTIAL EQUATION"), X EQS:MAKELIST(READ("DIFF(",Y[I],",T)="),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT(EQLAM), X WNULL:VFUN(W,0), X ZWNULL:VFUN(ZWLIST,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X DATABASE:APPEND(DATABASE,MAP(LAMBDA([U],DIFF(U,LAM)=0),W)), X NORM:AFUN.CFUN, X TEMP1:EV(EQLAM,SUB,DIFF), X PRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), X NULLANS:READ(" ARE ZERO, Y/N") X )$ X X XG_POLY(I,J):=BLOCK(/*PROVIDED ALL NECESSARY KNOWLEDGE ABOUT THE TAYLOR X COEFFICIENTS OF W(AMP,LAM) IS STORED IN THE LIST X DATABASE, G_POLY CALCULATES ANY SPECIFIC TAYLOR X COEFFICIENT OF THE BIFURCATION EQUATION X G(AMP,LAM)= 0 VIA THE PROJECTION ONTO CFUN.*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X BIFEQ[I,J]:RATSIMP(TEMP4.AFUN))$ X X X XW_POLY(I,J):=BLOCK(/*THE FUNCTION W_POLY ALLOWS TO CALCULATE ITERATIVELY, X I.E. STARTING WITH THE LOWEST ORDER, ALL TAYLOR X COEFFICIENTS OF W(AMP,LAM) BY PROJECTING ONTO THE X COMPLEMENT OF CFUN*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL), X TEMP5:EV(TEMP4,ZWNULL).AFUN, X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:SOLVE(TEMP6,ZWLIST), X TEMP8:EV(ZWLIST,TEMP7), X TEMP9:RATSIMP(TEMP8.AFUN), X IF TEMP9 = 0 THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),TEMP8) X ELSE X (TEMP10:RATSIMP(TEMP8 - TEMP9/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J), X U,1,N),TEMP10)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X X X X X X X X/**************************************************************************/ X/**************************************************************************/ X X X X/* PROGRAM NUMBER 12: A LIAPUNOV-SCHMIDT REDUCTION FOR STEADY STATE X BIFURCATIONS IN SYSTEMS OF PARTIAL DIFFERENTIAL EQUATIONS X DEPENDING ON ONE INDEPENDENT SPACE VARIABLE. SEE PAGE 189 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". X*/ X X X X X/*THIS FILE CONTAINS REDUCTION3(), A FUNCTION TO PERFORM A LIAPUNOV- XSCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS FROM N COUPLED XPARTIAL DIFFERENTIAL EQUATIONS ON ONE SPATIAL DIMENSION. XTHE FOLLOWING ADDITIONAL FUNCTIONS ARE SUPPLIED: X SETUP() ALLOWS THE PROBLEM TO BE ENTERED. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM). X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X SOLVE_ODE(EXP) SOLVES CERTAIN ORDINARY DIFFERENTIAL EQUATIONS VIA A FOURIER X MODE ANSATZ. X FEEDW(EXP) ENSURES THAT = 0 . X FIND_TRIG(EXP) IDENTIFIES FOURIER MODES. X SETIFY(LIST) TRANSFORMS A LIST INTO A SET. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE, ...] X DIFFNULL(I,J) SETS THE DERIVATIVE DIFF(W,AMP,I,LAM,J) TO ZERO.*/ X XREDUCTION3():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X XSETUP():=( /* THIS FUNCTION PERFORMS THE INPUT FOR THE LIAPUNOV-SCHMIDT X REDUCTION*/ XASSUME_POS:TRUE, XLS_LIST:[], XN:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), XY:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), XXVAR:READ("ENTER THE SPATIAL COORDINATE"), XALPHA:READ("ENTER THE BIFURCATION PARAMETER"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), XPRINT("WE DEFINE LAM = ",ALPHA - CAL), XCFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), XAFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS A LIST"), XKILL(W), XW:MAKELIST(CONCAT(W,I),I,1,N), XZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), XNULLIST:MAKELIST(0,I,1,N), XDEPENDS(APPEND(ZWLIST,W,Y),CONS(XVAR,[AMP,LAM])), XEQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), XEQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), XPRINT(EQLAM), XLEN:READ("WHAT IS THE LENGTH OF THE SPACE INTERVAL"), XWNULL:VFUN(W,0), XSUB:MAPLIST("=",Y,AMP*CFUN+W), XDATABASE:APPEND(DIFNULL(1,0),DIFNULL(0,1)), XZWNULL:VFUN(ZWLIST,0), XNORM:INTEGRATE(AFUN.CFUN,XVAR,0,LEN), XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICIENTS ARE 0"), XNULLANS:READ("Y,N") X)$ X X XG_POLY(I,J):=BLOCK( X/*THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) = 0. IT REQUIRES KOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS G_POLY(",I,",",J,")IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*SET THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZERO*/ X TEMP3:SUBST(DIFNULL(I,J),TEMP2), X /*ENTER ALL INFORMATION IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE)), X /*PROJECT ONTO AFUN*/ X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ[I,J]:INTEGRATE(TRIGREDUCE(TEMP5),XVAR,0,LEN))$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). IT RETURNS A DIFFERENTIAL EQUATION XFOR THE PARTICULAR COEFFICIENT OF INTEREST (CALLED ZW1,ZW2...). THIS D.E. XIS SOLVED VIA SOLVE_ODE AND THE RESULT IS FED INTO DATABASE FROM FEEDW*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W(AMP,",I,",LAM,",J,") IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN ( X ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*RENAME THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZW */ X TEMP3:SUBST(WMAX_DIFF,TEMP2), X /*ENTER ALL INFORMATION IN STORED IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE), X /*THIS IS THE PROJECTION Q ONTO THE RANGE OF THE X LINEAR DIFFERENTIAL OPERATOR IN THE PROBLEM*/ X TEMP5:INTEGRATE(EV(TEMP4,ZWNULL).AFUN,XVAR,0,LEN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:TRIGREDUCE(TEMP6), X /*THE SET OF O.D.E.'S TO SOLVE */ X W_DE:EXPAND(TEMP7), X TEMP8:EV(W_DE,VFUN(ZWLIST,0)), X /*IF THE PARTICULAR SOLUTION OF W_DE IS ZERO THEN W=0 */ X IF NULLIST=TEMP8 THEN ( ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE)), X TEMP9:SOLVE_ODE(TEMP8), X FEEDW(TEMP9) X )$ X X XSOLVE_ODE(EXP):=( X/*THIS FUNCTION SOLVE THE D.E. W_DE BY A FOURIER MODE ANSATZ*/ X TRIGFUN:[], X CONST:FALSE, X FOR I THRU N DO ( X /*DETERMINE THE FOURIER MODES*/ X TRIG1:EXP[I], X IF TRIG1 # 0 THEN ( X TRIG2:APPLY1(TRIG1,SINNULL,COSNULL), X IF TRIG2 # 0 THEN ( X CONST:TRUE, X TRIG1:TRIG1-TRIG2), X TRIGFUN:APPEND(FIND_TRIG(TRIG1),TRIGFUN))), X SOL1:DELETE(DUMMY,SETIFY(TRIGFUN)), X /*MAKE AN ANSATZ*/ X ANSATZ:MAKELIST(SUM(AM[I,J]*SOL1[I],I,1,LENGTH(SOL1)),J,1,N), X SOL2:EV(W_DE,MAP("=",ZWLIST,ANSATZ),DIFF), X SOL3:MAKELIST(RATCOEF(SOL2,I),I,SOL1), X EQLIST:[], X FOR I THRU LENGTH(SOL3) DO EQLIST:APPEND(EQLIST,SOL3[I]), X VARLIST:[], X FOR I THRU N DO X FOR J THRU LENGTH(SOL1) DO VARLIST:CONS(AM[J,I],VARLIST), X /*FIND THE AMPLITUDES OF THE FOURIER MODES*/ X SOL4:SOLVE(EQLIST,VARLIST), X /*SOLVE FOR THE CONSTANT FOURIER MODE IF NECESSARY*/ X CANSATZ:0, X IF CONST = TRUE THEN (CANSATZ:MAKELIST(CONCAT(CON,I),I,1,N), X CSOL1:EV(W_DE,MAP("=",ZWLIST,CANSATZ),DIFF), X CSOL2:APPLY1(CSOL1,SINNULL,COSNULL), X CSOL3:SOLVE(CSOL2,CANSATZ)), X EV(ANSATZ+CANSATZ,SOL4,CSOL3))$ X X XFEEDW(EXP):=( X/* THIS FUNCTION ALLOWS THE RESULT OF THE ODE-SOLVER TO BE ENTERED INTO X THE LIST DATABASE. IT CHECKS FOR ORTHOGONALITY TO THE CRITICAL ADJOINT X EIGENFUNCTION AND REMOVES SOLUTIONS OF THE HOMOGENEOUS EQUATION (I.E. X NONORTHOGONAL TERMS)*/ X F2:INTEGRATE(EXP.AFUN,XVAR,0,LEN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X EXP) X ELSE (ORTHO_RESULT:RATSIMP(EXP- F2/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X/*COLLECT ALL INFORMATION STORED IN BIFEQ AND ASSEMBLE THE BIFURCATION XEQUATION*/ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X XSETIFY(LIST):=(/*TRANSFORMS A LIST INTO A SET, I.E. REMOVES DUPLICATES*/ X SET:[LIST[1]], X FOR I :2 THRU LENGTH(LIST) DO X IF NOT MEMBER(LIST[I],SET) THEN X SET:CONS(LIST[I],SET), X SET)$ X XFIND_TRIG(EXP):=(/*FINDS THE FOURIER MODES*/ X F_A1:ARGS(EXP+DUMMY), X F_A2:APPLY1(F_A1,SINFIND,COSFIND) X )$ X X/*AUXILIARY FUNCTIONS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE)$ XDEFRULE(COSFIND,DUMMY1*COS(DUMMY2),COS(DUMMY2))$ XDEFRULE(SINFIND,DUMMY1*SIN(DUMMY2),SIN(DUMMY2))$ XDEFRULE(COSNULL,COS(DUMMY1),0)$ XDEFRULE(SINNULL,SIN(DUMMY1),0)$ X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ XDIFNULL(I,J):=MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J)=0),W) $ X X X X X X X X X/*************************************************************************/ X/*************************************************************************/ X X X X X/*PROGRAM NUMBER 13: CONTAINS THE BUILDING BLOCKS TO PERFORM STEADY X STATE BIFURCATIONS FOR MORE PARTIAL DIFFERENTIAL EQUATIONS. SEE PAGES X 198-202 IN "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER X ALGEBRA" */ X X X X X/*THE FOLLOWING FUNCTIONS CAN BE USED TO PERFORM A LIAPUNOV-SCHMIDT REDUCTION X FOR THE 2D BENARD PROBLEM */ X X/*SETUP() ALLOWS THE PROBLEM TO BE ENTERED, X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM), X W_POLY(I,J) DETERMINES A DIFFERENTIAL EQUATION, WHOSE SOLUTION IS THE X COEFFICIENT OF AMP^I LAM^J IN W(AMP,LAM), X FEEDW() ALLOWS THIS SOLUTION TO BE ENTERED INTO THE LIST DATABASE, X INT(EXP) FINDS MULTIPLE DEFINITE INTEGRALS EFFECTIVELY, X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE ...] X*/ X XSETUP():=( X N:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), X Y:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), X SPACE:READ("ENTER NUMBER OF SPATIAL COORDINATES"), X XVAR:READ("ENTER THE SPATIAL COORDINATES AS A LIST"), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), X PRINT("WE DEFINE LAM = ",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS ALIST"), X KILL(W), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W,Y),APPEND([AMP,LAM],XVAR)), X EQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT("ENTER THE LOWER LEFT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X LBOUND:READ("[",X[1],"=...,]"), X PRINT("ENTER THE UPPER RIGHT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X UBOUND:READ("[",X[1],"=...,]"), X WNULL:VFUN(W,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X ZWNULL:VFUN(ZWLIST,0), X NORM:INT(AFUN.CFUN), X TEMP1:EV(EQLAM,SUB,DIFF), X EQLAM X )$ X X X XG_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ:INT(TEMP5))$ X X X X XW_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:INT(EV(TEMP4,ZWNULL).AFUN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X FOR I THRU SPACE DO TEMP7:TRIGREDUCE(TEMP6,XVAR[I]), X W_DE:RATSIMP(TEMP7 ), X PRINT("NOW SOLVE THE EQUATIONS",W_DE,"=0! THEY ARE GIVEN IN W_DE"), X PRINT("CALL FEEDW() TO PROCEED"))$ X X X XFEEDW():=( X ADDBASE:[], X RESULT:MAKELIST((PRINT("ENTER",ZWLIST[K]),READ()),K,1,N), X F2:INT(RESULT.AFUN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X RESULT) X ELSE (ORTHO_RESULT:RATSIMP(RESULT- F2/NORM*CFUN), X/*PROJECTION Q*/ X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X ADDBASE)$ X X XINT(EXP):=(%INT:EXP,FOR I THRU SPACE DO %INT:INTEGRATE(TRIGREDUCE(%INT, X XVAR[I]),XVAR[I]), X RATSIMP(EV(%INT,UBOUND)-EV(%INT,LBOUND)))$ X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ END_OF_FILE if test 39932 -ne `wc -c <'program14'`; then echo shar: \"'program14'\" unpacked with wrong size! fi # end of 'program14' fi if test -f 'program15' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program15'\" else echo shar: Extracting \"'program15'\" \(39932 characters\) sed "s/^X//" >'program15' <<'END_OF_FILE' X/* PROGRAM NUMBER 6: TWOVAR(), A FUNCTION TO PERFORM A TWO VARIABLE X EXPANSION TO ORDER EPSILON FOR N WEAKLY COUPLED PERTURBED NONLINEAR X OSCILLATORS. SEE PAGE 94 IN "PERTURBATION METHODS, BIFURCATION X THEORY AND COMPUTER ALGEBRA". */ X X X XTWOVAR():=BLOCK( X XCHOICE:READ("DO YOU WANT TO ENTER NEW DATA (Y/N)"), X XIF CHOICE = N THEN GO(JUMP1), X X/* INPUT D.E.'S */ XN:READ("NUMBER OF D.E.'S"), X XPRINT("THE",N,"D.E.'S WILL BE IN THE FORM:"), XPRINT("X[I]'' + W[I]^2 X[I] = E F[I](X[1],...,X[",N,"],T)"), X XFOR I:1 THRU N DO X X[I]:READ("ENTER SYMBOL FOR X[",I,"]"), X XFOR I:1 THRU N DO X DEPENDS(X[I],T), X XFOR I:1 THRU N DO X W[I]:READ("ENTER W[",I,"]"), X XFOR I:1 THRU N DO X F[I]:READ("ENTER F[",I,"]"), X XJUMP1, X X/* UPDATE EQS FOR SUBSTITUTION OF RESONANT VALUES ON 2ND TIME THRU */ XFOR I:1 THRU N DO( X W[I]:EV(W[I]), X F[I]:EV(F[I])), X X/* ECHO BACK THE D.E.'S */ XPRINT("THE D.E.'S ARE ENTERED AS:"), XFOR I:1 THRU N DO X PRINT(X[I],"'' +",W[I]^2*X[I],"=",E*F[I]), X XPRINT("THE METHOD ASSUMES A SOLUTION IN THE FORM:"), XPRINT("X[I] = X0[I] + E X1[I]"), XPRINT("WHERE X0[I] = A[I](ETA) COS W[I] XI + B[I](ETA) SIN W[I] XI"), XPRINT("WHERE XI = T AND ETA = E T"), X X/* DON'T USE A OR B AS PARAMETERS IN THE GIVEN D.E.'S */ XDEPENDS([A,B],ETA), X XFOR I:1 THRU N DO X X0[I]:A[I]*COS(W[I]*XI)+B[I]*SIN(W[I]*XI), X XFOR I:1 THRU N DO X G[I]:EV(F[I],T=XI), X XFOR I:1 THRU N DO( X FOR J:1 THRU N DO X G[I]:EV(G[I],X[J]::X0[J])), X XFOR I:1 THRU N DO( X G[I]:G[I]-2*DIFF(X0[I],XI,1,ETA,1), X G[I]:EV(G[I],DIFF), X G[I]:EXPAND(TRIGREDUCE(EXPAND(G[I])))), X X/* COLLECT SECULAR TERMS */ XFOR I:1 THRU N DO( X S[I]:COEFF(G[I],SIN(W[I]*XI)), X C[I]:COEFF(G[I],COS(W[I]*XI))), X X/* DISPLAY SECULAR TERMS */ XPRINT("REMOVAL OF SECULAR TERMS IN THE X1[I] EQS. GIVES:"), X XFOR I:1 THRU N DO( X PRINT(S[I],"= 0"), X PRINT(C[I],"= 0")), X XABEQS:APPEND(MAKELIST(S[I],I,1,N),MAKELIST(C[I],I,1,N)), X XCHOICE2:READ("DO YOU WANT TO TRANSFORM TO POLAR COORDINATES (Y/N)"), X XIF CHOICE2=N THEN GO(JUMP2), X X/* TRANSFORM TO POLAR COORDINATES */ XDEPENDS([R,THETA],ETA), XTRANS:MAKELIST([A[I]=R[I]*COS(THETA[I]),B[I]=R[I]*SIN(THETA[I])],I,1,N), XINTEQS:EV(ABEQS,TRANS,DIFF), XRTHEQS:SOLVE(INTEQS,APPEND(MAKELIST(DIFF(R[I],ETA),I,1,N), X MAKELIST(DIFF(THETA[I],ETA),I,1,N))), XRTHEQS:TRIGSIMP(RTHEQS), XRTHEQS:EXPAND(TRIGREDUCE(RTHEQS)), XPRINT(RTHEQS), X XJUMP2, X XCHOICE3:READ("DO YOU WANT TO SEARCH FOR RESONANT PARAMETER VALUES (Y/N)"), X XIF CHOICE3=N THEN GO(END), X X/* OBTAIN FREQUENCIES WHICH APPEAR ON RHS'S OF D.E.'S */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE FREQS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR I:1 THRU N DO( X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(G[I]), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X TEMP3:APPLY1(TEMP2,COSARG,SINARG), X TEMP4:EV(TEMP3,XI=1), X/* REMOVE DUPLICATE ARGUMENTS */ X FREQ[I]:SETIFY(TEMP4)), X X/* DISPLAY FREQUENCIES */ XFOR I:1 THRU N DO( X PRINT(X[I],"EQ'S RESONANT FREQ =",W[I]), X PRINT("FREQS ON RHS =",FREQ[I])), X XJUMP3, X XPAR:READ("WHICH PARAMETER TO SEARCH FOR ?"), X X/* COMPUTE RESONANT VALUES OF DESIRED PARAMETER */ XRESVALS:[], XFOR I:1 THRU N DO( X FOR J:1 THRU LENGTH(FREQ[I]) DO( X RES:SOLVE(PART(FREQ[I],J)=W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES), X RES:SOLVE(PART(FREQ[I],J)=-W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES))), X X/* ELIMINATE DUPLICATE PARAMETER VALUES */ XRESVALUES:SETIFY(RESVALS), X X/* DISPLAY RESONANT PARAMETER VALUES */ XPRINT(RESVALUES), X XCHOICE4:READ("DO YOU WANT TO SEARCH FOR ANOTHER PARAMETER (Y/N) ?"), X XIF CHOICE4=Y THEN GO(JUMP3), X XEND," ")$ X X/* AUXILIARY FUNCTIONS */ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X X X X X X X X X/************************************************************************/ X/************************************************************************/ X X X X/*PROGRAM NUMBER 7: AVERAGE(), ALLOWS TO PERFORM FIRST AND SECOND ORDER X AVERAGING. SEE PAGE 114 IN "PERTURBATION METHODS, BIFURCATION THEORY X AND COMPUTER ALGEBRA". */ X X X X/* PROGRAM TO PERFORM 1ST OR 2ND ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERAGE():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("DO YOU WANT FIRST OR SECOND ORDER AVERAGING? (1/2)"), X COEFVFEPS1:COEFF(VF2,EPS), X COEFVFEPS2:COEFF(VF2,EPS,2), X G1BAR:DEMOIVRE(APPLY1(COEFVFEPS1,KILLEXP)), X IF M=1 THEN RESULT:EPS*G1BAR X ELSE( X G1TILDE:COEFVFEPS1-G1BAR, X W1:INTEGRATE(G1TILDE,T), X REMARRAY(JACOB), X JACOB[I,J] := DIFF(G1TILDE[I],X[J]), X JACG1TILDE:GENMATRIX(JACOB,N,N), X G2BAR:DEMOIVRE(APPLY1(EXPAND(JACG1TILDE.W1)+COEFVFEPS2,KILLEXP)), X RESULT:EPS*G1BAR+EPS^2*G2BAR), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS TO KILL EXPONENTIAL TERMS */ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X X X X X X X X X/*************************************************************************/ X/*************************************************************************/ X X X X/* PROGRAM NUMBER 8: AVERM(), ALLOWS MTH ORDER AVERAGING. SEE PAGE 128 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X/* PROGRAM TO PERFORM MTH ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERM():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("ENTER ORDER OF AVERAGING"), X IF M=1 X/* FIRST ORDER AVERAGING */ X THEN (TEMP0:DEMOIVRE(APPLY1(F,KILLEXP)), X RESULT:TAYLOR(TEMP0,EPS,0,1)) X/* AVERAGING OF ORDER M > 1 */ X ELSE X Y:MAKELIST(CONCAT('Y,I),I,1,N), X WLIST:MAKELIST(CONCAT('W,I),I,1,N), X DEPENDS(WLIST,CONS(T,Y)), X TRAFO:Y, X XSUB:MAPLIST("=",Y,X), X/* WNULL SETS WLIST TO ZERO */ X WNULL:VFUN(WLIST,0), X JACOB[I,J] := DIFF(WLIST[I],Y[J]), X JAC:GENMATRIX(JACOB,N,N), X/* MAIN LOOP */ X FOR I :1 THRU M-1 DO ( X TEMP1:MAPLIST("=",X,Y+EPS^I*WLIST), X/* HERE X IS THE LIST [X1,X2,...,XN], Y IS THE LIST [Y1,Y2,...,YN], X WLIST IS THE TRANSFORMATION VECTOR [W1,W2,...,WN] */ X TEMP2:TAYLOR(SUM((-EPS)^(I*J)*JAC^^J,J,0,M-1). X (EV(VF2,TEMP1) - DIFF(WLIST,T)*EPS^I),EPS,0,M), X/* JAC IS THE JACOBIAN D WLIST/DY OF THE TRANSFORMATION WLIST */ X TEMP3:MATTOLIST(TEMP2,N), X TEMP4:MAP(EXPAND,TAYLOR(EV(TEMP3,WNULL,DIFF),EPS,0,I)), X/* THE ITH ORDER MEAN */ X MEAN:APPLY1(TEMP4,KILLEXP), X TEMP6:EXPAND(TEMP4-MEAN), X TEMP7:EV(INTEGRATE(TEMP6,T),EPS=1), X/* THE ITH ORDER TRANSFORMATION */ X TEMP8:MAPLIST("=",WLIST,TEMP7), X TEMP9:RATSIMP(TEMP8), X/* THE TRANSFORMED DE */ X VF2:EXPAND(EV(TEMP3,TEMP9,DIFF,XSUB,INFEVAL)), X/* THE SUM OF ALL TRANSFORMATIONS */ X TRAFO:EXPAND(TRAFO+EV(EPS^I*WLIST,TEMP9))), X/* END OF MAIN LOOP */ X PRINT("THE TRANSFORMATION: ",X,"="), X PRINT(RATSIMP(DEMOIVRE(TRAFO))), X/* THE FINAL AVERAGING */ X RESULT:APPLY1(VF2,KILLEXP), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS */ X VFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X MATTOLIST(MAT,DIM):=IF DIM>1 THEN MAKELIST(MAT[I,1],I,1,DIM) ELSE [MAT]$ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X X X X X X X X X/***********************************************************************/ X/***********************************************************************/ X X X X/* PROGRAM NUMBER 9: LIE(), LIE TRANSFORMATIONS FOR HAMILTONIAN SYSTEMS. X SEE PAGE 144 IN "PERTURBATION METHODS, BIFURCATION THEORY AND X COMPUTER ALGEBRA". */ X X XLIE():=BLOCK( X X/* INPUT PROBLEM ? */ XCHOICE1:READ("DO YOU WANT TO INPUT A NEW PROBLEM (Y/N) ?"), X XIF CHOICE1=N THEN GO(JUMP1), X X/* INPUT PROBLEM */ XN:READ("ENTER NUMBER OF DEGREES OF FREEDOM"), XFOR II:1 THRU N DO ( X Q[II]:READ("ENTER SYMBOL FOR Q[",II,"]"), X P[II]:READ("ENTER SYMBOL FOR P[",II,"]")), XKILL(W), XPRINT("THE HAMILTONIAN DEPENDS ON THE Q'S, P'S, T AND E (SMALL PARAMETER)"), XPRINT("THE E=0 PROBLEM MUST BE OF THE FORM:"), XPRINT("H =",SUM((P[II]^2+W[II]^2*Q[II]^2)/2,II,1,N)), XHORIGINAL:READ("ENTER THE HAMILTONIAN"), XPRINT("H =",HORIGINAL), X X/* TRANSFORM TO ACTION-ANGLE VARIABLES */ X/* FIND THE W[II]'S */ XH0:EV(HORIGINAL,E=0), XFOR II:1 THRU N DO X W[II]:SQRT(DIFF(H0,Q[II],2)), XPRINT("THE ACTION-ANGLE VARIABLES ARE J'S FOR ACTION, PHI'S FOR ANGLE"), XFOR II:1 THRU N DO X TR[II]:[Q[II]=SQRT(2*J[II]/W[II])*SIN(PHI[II]), X P[II]=SQRT(2*J[II]*W[II])*COS(PHI[II])], XTRAN:MAKELIST(TR[II],II,1,N), XH:EV(HORIGINAL,TRAN,ASSUME_POS:TRUE,INFEVAL), XH:TRIGSIMP(H), XH:EXPAND(TRIGREDUCE(EXPAND(H))), XPRINT("H =",H), X XJUMP1, X X/* INPUT TRUNCATION ORDER */ XNTRUNC:READ("ENTER HIGHEST ORDER TERMS IN E TO BE KEPT"), XFOR II:0 THRU NTRUNC DO X H[II]:RATCOEF(H,E,II), X X/* LIE TRANSFORMS */ X/* NEAR IDENTITY TRANSFORMATION FROM (J,PHI)'S TO (I,PSI)'S */ X X/* UPDATE VARIABLES */ XFOR II:1 THRU N DO( X P[II]:I[II], X Q[II]:PSI[II]), X X/* REPLACE J AND PHI BY I AND PSI IN H'S */ XFOR II:0 THRU NTRUNC DO X H[II]:EV(H[II],MAKELIST(J[III]=I[III],III,1,N), X MAKELIST(PHI[III]=PSI[III],III,1,N)), X XK[0]:H[0], X X/* DECLARE WGEN[I] TO BE A FN OF T, Q'S AND P'S */ XKILL(WGEN), XDEPENDS(WGEN,[T]), XFOR II:1 THRU N DO X DEPENDS(WGEN,[Q[II],P[II]]), X X/* E=0 PROBLEM IS OF FORM SUM(W[II]*I[II]) */ X/* CHOOSE WGEN[II] TO KILL AS MUCH AS POSSIBLE IN EQ(II) */ X/* EQUATE K[II] TO UNREMOVABLE TERMS */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE ARGS OF TRIG TERMS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR II:1 THRU NTRUNC DO( X X EQN[II]:EXPAND(TRIGREDUCE(EXPAND(EQ(II)))), X TEMP:EXPAND(EV(EQN[II],WGEN[II]=0)), X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(TEMP), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X ARG1:APPLY1(TEMP2,COSARG,SINARG), X/* REMOVE DUPLICATE ARGUMENTS */ X ARG2:SETIFY(ARG1), X/* REMOVE RESONANT ARGUMENTS */ X ARG3:SUBLIST(ARG2,NOTRESP), X X/* CHOOSE WGEN TO ELIMINATE NONRESONANT TERMS */ X LENG:LENGTH(ARG3), X WGENTEMP1:0, X FOR JJ:1 THRU LENG DO( X WGENTEMP2:AAA*COS(PART(ARG3,JJ))+BBB*SIN(PART(ARG3,JJ)), X TEMP4:EV(EQN[II],WGEN[II]=WGENTEMP2,DIFF), X TEMP5:SOLVE([RATCOEF(TEMP4,COS(PART(ARG3,JJ))), X RATCOEF(TEMP4,SIN(PART(ARG3,JJ)))],[AAA,BBB]), X WGENTEMP1:WGENTEMP1+EV(WGENTEMP2,TEMP5)), X WGEN[II]:WGENTEMP1, X PRINT("WGEN[",II,"] ="), X PRINT(WGEN[II]), X K[II]:EXPAND(EV(EQN[II],DIFF)), X K[II]:EXPAND(RATSIMP(K[II])), X PRINT("THE TRANSFORMED HAMILTONIAN K[",II,"] ="), X PRINT(K[II])), X XKAMILTONIAN:SUM(K[II]*E^II,II,0,NTRUNC), XPRINT ("THE TRANSFORMED HAMILTONIAN IS "), XPRINT ("K =",KAMILTONIAN), X XCHOICE2:READ("DO YOU WANT TO SEE THE NEAR IDENTITY TRANSFORMATION (Y/N) ?"), XIF CHOICE2=N THEN GO(END), X X/* THE NEAR IDENTITY TRANSFORMATION */ XFOR II:1 THRU N DO( X JTRANS[II]:SUM(S(III,P[II])*E^III,III,0,NTRUNC), X PHITRANS[II]:SUM(S(III,Q[II])*E^III,III,0,NTRUNC)), XFOR II:1 THRU N DO( X PRINT(J[II],"=",JTRANS[II]), X PRINT(PHI[II],"=",PHITRANS[II])), X XEND, XKAMILTONIAN)$ X X XPOISSON(F,G):= X SUM(DIFF(F,Q[II])*DIFF(G,P[II])-DIFF(F,P[II])*DIFF(G,Q[II]),II,1,N)$ X XL(II,F):=POISSON(WGEN[II],F)$ X XS(II,F):=(IF II=0 THEN F ELSE SUM(L(II-M,S(M,F)),M,0,II-1)/II)$ X XEQ(II):=(H[II]+(DIFF(WGEN[II],T)+POISSON(WGEN[II],H[0]))/II X +SUM(L(II-M,K[M])+M*S(II-M,H[M]),M,1,II-1)/II)$ X XLZAP(ANY):=DIFF(ANY,T)+POISSON(ANY,H[0])$ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XNOTRESP(EXP):=NOT IS(LZAP(EXP) = 0)$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X X X X X X X X X X/**************************************************************************/ X/**************************************************************************/ X X X/*PROGRAM NUMBER 10: REDUCTION1(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN ONE DIFFERENTIAL EQUATION DEPENDING X ON ONE INDEPENDENT VARIABLE. SEE PAGE 170 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X X X/*THIS FILE CONTAINS REDUCTION1(), A FUNCTION TO PERFORM A LIAPUNOV- X SCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS OF NONLINEAR D.E.'S X OF THE FORM Y'' + F(Y,Y',ALPHA) = 0. Y = Y(X) IS DEFINED ON A REAL X INTERVAL WITH DIRICHLET OR NEUMANN BOUNDARY CONDITIONS AND F DEPENDS X ONLY LINEARLY ON ALPHA. X IT ALSO CONTAINS THESE ADDITIONAL FUNCTIONS: X SETUP() ALLOWS TO ENTER THE PROBLEM. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION. X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X PROJECT(EXP) ENSURES THAT =0. X NEUMANNBC(EXP) ACCOUNTS FOR NEUMANN BOUNDARY CONDITIONS. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. */ X X XREDUCTION1():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X XSETUP():=( X/*THE FUNCTION SETUP ASKS FOR THE VARIABLES OF THE D.E., THE BIFURCATION XPOINT, THE CRITICAL EIGENFUNCTION, THE X-INTERVAL, THE BOUNDARY XCONDITIONS AND THE DIFFERENTIAL EQUATION. */ XASSUME_POS:TRUE, XLS_LIST:[], XY:READ("ENTER DEPENDENT VARIABLE"), XPRINT("USE X AS THE INDEPENDENT VARIABLE AND ALPHA AS A PARAMETER TO VARY"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE ALPHA"), XPRINT("WE DEFINE LAM = ALPHA - ",CAL), XCFUN: READ("ENTER THE CRITICAL EIGENFUNCTION"), XDEPENDS([ZW,Y,W],[AMP,LAM,X]), XLEN:READ("WHAT IS THE LENGTH OF THE X-INTERVAL"), XNORM:INTEGRATE(CFUN^2,X,0,LEN), XPRINT("SPECIFY THE BOUNDARY CONDITIONS"), XPRINT("YOUR CHOICE FOR THE B.C. ON Y AT X=0 AND X=",LEN), XPRINT("ENTER 1 FOR Y=0, 2 FOR Y'=0"), XBCL:READ("B.C. AT 0?"), XBCR:READ("B.C. AT",LEN,"?"), XEQ:DIFF(Y,X,2) + X READ("THE D.E. IS OF THE FORM Y'' + F(Y,Y',ALPHA) = 0,ENTER F"), XEQLAM:EV(EQ,ALPHA=LAM+CAL), XPRINT(EQLAM), XDATABASE:[DIFF(W,AMP)=0,DIFF(W,LAM)=0], XSUB:Y=AMP*CFUN+W, XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), XNULLANS:READ(" ARE ZERO, Y/N") X)$ X X XG_POLY(I,J):=BLOCK( X/* THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) =0. IT REQUIRES KNOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE */ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, X Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*THIS DIFFERENTIAL OF W WILL BE ANNIHILATED BY THE PROJECTION ONTO THE XCRITICAL EIGENFUNCTION, HENCE WE SET IT TO ZERO ALREADY HERE */ X TEMP3:SUBST('DIFF(W,AMP,I,LAM,J)=0,TEMP2), X/* HERE WE INSERT ALL KNOWLEDGE GAINED THROUGH W_POLY */ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X/*PROJECTION ONTO CFUN, THE CRITICAL EIGENFUNCTION */ X TEMP5:TRIGREDUCE(CFUN*TEMP4), X BIFEQ[I,J]:RATSIMP(INTEGRATE(TEMP5,X,0,LEN)) X )$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). THE RESULT IS FED INTO DATABASE.*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*WE SUBSTITUTE ZW FOR THE UNKNOWN TAYLOR COEFFICIENT AND SOLVE XAN O.D.E. FOR ZW IN TEMP7 */ X TEMP3:SUBST(DIFF(W,AMP,I,LAM,J)=ZW,TEMP2), X/*NOW WE INSERT ALL KNOWLEDGE GAINED THROUGH PREVIOUS ITERATIONS*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X TEMP5:TRIGREDUCE(TEMP4), X/*THIS ENSURES THAT THE R.H.S. FO THE D.E. TEMP6 IS ORTHOGONAL TO XTHE SOLUTION OF THE HOMOGENEOUS EQUATION */ X TEMP6:PROJECT(TEMP5), X TEMP7:ODE2(TEMP6,ZW,X), X/*SATISFY BOUNDARY CONDITIONS*/ X IF BCL*BCR=1 THEN TEMP8:BC2(TEMP7,X=0,ZW=0,X=LEN,ZW=0) X ELSE X TEMP8:NEUMANNBC(TEMP7), X/*MAKE SURE THAT = 0*/ X TEMP9:PROJECT(TEMP8), X ADDBASE:['DIFF(W,AMP,I,LAM,J)=RHS(TEMP9)], X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X XPROJECT(EXP):=( X PRO1:EV(EXP,ZW=0), X PRO2:INTEGRATE(PRO1*CFUN,X,0,LEN)/NORM, X EXPAND(EXP-PRO2*CFUN))$ XNEUMANNBC(EXP):=( X REXP:RHS(EXP), X NBC1:DIFF(REXP,X), X IF BCL=1 AND BCR=2 THEN NBC2:SOLVE([EV(REXP,X=0),EV(NBC1,X=LEN)], X [%K1,%K2]), X IF BCL=2 AND BCR=1 THEN NBC2:SOLVE([EV(REXP,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X IF BCL=2 AND BCR=2 THEN NBC2:SOLVE([EV(NBC1,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X EV(EXP,NBC2))$ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X X X X X X X X/**********************************************************************/ X/**********************************************************************/ X X X X X/* PROGRAM NUMBER 11: REDUCTION2(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN SYSTEMS OF ORDINARY DIFFERENTIAL X EQUATIONS. SEE PAGE 176 AND PAGE 211 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA".*/ X X X X X XREDUCTION2():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X X XSETUP():=(/*THE FUNCTION SETUP NEEDS THE DIMENSION N OF THE SYSTEM OF X O.D.E., THE N-DIMENSIONAL LIST OF VARIABLES Y, THE SYSTEM X OF O.D.E'S CALLED EQS, THE CRITICAL EIGENVECTOR CFUN AND X THE ADJOINT CRITICAL EIGENVECTOR AFUN.*/ X PRINTFLAG:TRUE, X LS_LIST:[], X N:READ("NUMBER OF EQUATIONS"), X Y:MAKELIST(READ("ENTER VARIABLE NUMBER",I),I,1,N), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE",ALPHA), X PRINT("WE DEFINE LAM =",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENVECTOR AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENVECTOR"), X KILL(W,ZWLIST), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W),[AMP,LAM]), X PRINT("ENTER THE DIFFERENTIAL EQUATION"), X EQS:MAKELIST(READ("DIFF(",Y[I],",T)="),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT(EQLAM), X WNULL:VFUN(W,0), X ZWNULL:VFUN(ZWLIST,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X DATABASE:APPEND(DATABASE,MAP(LAMBDA([U],DIFF(U,LAM)=0),W)), X NORM:AFUN.CFUN, X TEMP1:EV(EQLAM,SUB,DIFF), X PRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), X NULLANS:READ(" ARE ZERO, Y/N") X )$ X X XG_POLY(I,J):=BLOCK(/*PROVIDED ALL NECESSARY KNOWLEDGE ABOUT THE TAYLOR X COEFFICIENTS OF W(AMP,LAM) IS STORED IN THE LIST X DATABASE, G_POLY CALCULATES ANY SPECIFIC TAYLOR X COEFFICIENT OF THE BIFURCATION EQUATION X G(AMP,LAM)= 0 VIA THE PROJECTION ONTO CFUN.*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X BIFEQ[I,J]:RATSIMP(TEMP4.AFUN))$ X X X XW_POLY(I,J):=BLOCK(/*THE FUNCTION W_POLY ALLOWS TO CALCULATE ITERATIVELY, X I.E. STARTING WITH THE LOWEST ORDER, ALL TAYLOR X COEFFICIENTS OF W(AMP,LAM) BY PROJECTING ONTO THE X COMPLEMENT OF CFUN*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL), X TEMP5:EV(TEMP4,ZWNULL).AFUN, X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:SOLVE(TEMP6,ZWLIST), X TEMP8:EV(ZWLIST,TEMP7), X TEMP9:RATSIMP(TEMP8.AFUN), X IF TEMP9 = 0 THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),TEMP8) X ELSE X (TEMP10:RATSIMP(TEMP8 - TEMP9/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J), X U,1,N),TEMP10)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X X X X X X X X/**************************************************************************/ X/**************************************************************************/ X X X X/* PROGRAM NUMBER 12: A LIAPUNOV-SCHMIDT REDUCTION FOR STEADY STATE X BIFURCATIONS IN SYSTEMS OF PARTIAL DIFFERENTIAL EQUATIONS X DEPENDING ON ONE INDEPENDENT SPACE VARIABLE. SEE PAGE 189 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". X*/ X X X X X/*THIS FILE CONTAINS REDUCTION3(), A FUNCTION TO PERFORM A LIAPUNOV- XSCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS FROM N COUPLED XPARTIAL DIFFERENTIAL EQUATIONS ON ONE SPATIAL DIMENSION. XTHE FOLLOWING ADDITIONAL FUNCTIONS ARE SUPPLIED: X SETUP() ALLOWS THE PROBLEM TO BE ENTERED. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM). X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X SOLVE_ODE(EXP) SOLVES CERTAIN ORDINARY DIFFERENTIAL EQUATIONS VIA A FOURIER X MODE ANSATZ. X FEEDW(EXP) ENSURES THAT = 0 . X FIND_TRIG(EXP) IDENTIFIES FOURIER MODES. X SETIFY(LIST) TRANSFORMS A LIST INTO A SET. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE, ...] X DIFFNULL(I,J) SETS THE DERIVATIVE DIFF(W,AMP,I,LAM,J) TO ZERO.*/ X XREDUCTION3():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X XSETUP():=( /* THIS FUNCTION PERFORMS THE INPUT FOR THE LIAPUNOV-SCHMIDT X REDUCTION*/ XASSUME_POS:TRUE, XLS_LIST:[], XN:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), XY:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), XXVAR:READ("ENTER THE SPATIAL COORDINATE"), XALPHA:READ("ENTER THE BIFURCATION PARAMETER"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), XPRINT("WE DEFINE LAM = ",ALPHA - CAL), XCFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), XAFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS A LIST"), XKILL(W), XW:MAKELIST(CONCAT(W,I),I,1,N), XZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), XNULLIST:MAKELIST(0,I,1,N), XDEPENDS(APPEND(ZWLIST,W,Y),CONS(XVAR,[AMP,LAM])), XEQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), XEQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), XPRINT(EQLAM), XLEN:READ("WHAT IS THE LENGTH OF THE SPACE INTERVAL"), XWNULL:VFUN(W,0), XSUB:MAPLIST("=",Y,AMP*CFUN+W), XDATABASE:APPEND(DIFNULL(1,0),DIFNULL(0,1)), XZWNULL:VFUN(ZWLIST,0), XNORM:INTEGRATE(AFUN.CFUN,XVAR,0,LEN), XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICIENTS ARE 0"), XNULLANS:READ("Y,N") X)$ X X XG_POLY(I,J):=BLOCK( X/*THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) = 0. IT REQUIRES KOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS G_POLY(",I,",",J,")IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*SET THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZERO*/ X TEMP3:SUBST(DIFNULL(I,J),TEMP2), X /*ENTER ALL INFORMATION IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE)), X /*PROJECT ONTO AFUN*/ X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ[I,J]:INTEGRATE(TRIGREDUCE(TEMP5),XVAR,0,LEN))$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). IT RETURNS A DIFFERENTIAL EQUATION XFOR THE PARTICULAR COEFFICIENT OF INTEREST (CALLED ZW1,ZW2...). THIS D.E. XIS SOLVED VIA SOLVE_ODE AND THE RESULT IS FED INTO DATABASE FROM FEEDW*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W(AMP,",I,",LAM,",J,") IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN ( X ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*RENAME THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZW */ X TEMP3:SUBST(WMAX_DIFF,TEMP2), X /*ENTER ALL INFORMATION IN STORED IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE), X /*THIS IS THE PROJECTION Q ONTO THE RANGE OF THE X LINEAR DIFFERENTIAL OPERATOR IN THE PROBLEM*/ X TEMP5:INTEGRATE(EV(TEMP4,ZWNULL).AFUN,XVAR,0,LEN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:TRIGREDUCE(TEMP6), X /*THE SET OF O.D.E.'S TO SOLVE */ X W_DE:EXPAND(TEMP7), X TEMP8:EV(W_DE,VFUN(ZWLIST,0)), X /*IF THE PARTICULAR SOLUTION OF W_DE IS ZERO THEN W=0 */ X IF NULLIST=TEMP8 THEN ( ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE)), X TEMP9:SOLVE_ODE(TEMP8), X FEEDW(TEMP9) X )$ X X XSOLVE_ODE(EXP):=( X/*THIS FUNCTION SOLVE THE D.E. W_DE BY A FOURIER MODE ANSATZ*/ X TRIGFUN:[], X CONST:FALSE, X FOR I THRU N DO ( X /*DETERMINE THE FOURIER MODES*/ X TRIG1:EXP[I], X IF TRIG1 # 0 THEN ( X TRIG2:APPLY1(TRIG1,SINNULL,COSNULL), X IF TRIG2 # 0 THEN ( X CONST:TRUE, X TRIG1:TRIG1-TRIG2), X TRIGFUN:APPEND(FIND_TRIG(TRIG1),TRIGFUN))), X SOL1:DELETE(DUMMY,SETIFY(TRIGFUN)), X /*MAKE AN ANSATZ*/ X ANSATZ:MAKELIST(SUM(AM[I,J]*SOL1[I],I,1,LENGTH(SOL1)),J,1,N), X SOL2:EV(W_DE,MAP("=",ZWLIST,ANSATZ),DIFF), X SOL3:MAKELIST(RATCOEF(SOL2,I),I,SOL1), X EQLIST:[], X FOR I THRU LENGTH(SOL3) DO EQLIST:APPEND(EQLIST,SOL3[I]), X VARLIST:[], X FOR I THRU N DO X FOR J THRU LENGTH(SOL1) DO VARLIST:CONS(AM[J,I],VARLIST), X /*FIND THE AMPLITUDES OF THE FOURIER MODES*/ X SOL4:SOLVE(EQLIST,VARLIST), X /*SOLVE FOR THE CONSTANT FOURIER MODE IF NECESSARY*/ X CANSATZ:0, X IF CONST = TRUE THEN (CANSATZ:MAKELIST(CONCAT(CON,I),I,1,N), X CSOL1:EV(W_DE,MAP("=",ZWLIST,CANSATZ),DIFF), X CSOL2:APPLY1(CSOL1,SINNULL,COSNULL), X CSOL3:SOLVE(CSOL2,CANSATZ)), X EV(ANSATZ+CANSATZ,SOL4,CSOL3))$ X X XFEEDW(EXP):=( X/* THIS FUNCTION ALLOWS THE RESULT OF THE ODE-SOLVER TO BE ENTERED INTO X THE LIST DATABASE. IT CHECKS FOR ORTHOGONALITY TO THE CRITICAL ADJOINT X EIGENFUNCTION AND REMOVES SOLUTIONS OF THE HOMOGENEOUS EQUATION (I.E. X NONORTHOGONAL TERMS)*/ X F2:INTEGRATE(EXP.AFUN,XVAR,0,LEN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X EXP) X ELSE (ORTHO_RESULT:RATSIMP(EXP- F2/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X/*COLLECT ALL INFORMATION STORED IN BIFEQ AND ASSEMBLE THE BIFURCATION XEQUATION*/ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X XSETIFY(LIST):=(/*TRANSFORMS A LIST INTO A SET, I.E. REMOVES DUPLICATES*/ X SET:[LIST[1]], X FOR I :2 THRU LENGTH(LIST) DO X IF NOT MEMBER(LIST[I],SET) THEN X SET:CONS(LIST[I],SET), X SET)$ X XFIND_TRIG(EXP):=(/*FINDS THE FOURIER MODES*/ X F_A1:ARGS(EXP+DUMMY), X F_A2:APPLY1(F_A1,SINFIND,COSFIND) X )$ X X/*AUXILIARY FUNCTIONS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE)$ XDEFRULE(COSFIND,DUMMY1*COS(DUMMY2),COS(DUMMY2))$ XDEFRULE(SINFIND,DUMMY1*SIN(DUMMY2),SIN(DUMMY2))$ XDEFRULE(COSNULL,COS(DUMMY1),0)$ XDEFRULE(SINNULL,SIN(DUMMY1),0)$ X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ XDIFNULL(I,J):=MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J)=0),W) $ X X X X X X X X X/*************************************************************************/ X/*************************************************************************/ X X X X X/*PROGRAM NUMBER 13: CONTAINS THE BUILDING BLOCKS TO PERFORM STEADY X STATE BIFURCATIONS FOR MORE PARTIAL DIFFERENTIAL EQUATIONS. SEE PAGES X 198-202 IN "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER X ALGEBRA" */ X X X X X/*THE FOLLOWING FUNCTIONS CAN BE USED TO PERFORM A LIAPUNOV-SCHMIDT REDUCTION X FOR THE 2D BENARD PROBLEM */ X X/*SETUP() ALLOWS THE PROBLEM TO BE ENTERED, X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM), X W_POLY(I,J) DETERMINES A DIFFERENTIAL EQUATION, WHOSE SOLUTION IS THE X COEFFICIENT OF AMP^I LAM^J IN W(AMP,LAM), X FEEDW() ALLOWS THIS SOLUTION TO BE ENTERED INTO THE LIST DATABASE, X INT(EXP) FINDS MULTIPLE DEFINITE INTEGRALS EFFECTIVELY, X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE ...] X*/ X XSETUP():=( X N:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), X Y:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), X SPACE:READ("ENTER NUMBER OF SPATIAL COORDINATES"), X XVAR:READ("ENTER THE SPATIAL COORDINATES AS A LIST"), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), X PRINT("WE DEFINE LAM = ",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS ALIST"), X KILL(W), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W,Y),APPEND([AMP,LAM],XVAR)), X EQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT("ENTER THE LOWER LEFT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X LBOUND:READ("[",X[1],"=...,]"), X PRINT("ENTER THE UPPER RIGHT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X UBOUND:READ("[",X[1],"=...,]"), X WNULL:VFUN(W,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X ZWNULL:VFUN(ZWLIST,0), X NORM:INT(AFUN.CFUN), X TEMP1:EV(EQLAM,SUB,DIFF), X EQLAM X )$ X X X XG_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ:INT(TEMP5))$ X X X X XW_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:INT(EV(TEMP4,ZWNULL).AFUN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X FOR I THRU SPACE DO TEMP7:TRIGREDUCE(TEMP6,XVAR[I]), X W_DE:RATSIMP(TEMP7 ), X PRINT("NOW SOLVE THE EQUATIONS",W_DE,"=0! THEY ARE GIVEN IN W_DE"), X PRINT("CALL FEEDW() TO PROCEED"))$ X X X XFEEDW():=( X ADDBASE:[], X RESULT:MAKELIST((PRINT("ENTER",ZWLIST[K]),READ()),K,1,N), X F2:INT(RESULT.AFUN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X RESULT) X ELSE (ORTHO_RESULT:RATSIMP(RESULT- F2/NORM*CFUN), X/*PROJECTION Q*/ X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X ADDBASE)$ X X XINT(EXP):=(%INT:EXP,FOR I THRU SPACE DO %INT:INTEGRATE(TRIGREDUCE(%INT, X XVAR[I]),XVAR[I]), X RATSIMP(EV(%INT,UBOUND)-EV(%INT,LBOUND)))$ X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ END_OF_FILE if test 39932 -ne `wc -c <'program15'`; then echo shar: \"'program15'\" unpacked with wrong size! fi # end of 'program15' fi if test -f 'program16' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'program16'\" else echo shar: Extracting \"'program16'\" \(39932 characters\) sed "s/^X//" >'program16' <<'END_OF_FILE' X/* PROGRAM NUMBER 6: TWOVAR(), A FUNCTION TO PERFORM A TWO VARIABLE X EXPANSION TO ORDER EPSILON FOR N WEAKLY COUPLED PERTURBED NONLINEAR X OSCILLATORS. SEE PAGE 94 IN "PERTURBATION METHODS, BIFURCATION X THEORY AND COMPUTER ALGEBRA". */ X X X XTWOVAR():=BLOCK( X XCHOICE:READ("DO YOU WANT TO ENTER NEW DATA (Y/N)"), X XIF CHOICE = N THEN GO(JUMP1), X X/* INPUT D.E.'S */ XN:READ("NUMBER OF D.E.'S"), X XPRINT("THE",N,"D.E.'S WILL BE IN THE FORM:"), XPRINT("X[I]'' + W[I]^2 X[I] = E F[I](X[1],...,X[",N,"],T)"), X XFOR I:1 THRU N DO X X[I]:READ("ENTER SYMBOL FOR X[",I,"]"), X XFOR I:1 THRU N DO X DEPENDS(X[I],T), X XFOR I:1 THRU N DO X W[I]:READ("ENTER W[",I,"]"), X XFOR I:1 THRU N DO X F[I]:READ("ENTER F[",I,"]"), X XJUMP1, X X/* UPDATE EQS FOR SUBSTITUTION OF RESONANT VALUES ON 2ND TIME THRU */ XFOR I:1 THRU N DO( X W[I]:EV(W[I]), X F[I]:EV(F[I])), X X/* ECHO BACK THE D.E.'S */ XPRINT("THE D.E.'S ARE ENTERED AS:"), XFOR I:1 THRU N DO X PRINT(X[I],"'' +",W[I]^2*X[I],"=",E*F[I]), X XPRINT("THE METHOD ASSUMES A SOLUTION IN THE FORM:"), XPRINT("X[I] = X0[I] + E X1[I]"), XPRINT("WHERE X0[I] = A[I](ETA) COS W[I] XI + B[I](ETA) SIN W[I] XI"), XPRINT("WHERE XI = T AND ETA = E T"), X X/* DON'T USE A OR B AS PARAMETERS IN THE GIVEN D.E.'S */ XDEPENDS([A,B],ETA), X XFOR I:1 THRU N DO X X0[I]:A[I]*COS(W[I]*XI)+B[I]*SIN(W[I]*XI), X XFOR I:1 THRU N DO X G[I]:EV(F[I],T=XI), X XFOR I:1 THRU N DO( X FOR J:1 THRU N DO X G[I]:EV(G[I],X[J]::X0[J])), X XFOR I:1 THRU N DO( X G[I]:G[I]-2*DIFF(X0[I],XI,1,ETA,1), X G[I]:EV(G[I],DIFF), X G[I]:EXPAND(TRIGREDUCE(EXPAND(G[I])))), X X/* COLLECT SECULAR TERMS */ XFOR I:1 THRU N DO( X S[I]:COEFF(G[I],SIN(W[I]*XI)), X C[I]:COEFF(G[I],COS(W[I]*XI))), X X/* DISPLAY SECULAR TERMS */ XPRINT("REMOVAL OF SECULAR TERMS IN THE X1[I] EQS. GIVES:"), X XFOR I:1 THRU N DO( X PRINT(S[I],"= 0"), X PRINT(C[I],"= 0")), X XABEQS:APPEND(MAKELIST(S[I],I,1,N),MAKELIST(C[I],I,1,N)), X XCHOICE2:READ("DO YOU WANT TO TRANSFORM TO POLAR COORDINATES (Y/N)"), X XIF CHOICE2=N THEN GO(JUMP2), X X/* TRANSFORM TO POLAR COORDINATES */ XDEPENDS([R,THETA],ETA), XTRANS:MAKELIST([A[I]=R[I]*COS(THETA[I]),B[I]=R[I]*SIN(THETA[I])],I,1,N), XINTEQS:EV(ABEQS,TRANS,DIFF), XRTHEQS:SOLVE(INTEQS,APPEND(MAKELIST(DIFF(R[I],ETA),I,1,N), X MAKELIST(DIFF(THETA[I],ETA),I,1,N))), XRTHEQS:TRIGSIMP(RTHEQS), XRTHEQS:EXPAND(TRIGREDUCE(RTHEQS)), XPRINT(RTHEQS), X XJUMP2, X XCHOICE3:READ("DO YOU WANT TO SEARCH FOR RESONANT PARAMETER VALUES (Y/N)"), X XIF CHOICE3=N THEN GO(END), X X/* OBTAIN FREQUENCIES WHICH APPEAR ON RHS'S OF D.E.'S */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE FREQS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR I:1 THRU N DO( X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(G[I]), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X TEMP3:APPLY1(TEMP2,COSARG,SINARG), X TEMP4:EV(TEMP3,XI=1), X/* REMOVE DUPLICATE ARGUMENTS */ X FREQ[I]:SETIFY(TEMP4)), X X/* DISPLAY FREQUENCIES */ XFOR I:1 THRU N DO( X PRINT(X[I],"EQ'S RESONANT FREQ =",W[I]), X PRINT("FREQS ON RHS =",FREQ[I])), X XJUMP3, X XPAR:READ("WHICH PARAMETER TO SEARCH FOR ?"), X X/* COMPUTE RESONANT VALUES OF DESIRED PARAMETER */ XRESVALS:[], XFOR I:1 THRU N DO( X FOR J:1 THRU LENGTH(FREQ[I]) DO( X RES:SOLVE(PART(FREQ[I],J)=W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES), X RES:SOLVE(PART(FREQ[I],J)=-W[I],PAR), X IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES))), X X/* ELIMINATE DUPLICATE PARAMETER VALUES */ XRESVALUES:SETIFY(RESVALS), X X/* DISPLAY RESONANT PARAMETER VALUES */ XPRINT(RESVALUES), X XCHOICE4:READ("DO YOU WANT TO SEARCH FOR ANOTHER PARAMETER (Y/N) ?"), X XIF CHOICE4=Y THEN GO(JUMP3), X XEND," ")$ X X/* AUXILIARY FUNCTIONS */ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X X X X X X X X X/************************************************************************/ X/************************************************************************/ X X X X/*PROGRAM NUMBER 7: AVERAGE(), ALLOWS TO PERFORM FIRST AND SECOND ORDER X AVERAGING. SEE PAGE 114 IN "PERTURBATION METHODS, BIFURCATION THEORY X AND COMPUTER ALGEBRA". */ X X X X/* PROGRAM TO PERFORM 1ST OR 2ND ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERAGE():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("DO YOU WANT FIRST OR SECOND ORDER AVERAGING? (1/2)"), X COEFVFEPS1:COEFF(VF2,EPS), X COEFVFEPS2:COEFF(VF2,EPS,2), X G1BAR:DEMOIVRE(APPLY1(COEFVFEPS1,KILLEXP)), X IF M=1 THEN RESULT:EPS*G1BAR X ELSE( X G1TILDE:COEFVFEPS1-G1BAR, X W1:INTEGRATE(G1TILDE,T), X REMARRAY(JACOB), X JACOB[I,J] := DIFF(G1TILDE[I],X[J]), X JACG1TILDE:GENMATRIX(JACOB,N,N), X G2BAR:DEMOIVRE(APPLY1(EXPAND(JACG1TILDE.W1)+COEFVFEPS2,KILLEXP)), X RESULT:EPS*G1BAR+EPS^2*G2BAR), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS TO KILL EXPONENTIAL TERMS */ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X X X X X X X X X/*************************************************************************/ X/*************************************************************************/ X X X X/* PROGRAM NUMBER 8: AVERM(), ALLOWS MTH ORDER AVERAGING. SEE PAGE 128 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X/* PROGRAM TO PERFORM MTH ORDER AVERAGING X ON AN N-DIMENSIONAL SYSTEM OF NONAUTONOMOUS ODE'S */ X X/* AVERAGING IS PERFORMED BY CONVERTING TRIG TERMS TO X COMPLEX EXPONENTIALS, THEN KILLING EXPONENTIALS */ X XAVERM():=BLOCK( X CHOICE1:READ("DO YOU WANT TO ENTER A NEW PROBLEM? (Y/N)"), X IF CHOICE1 = N THEN GO(JUMP1), X KILL(X), X PRINT("ARE YOU CONSIDERING A WEAKLY NONLINEAR OSCILLATOR OF THE FORM"), X CHOICE2:READ("Z'' + OMEGA0^2 Z = EPS F(Z,ZDOT,T) ? (Y/N)"), X IF CHOICE2 = N THEN GO(JUMP2), X X/* ENTER DATA FOR SINGLE OSCILLATOR PROBLEM */ X N:2, X OMEGA0:READ("ENTER OMEGA0"), X F:READ("ENTER F(Z,ZDOT,T)")*EPS, X/* DOES F(Z,ZDOT,T) DEPEND EXPLICITLY ON T? */ X TEST:DIFF(F,T), X IF TEST=0 THEN OMEGA1:OMEGA0 X ELSE( X OMEGA:READ("ENTER THE FORCING FREQUENCY"), X K:READ("ENTER THE ORDER OF THE SUBHARMONIC RESONANCE"), X OMEGA1:OMEGA/K), X/* VAN DER POL TRANSFORMATION */ X ZSUB:[Z=COS(OMEGA1*T)*X1-SIN(OMEGA1*T)*X2, X ZDOT=-OMEGA1*SIN(OMEGA1*T)*X1-OMEGA1*COS(OMEGA1*T)*X2], X/* SUBSTITUTE ZSUB INTO TRANSFORMED EQS */ X/* SCALE TIME SO THAT AVERAGING OCCURS OVER PERIOD 2 PI */ X VF:EV([-1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*SIN(OMEGA1*T), X -1/OMEGA1^2*(EPS*KAPOMEGA/K^2*Z + F)*COS(OMEGA1*T)], X ZSUB,T=TAU/OMEGA1,INFEVAL), X VF:EV(VF,TAU=T), X IF OMEGA1 # OMEGA0 X THEN PRINT("WE WRITE EPS*KAPOMEGA =",OMEGA^2-K^2*OMEGA0^2) X ELSE VF:EV(VF,KAPOMEGA=0), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU 2 DO VF2[I]:MAP(FACTOR,VF2[I]), X X:[X1,X2], X GO(JUMP1), X X XJUMP2, X/* ENTER DATA FOR GENERAL PROBLEM OF N ODE'S */ X N:READ("ENTER NUMBER OF DIFFERENTIAL EQUATIONS"), X X:MAKELIST(CONCAT('X,I),I,1,N), X PRINT("THE ODE'S ARE OF THE FORM:"), X PRINT("DX/DT = EPS F(X,T)"), X PRINT("WHERE X =",X), X PRINT("SCALE TIME T SUCH THAT AVERAGING OCCURS OVER 2 PI"), X VF:MAKELIST(READ("ENTER RHS(",I,")=EPS*...")*EPS,I,1,N), X FOR I:1 THRU N DO PRINT("D",X[I],"/DT =",VF[I]), X VF2:EXPAND(EXPONENTIALIZE(VF)), X FOR I:1 THRU N DO VF2[I]:MAP(FACTOR,VF2[I]), X X XJUMP1, X/* AVERAGING */ X M:READ("ENTER ORDER OF AVERAGING"), X IF M=1 X/* FIRST ORDER AVERAGING */ X THEN (TEMP0:DEMOIVRE(APPLY1(F,KILLEXP)), X RESULT:TAYLOR(TEMP0,EPS,0,1)) X/* AVERAGING OF ORDER M > 1 */ X ELSE X Y:MAKELIST(CONCAT('Y,I),I,1,N), X WLIST:MAKELIST(CONCAT('W,I),I,1,N), X DEPENDS(WLIST,CONS(T,Y)), X TRAFO:Y, X XSUB:MAPLIST("=",Y,X), X/* WNULL SETS WLIST TO ZERO */ X WNULL:VFUN(WLIST,0), X JACOB[I,J] := DIFF(WLIST[I],Y[J]), X JAC:GENMATRIX(JACOB,N,N), X/* MAIN LOOP */ X FOR I :1 THRU M-1 DO ( X TEMP1:MAPLIST("=",X,Y+EPS^I*WLIST), X/* HERE X IS THE LIST [X1,X2,...,XN], Y IS THE LIST [Y1,Y2,...,YN], X WLIST IS THE TRANSFORMATION VECTOR [W1,W2,...,WN] */ X TEMP2:TAYLOR(SUM((-EPS)^(I*J)*JAC^^J,J,0,M-1). X (EV(VF2,TEMP1) - DIFF(WLIST,T)*EPS^I),EPS,0,M), X/* JAC IS THE JACOBIAN D WLIST/DY OF THE TRANSFORMATION WLIST */ X TEMP3:MATTOLIST(TEMP2,N), X TEMP4:MAP(EXPAND,TAYLOR(EV(TEMP3,WNULL,DIFF),EPS,0,I)), X/* THE ITH ORDER MEAN */ X MEAN:APPLY1(TEMP4,KILLEXP), X TEMP6:EXPAND(TEMP4-MEAN), X TEMP7:EV(INTEGRATE(TEMP6,T),EPS=1), X/* THE ITH ORDER TRANSFORMATION */ X TEMP8:MAPLIST("=",WLIST,TEMP7), X TEMP9:RATSIMP(TEMP8), X/* THE TRANSFORMED DE */ X VF2:EXPAND(EV(TEMP3,TEMP9,DIFF,XSUB,INFEVAL)), X/* THE SUM OF ALL TRANSFORMATIONS */ X TRAFO:EXPAND(TRAFO+EV(EPS^I*WLIST,TEMP9))), X/* END OF MAIN LOOP */ X PRINT("THE TRANSFORMATION: ",X,"="), X PRINT(RATSIMP(DEMOIVRE(TRAFO))), X/* THE FINAL AVERAGING */ X RESULT:APPLY1(VF2,KILLEXP), X/* REPLACE X BY Y */ X FOR I:1 THRU N DO RESULT:SUBST(CONCAT('Y,I),CONCAT('X,I),RESULT), X RESULT)$ X X X/* AUXILIARY FUNCTIONS */ X VFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X MATTOLIST(MAT,DIM):=IF DIM>1 THEN MAKELIST(MAT[I,1],I,1,DIM) ELSE [MAT]$ X CONTAINS_T(EXP):= NOT FREEOF(T,EXP)$ X MATCHDECLARE(Q,CONTAINS_T)$ X DEFRULE(KILLEXP,EXP(Q),0)$ X X X X X X X X X/***********************************************************************/ X/***********************************************************************/ X X X X/* PROGRAM NUMBER 9: LIE(), LIE TRANSFORMATIONS FOR HAMILTONIAN SYSTEMS. X SEE PAGE 144 IN "PERTURBATION METHODS, BIFURCATION THEORY AND X COMPUTER ALGEBRA". */ X X XLIE():=BLOCK( X X/* INPUT PROBLEM ? */ XCHOICE1:READ("DO YOU WANT TO INPUT A NEW PROBLEM (Y/N) ?"), X XIF CHOICE1=N THEN GO(JUMP1), X X/* INPUT PROBLEM */ XN:READ("ENTER NUMBER OF DEGREES OF FREEDOM"), XFOR II:1 THRU N DO ( X Q[II]:READ("ENTER SYMBOL FOR Q[",II,"]"), X P[II]:READ("ENTER SYMBOL FOR P[",II,"]")), XKILL(W), XPRINT("THE HAMILTONIAN DEPENDS ON THE Q'S, P'S, T AND E (SMALL PARAMETER)"), XPRINT("THE E=0 PROBLEM MUST BE OF THE FORM:"), XPRINT("H =",SUM((P[II]^2+W[II]^2*Q[II]^2)/2,II,1,N)), XHORIGINAL:READ("ENTER THE HAMILTONIAN"), XPRINT("H =",HORIGINAL), X X/* TRANSFORM TO ACTION-ANGLE VARIABLES */ X/* FIND THE W[II]'S */ XH0:EV(HORIGINAL,E=0), XFOR II:1 THRU N DO X W[II]:SQRT(DIFF(H0,Q[II],2)), XPRINT("THE ACTION-ANGLE VARIABLES ARE J'S FOR ACTION, PHI'S FOR ANGLE"), XFOR II:1 THRU N DO X TR[II]:[Q[II]=SQRT(2*J[II]/W[II])*SIN(PHI[II]), X P[II]=SQRT(2*J[II]*W[II])*COS(PHI[II])], XTRAN:MAKELIST(TR[II],II,1,N), XH:EV(HORIGINAL,TRAN,ASSUME_POS:TRUE,INFEVAL), XH:TRIGSIMP(H), XH:EXPAND(TRIGREDUCE(EXPAND(H))), XPRINT("H =",H), X XJUMP1, X X/* INPUT TRUNCATION ORDER */ XNTRUNC:READ("ENTER HIGHEST ORDER TERMS IN E TO BE KEPT"), XFOR II:0 THRU NTRUNC DO X H[II]:RATCOEF(H,E,II), X X/* LIE TRANSFORMS */ X/* NEAR IDENTITY TRANSFORMATION FROM (J,PHI)'S TO (I,PSI)'S */ X X/* UPDATE VARIABLES */ XFOR II:1 THRU N DO( X P[II]:I[II], X Q[II]:PSI[II]), X X/* REPLACE J AND PHI BY I AND PSI IN H'S */ XFOR II:0 THRU NTRUNC DO X H[II]:EV(H[II],MAKELIST(J[III]=I[III],III,1,N), X MAKELIST(PHI[III]=PSI[III],III,1,N)), X XK[0]:H[0], X X/* DECLARE WGEN[I] TO BE A FN OF T, Q'S AND P'S */ XKILL(WGEN), XDEPENDS(WGEN,[T]), XFOR II:1 THRU N DO X DEPENDS(WGEN,[Q[II],P[II]]), X X/* E=0 PROBLEM IS OF FORM SUM(W[II]*I[II]) */ X/* CHOOSE WGEN[II] TO KILL AS MUCH AS POSSIBLE IN EQ(II) */ X/* EQUATE K[II] TO UNREMOVABLE TERMS */ X X/* DEFINE PATTERN MATCHING RULES TO ISOLATE ARGS OF TRIG TERMS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE), XDEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2), XDEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2), X XFOR II:1 THRU NTRUNC DO( X X EQN[II]:EXPAND(TRIGREDUCE(EXPAND(EQ(II)))), X TEMP:EXPAND(EV(EQN[II],WGEN[II]=0)), X/* CHANGE SUM TO A LIST */ X TEMP1:ARGS(TEMP), X/* REMOVE CONSTANT TERMS */ X TEMP2:MAP(TRIGIDENTIFY,TEMP1), X/* ISOLATE ARGUMENTS OF TRIG TERMS */ X ARG1:APPLY1(TEMP2,COSARG,SINARG), X/* REMOVE DUPLICATE ARGUMENTS */ X ARG2:SETIFY(ARG1), X/* REMOVE RESONANT ARGUMENTS */ X ARG3:SUBLIST(ARG2,NOTRESP), X X/* CHOOSE WGEN TO ELIMINATE NONRESONANT TERMS */ X LENG:LENGTH(ARG3), X WGENTEMP1:0, X FOR JJ:1 THRU LENG DO( X WGENTEMP2:AAA*COS(PART(ARG3,JJ))+BBB*SIN(PART(ARG3,JJ)), X TEMP4:EV(EQN[II],WGEN[II]=WGENTEMP2,DIFF), X TEMP5:SOLVE([RATCOEF(TEMP4,COS(PART(ARG3,JJ))), X RATCOEF(TEMP4,SIN(PART(ARG3,JJ)))],[AAA,BBB]), X WGENTEMP1:WGENTEMP1+EV(WGENTEMP2,TEMP5)), X WGEN[II]:WGENTEMP1, X PRINT("WGEN[",II,"] ="), X PRINT(WGEN[II]), X K[II]:EXPAND(EV(EQN[II],DIFF)), X K[II]:EXPAND(RATSIMP(K[II])), X PRINT("THE TRANSFORMED HAMILTONIAN K[",II,"] ="), X PRINT(K[II])), X XKAMILTONIAN:SUM(K[II]*E^II,II,0,NTRUNC), XPRINT ("THE TRANSFORMED HAMILTONIAN IS "), XPRINT ("K =",KAMILTONIAN), X XCHOICE2:READ("DO YOU WANT TO SEE THE NEAR IDENTITY TRANSFORMATION (Y/N) ?"), XIF CHOICE2=N THEN GO(END), X X/* THE NEAR IDENTITY TRANSFORMATION */ XFOR II:1 THRU N DO( X JTRANS[II]:SUM(S(III,P[II])*E^III,III,0,NTRUNC), X PHITRANS[II]:SUM(S(III,Q[II])*E^III,III,0,NTRUNC)), XFOR II:1 THRU N DO( X PRINT(J[II],"=",JTRANS[II]), X PRINT(PHI[II],"=",PHITRANS[II])), X XEND, XKAMILTONIAN)$ X X XPOISSON(F,G):= X SUM(DIFF(F,Q[II])*DIFF(G,P[II])-DIFF(F,P[II])*DIFF(G,Q[II]),II,1,N)$ X XL(II,F):=POISSON(WGEN[II],F)$ X XS(II,F):=(IF II=0 THEN F ELSE SUM(L(II-M,S(M,F)),M,0,II-1)/II)$ X XEQ(II):=(H[II]+(DIFF(WGEN[II],T)+POISSON(WGEN[II],H[0]))/II X +SUM(L(II-M,K[M])+M*S(II-M,H[M]),M,1,II-1)/II)$ X XLZAP(ANY):=DIFF(ANY,T)+POISSON(ANY,H[0])$ X XTRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$ X XNOTRESP(EXP):=NOT IS(LZAP(EXP) = 0)$ X XSETIFY(LIST):=( X SET:[LIST[1]], X FOR I:2 THRU LENGTH(LIST) DO( X IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)), X SET)$ X X X X X X X X X X/**************************************************************************/ X/**************************************************************************/ X X X/*PROGRAM NUMBER 10: REDUCTION1(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN ONE DIFFERENTIAL EQUATION DEPENDING X ON ONE INDEPENDENT VARIABLE. SEE PAGE 170 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA". */ X X X X X/*THIS FILE CONTAINS REDUCTION1(), A FUNCTION TO PERFORM A LIAPUNOV- X SCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS OF NONLINEAR D.E.'S X OF THE FORM Y'' + F(Y,Y',ALPHA) = 0. Y = Y(X) IS DEFINED ON A REAL X INTERVAL WITH DIRICHLET OR NEUMANN BOUNDARY CONDITIONS AND F DEPENDS X ONLY LINEARLY ON ALPHA. X IT ALSO CONTAINS THESE ADDITIONAL FUNCTIONS: X SETUP() ALLOWS TO ENTER THE PROBLEM. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION. X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X PROJECT(EXP) ENSURES THAT =0. X NEUMANNBC(EXP) ACCOUNTS FOR NEUMANN BOUNDARY CONDITIONS. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. */ X X XREDUCTION1():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X XSETUP():=( X/*THE FUNCTION SETUP ASKS FOR THE VARIABLES OF THE D.E., THE BIFURCATION XPOINT, THE CRITICAL EIGENFUNCTION, THE X-INTERVAL, THE BOUNDARY XCONDITIONS AND THE DIFFERENTIAL EQUATION. */ XASSUME_POS:TRUE, XLS_LIST:[], XY:READ("ENTER DEPENDENT VARIABLE"), XPRINT("USE X AS THE INDEPENDENT VARIABLE AND ALPHA AS A PARAMETER TO VARY"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE ALPHA"), XPRINT("WE DEFINE LAM = ALPHA - ",CAL), XCFUN: READ("ENTER THE CRITICAL EIGENFUNCTION"), XDEPENDS([ZW,Y,W],[AMP,LAM,X]), XLEN:READ("WHAT IS THE LENGTH OF THE X-INTERVAL"), XNORM:INTEGRATE(CFUN^2,X,0,LEN), XPRINT("SPECIFY THE BOUNDARY CONDITIONS"), XPRINT("YOUR CHOICE FOR THE B.C. ON Y AT X=0 AND X=",LEN), XPRINT("ENTER 1 FOR Y=0, 2 FOR Y'=0"), XBCL:READ("B.C. AT 0?"), XBCR:READ("B.C. AT",LEN,"?"), XEQ:DIFF(Y,X,2) + X READ("THE D.E. IS OF THE FORM Y'' + F(Y,Y',ALPHA) = 0,ENTER F"), XEQLAM:EV(EQ,ALPHA=LAM+CAL), XPRINT(EQLAM), XDATABASE:[DIFF(W,AMP)=0,DIFF(W,LAM)=0], XSUB:Y=AMP*CFUN+W, XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), XNULLANS:READ(" ARE ZERO, Y/N") X)$ X X XG_POLY(I,J):=BLOCK( X/* THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) =0. IT REQUIRES KNOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE */ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, X Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*THIS DIFFERENTIAL OF W WILL BE ANNIHILATED BY THE PROJECTION ONTO THE XCRITICAL EIGENFUNCTION, HENCE WE SET IT TO ZERO ALREADY HERE */ X TEMP3:SUBST('DIFF(W,AMP,I,LAM,J)=0,TEMP2), X/* HERE WE INSERT ALL KNOWLEDGE GAINED THROUGH W_POLY */ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X/*PROJECTION ONTO CFUN, THE CRITICAL EIGENFUNCTION */ X TEMP5:TRIGREDUCE(CFUN*TEMP4), X BIFEQ[I,J]:RATSIMP(INTEGRATE(TEMP5,X,0,LEN)) X )$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). THE RESULT IS FED INTO DATABASE.*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X/*WE SUBSTITUTE ZW FOR THE UNKNOWN TAYLOR COEFFICIENT AND SOLVE XAN O.D.E. FOR ZW IN TEMP7 */ X TEMP3:SUBST(DIFF(W,AMP,I,LAM,J)=ZW,TEMP2), X/*NOW WE INSERT ALL KNOWLEDGE GAINED THROUGH PREVIOUS ITERATIONS*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,W=0), X TEMP5:TRIGREDUCE(TEMP4), X/*THIS ENSURES THAT THE R.H.S. FO THE D.E. TEMP6 IS ORTHOGONAL TO XTHE SOLUTION OF THE HOMOGENEOUS EQUATION */ X TEMP6:PROJECT(TEMP5), X TEMP7:ODE2(TEMP6,ZW,X), X/*SATISFY BOUNDARY CONDITIONS*/ X IF BCL*BCR=1 THEN TEMP8:BC2(TEMP7,X=0,ZW=0,X=LEN,ZW=0) X ELSE X TEMP8:NEUMANNBC(TEMP7), X/*MAKE SURE THAT = 0*/ X TEMP9:PROJECT(TEMP8), X ADDBASE:['DIFF(W,AMP,I,LAM,J)=RHS(TEMP9)], X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X XPROJECT(EXP):=( X PRO1:EV(EXP,ZW=0), X PRO2:INTEGRATE(PRO1*CFUN,X,0,LEN)/NORM, X EXPAND(EXP-PRO2*CFUN))$ XNEUMANNBC(EXP):=( X REXP:RHS(EXP), X NBC1:DIFF(REXP,X), X IF BCL=1 AND BCR=2 THEN NBC2:SOLVE([EV(REXP,X=0),EV(NBC1,X=LEN)], X [%K1,%K2]), X IF BCL=2 AND BCR=1 THEN NBC2:SOLVE([EV(REXP,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X IF BCL=2 AND BCR=2 THEN NBC2:SOLVE([EV(NBC1,X=LEN),EV(NBC1,X=0)], X [%K1,%K2]), X EV(EXP,NBC2))$ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X X X X X X X X/**********************************************************************/ X/**********************************************************************/ X X X X X/* PROGRAM NUMBER 11: REDUCTION2(), A LIAPUNOV-SCHMIDT REDUCTION FOR X STEADY STATE BIFURCATIONS IN SYSTEMS OF ORDINARY DIFFERENTIAL X EQUATIONS. SEE PAGE 176 AND PAGE 211 IN "PERTURBATION METHODS, X BIFURCATION THEORY AND COMPUTER ALGEBRA".*/ X X X X X XREDUCTION2():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X X X XSETUP():=(/*THE FUNCTION SETUP NEEDS THE DIMENSION N OF THE SYSTEM OF X O.D.E., THE N-DIMENSIONAL LIST OF VARIABLES Y, THE SYSTEM X OF O.D.E'S CALLED EQS, THE CRITICAL EIGENVECTOR CFUN AND X THE ADJOINT CRITICAL EIGENVECTOR AFUN.*/ X PRINTFLAG:TRUE, X LS_LIST:[], X N:READ("NUMBER OF EQUATIONS"), X Y:MAKELIST(READ("ENTER VARIABLE NUMBER",I),I,1,N), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE",ALPHA), X PRINT("WE DEFINE LAM =",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENVECTOR AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENVECTOR"), X KILL(W,ZWLIST), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W),[AMP,LAM]), X PRINT("ENTER THE DIFFERENTIAL EQUATION"), X EQS:MAKELIST(READ("DIFF(",Y[I],",T)="),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT(EQLAM), X WNULL:VFUN(W,0), X ZWNULL:VFUN(ZWLIST,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X DATABASE:APPEND(DATABASE,MAP(LAMBDA([U],DIFF(U,LAM)=0),W)), X NORM:AFUN.CFUN, X TEMP1:EV(EQLAM,SUB,DIFF), X PRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICENTS"), X NULLANS:READ(" ARE ZERO, Y/N") X )$ X X XG_POLY(I,J):=BLOCK(/*PROVIDED ALL NECESSARY KNOWLEDGE ABOUT THE TAYLOR X COEFFICIENTS OF W(AMP,LAM) IS STORED IN THE LIST X DATABASE, G_POLY CALCULATES ANY SPECIFIC TAYLOR X COEFFICIENT OF THE BIFURCATION EQUATION X G(AMP,LAM)= 0 VIA THE PROJECTION ONTO CFUN.*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS LS(",I,J,") IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X BIFEQ[I,J]:RATSIMP(TEMP4.AFUN))$ X X X XW_POLY(I,J):=BLOCK(/*THE FUNCTION W_POLY ALLOWS TO CALCULATE ITERATIVELY, X I.E. STARTING WITH THE LOWEST ORDER, ALL TAYLOR X COEFFICIENTS OF W(AMP,LAM) BY PROJECTING ONTO THE X COMPLEMENT OF CFUN*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W,AMP,",I,"LAM,",J," X IDENTICALLY ZERO, Y/N"), X IF ZEROANS = Y THEN X (ADDBASE:['DIFF(W,AMP,I,LAM,J)=0], X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X TEMP2:SUBST(WMAX_DIFF,TEMP2), X TEMP3:SUBST(DATABASE,TEMP2), X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL), X TEMP5:EV(TEMP4,ZWNULL).AFUN, X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:SOLVE(TEMP6,ZWLIST), X TEMP8:EV(ZWLIST,TEMP7), X TEMP9:RATSIMP(TEMP8.AFUN), X IF TEMP9 = 0 THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N),TEMP8) X ELSE X (TEMP10:RATSIMP(TEMP8 - TEMP9/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J), X U,1,N),TEMP10)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ X XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X X X X X X X X/**************************************************************************/ X/**************************************************************************/ X X X X/* PROGRAM NUMBER 12: A LIAPUNOV-SCHMIDT REDUCTION FOR STEADY STATE X BIFURCATIONS IN SYSTEMS OF PARTIAL DIFFERENTIAL EQUATIONS X DEPENDING ON ONE INDEPENDENT SPACE VARIABLE. SEE PAGE 189 IN X "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER ALGEBRA". X*/ X X X X X/*THIS FILE CONTAINS REDUCTION3(), A FUNCTION TO PERFORM A LIAPUNOV- XSCHMIDT REDUCTION FOR STEADY STATE BIFURCATIONS FROM N COUPLED XPARTIAL DIFFERENTIAL EQUATIONS ON ONE SPATIAL DIMENSION. XTHE FOLLOWING ADDITIONAL FUNCTIONS ARE SUPPLIED: X SETUP() ALLOWS THE PROBLEM TO BE ENTERED. X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM). X W_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN W(X;AMP,LAM). X SOLVE_ODE(EXP) SOLVES CERTAIN ORDINARY DIFFERENTIAL EQUATIONS VIA A FOURIER X MODE ANSATZ. X FEEDW(EXP) ENSURES THAT = 0 . X FIND_TRIG(EXP) IDENTIFIES FOURIER MODES. X SETIFY(LIST) TRANSFORMS A LIST INTO A SET. X G_EQ() ASSEMBLES THE BIFURCATION EQUATION. X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE, ...] X DIFFNULL(I,J) SETS THE DERIVATIVE DIFF(W,AMP,I,LAM,J) TO ZERO.*/ X XREDUCTION3():=BLOCK( X SETUP(), X ORDER:READ("TO WHICH ORDER DO YOU WANT TO CALCULATE"), X FOR I:2 THRU ORDER-1 DO W_POLY(I,0), X FOR I:1 THRU ORDER-2 DO W_POLY(I,1), X FOR I:1 THRU ORDER DO G_POLY(I,0), X FOR I:1 THRU ORDER-1 DO G_POLY(I,1), X G_EQ())$ X X XSETUP():=( /* THIS FUNCTION PERFORMS THE INPUT FOR THE LIAPUNOV-SCHMIDT X REDUCTION*/ XASSUME_POS:TRUE, XLS_LIST:[], XN:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), XY:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), XXVAR:READ("ENTER THE SPATIAL COORDINATE"), XALPHA:READ("ENTER THE BIFURCATION PARAMETER"), XCAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), XPRINT("WE DEFINE LAM = ",ALPHA - CAL), XCFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), XAFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS A LIST"), XKILL(W), XW:MAKELIST(CONCAT(W,I),I,1,N), XZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), XNULLIST:MAKELIST(0,I,1,N), XDEPENDS(APPEND(ZWLIST,W,Y),CONS(XVAR,[AMP,LAM])), XEQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), XEQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), XPRINT(EQLAM), XLEN:READ("WHAT IS THE LENGTH OF THE SPACE INTERVAL"), XWNULL:VFUN(W,0), XSUB:MAPLIST("=",Y,AMP*CFUN+W), XDATABASE:APPEND(DIFNULL(1,0),DIFNULL(0,1)), XZWNULL:VFUN(ZWLIST,0), XNORM:INTEGRATE(AFUN.CFUN,XVAR,0,LEN), XTEMP1:EV(EQLAM,SUB,DIFF), XPRINT("DO YOU KNOW APRIORI THAT SOME TAYLOR COEFFICIENTS ARE 0"), XNULLANS:READ("Y,N") X)$ X X XG_POLY(I,J):=BLOCK( X/*THIS IS A FUNCTION TO DETERMINE A PARTICULAR TAYLOR COEFFICIENT OF THE XBIFURCATION EQUATION G(AMP,LAM) = 0. IT REQUIRES KOWLEDGE ABOUT THE TAYLOR XCOEFFICIENTS OF W(AMP,LAM). THIS KNOWLEDGE IS STORED IN THE LIST DATABASE*/ X LS_LIST:CONS([K=I,L=J],LS_LIST), X IF NULLANS = Y THEN ( X ZEROANS:READ("IS G_POLY(",I,",",J,")IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN RETURN(BIFEQ[I,J]:0)), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*SET THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZERO*/ X TEMP3:SUBST(DIFNULL(I,J),TEMP2), X /*ENTER ALL INFORMATION IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE)), X /*PROJECT ONTO AFUN*/ X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ[I,J]:INTEGRATE(TRIGREDUCE(TEMP5),XVAR,0,LEN))$ X X X X XW_POLY(I,J):=BLOCK( X/* THIS FUNCTION ALLOWS TO DETERMINE ITERATIVELY ANY PARTICULAR TAYLOR XCOEFFICIENT OF THE FUNCTION W(AMP,LAM). IT RETURNS A DIFFERENTIAL EQUATION XFOR THE PARTICULAR COEFFICIENT OF INTEREST (CALLED ZW1,ZW2...). THIS D.E. XIS SOLVED VIA SOLVE_ODE AND THE RESULT IS FED INTO DATABASE FROM FEEDW*/ X IF NULLANS = Y THEN ( X ZEROANS:READ("IS DIFF(W(AMP,",I,",LAM,",J,") IDENTICALLY X ZERO, Y/N"), X IF ZEROANS = Y THEN ( X ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE))), X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X /*RENAME THE DERIVATIVES DIFF(W,AMP,I,LAM,J) TO ZW */ X TEMP3:SUBST(WMAX_DIFF,TEMP2), X /*ENTER ALL INFORMATION IN STORED IN DATABASE*/ X D_BASE_LENGTH:LENGTH(DATABASE), X FOR II THRU D_BASE_LENGTH DO X TEMP3:EV(SUBST(DATABASE[D_BASE_LENGTH+1-II],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EV(TEMP3,AMP=0,LAM=0,WNULL,INTEGRATE), X /*THIS IS THE PROJECTION Q ONTO THE RANGE OF THE X LINEAR DIFFERENTIAL OPERATOR IN THE PROBLEM*/ X TEMP5:INTEGRATE(EV(TEMP4,ZWNULL).AFUN,XVAR,0,LEN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X TEMP7:TRIGREDUCE(TEMP6), X /*THE SET OF O.D.E.'S TO SOLVE */ X W_DE:EXPAND(TEMP7), X TEMP8:EV(W_DE,VFUN(ZWLIST,0)), X /*IF THE PARTICULAR SOLUTION OF W_DE IS ZERO THEN W=0 */ X IF NULLIST=TEMP8 THEN ( ADDBASE:DIFNULL(I,J), X DATABASE:APPEND(DATABASE,ADDBASE), X RETURN(ADDBASE)), X TEMP9:SOLVE_ODE(TEMP8), X FEEDW(TEMP9) X )$ X X XSOLVE_ODE(EXP):=( X/*THIS FUNCTION SOLVE THE D.E. W_DE BY A FOURIER MODE ANSATZ*/ X TRIGFUN:[], X CONST:FALSE, X FOR I THRU N DO ( X /*DETERMINE THE FOURIER MODES*/ X TRIG1:EXP[I], X IF TRIG1 # 0 THEN ( X TRIG2:APPLY1(TRIG1,SINNULL,COSNULL), X IF TRIG2 # 0 THEN ( X CONST:TRUE, X TRIG1:TRIG1-TRIG2), X TRIGFUN:APPEND(FIND_TRIG(TRIG1),TRIGFUN))), X SOL1:DELETE(DUMMY,SETIFY(TRIGFUN)), X /*MAKE AN ANSATZ*/ X ANSATZ:MAKELIST(SUM(AM[I,J]*SOL1[I],I,1,LENGTH(SOL1)),J,1,N), X SOL2:EV(W_DE,MAP("=",ZWLIST,ANSATZ),DIFF), X SOL3:MAKELIST(RATCOEF(SOL2,I),I,SOL1), X EQLIST:[], X FOR I THRU LENGTH(SOL3) DO EQLIST:APPEND(EQLIST,SOL3[I]), X VARLIST:[], X FOR I THRU N DO X FOR J THRU LENGTH(SOL1) DO VARLIST:CONS(AM[J,I],VARLIST), X /*FIND THE AMPLITUDES OF THE FOURIER MODES*/ X SOL4:SOLVE(EQLIST,VARLIST), X /*SOLVE FOR THE CONSTANT FOURIER MODE IF NECESSARY*/ X CANSATZ:0, X IF CONST = TRUE THEN (CANSATZ:MAKELIST(CONCAT(CON,I),I,1,N), X CSOL1:EV(W_DE,MAP("=",ZWLIST,CANSATZ),DIFF), X CSOL2:APPLY1(CSOL1,SINNULL,COSNULL), X CSOL3:SOLVE(CSOL2,CANSATZ)), X EV(ANSATZ+CANSATZ,SOL4,CSOL3))$ X X XFEEDW(EXP):=( X/* THIS FUNCTION ALLOWS THE RESULT OF THE ODE-SOLVER TO BE ENTERED INTO X THE LIST DATABASE. IT CHECKS FOR ORTHOGONALITY TO THE CRITICAL ADJOINT X EIGENFUNCTION AND REMOVES SOLUTIONS OF THE HOMOGENEOUS EQUATION (I.E. X NONORTHOGONAL TERMS)*/ X F2:INTEGRATE(EXP.AFUN,XVAR,0,LEN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X EXP) X ELSE (ORTHO_RESULT:RATSIMP(EXP- F2/NORM*CFUN), X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X PRINT(ADDBASE))$ X/*COLLECT ALL INFORMATION STORED IN BIFEQ AND ASSEMBLE THE BIFURCATION XEQUATION*/ XG_EQ():= SUM(EV(AMP^K*LAM^L/K!*BIFEQ[K,L],LS_LIST[N]),N,1,LENGTH(LS_LIST))$ X X XSETIFY(LIST):=(/*TRANSFORMS A LIST INTO A SET, I.E. REMOVES DUPLICATES*/ X SET:[LIST[1]], X FOR I :2 THRU LENGTH(LIST) DO X IF NOT MEMBER(LIST[I],SET) THEN X SET:CONS(LIST[I],SET), X SET)$ X XFIND_TRIG(EXP):=(/*FINDS THE FOURIER MODES*/ X F_A1:ARGS(EXP+DUMMY), X F_A2:APPLY1(F_A1,SINFIND,COSFIND) X )$ X X/*AUXILIARY FUNCTIONS */ XMATCHDECLARE([DUMMY1,DUMMY2],TRUE)$ XDEFRULE(COSFIND,DUMMY1*COS(DUMMY2),COS(DUMMY2))$ XDEFRULE(SINFIND,DUMMY1*SIN(DUMMY2),SIN(DUMMY2))$ XDEFRULE(COSNULL,COS(DUMMY1),0)$ XDEFRULE(SINNULL,SIN(DUMMY1),0)$ X X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ XDIFNULL(I,J):=MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J)=0),W) $ X X X X X X X X X/*************************************************************************/ X/*************************************************************************/ X X X X X/*PROGRAM NUMBER 13: CONTAINS THE BUILDING BLOCKS TO PERFORM STEADY X STATE BIFURCATIONS FOR MORE PARTIAL DIFFERENTIAL EQUATIONS. SEE PAGES X 198-202 IN "PERTURBATION METHODS, BIFURCATION THEORY AND COMPUTER X ALGEBRA" */ X X X X X/*THE FOLLOWING FUNCTIONS CAN BE USED TO PERFORM A LIAPUNOV-SCHMIDT REDUCTION X FOR THE 2D BENARD PROBLEM */ X X/*SETUP() ALLOWS THE PROBLEM TO BE ENTERED, X G_POLY(I,J) CALCULATES THE COEFFICIENT OF AMP^I LAM^J IN THE BIFURCATION X EQUATION G(AMP,LAM), X W_POLY(I,J) DETERMINES A DIFFERENTIAL EQUATION, WHOSE SOLUTION IS THE X COEFFICIENT OF AMP^I LAM^J IN W(AMP,LAM), X FEEDW() ALLOWS THIS SOLUTION TO BE ENTERED INTO THE LIST DATABASE, X INT(EXP) FINDS MULTIPLE DEFINITE INTEGRALS EFFECTIVELY, X VFUN(LIST,VALUE) CREATES THE SUBSTITUTION LIST X [LIST[1] = VALUE, LIST[2] = VALUE ...] X*/ X XSETUP():=( X N:READ("ENTER THE NUMBER OF DIFFERENTIAL EQUATIONS"), X Y:READ("ENTER THE DEPENDENT VARIABLES AS A LIST"), X SPACE:READ("ENTER NUMBER OF SPATIAL COORDINATES"), X XVAR:READ("ENTER THE SPATIAL COORDINATES AS A LIST"), X ALPHA:READ("ENTER THE BIFURCATION PARAMETER"), X CAL:READ("ENTER THE CRITICAL BIFURCATION VALUE"), X PRINT("WE DEFINE LAM = ",ALPHA - CAL), X CFUN:READ("ENTER THE CRITICAL EIGENFUNCTION AS A LIST"), X AFUN:READ("ENTER THE ADJOINT CRITICAL EIGENFUNCTION AS ALIST"), X KILL(W), X W:MAKELIST(CONCAT(W,I),I,1,N), X ZWLIST:MAKELIST(CONCAT(ZW,I),I,1,N), X DEPENDS(APPEND(ZWLIST,W,Y),APPEND([AMP,LAM],XVAR)), X EQS:MAKELIST(READ("ENTER THE DIFFERENTIAL EQUATION NUMBER",I),I,1,N), X EQLAM:EV(EQS,EV(ALPHA) = LAM + CAL), X PRINT("ENTER THE LOWER LEFT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X LBOUND:READ("[",X[1],"=...,]"), X PRINT("ENTER THE UPPER RIGHT CORNER OF THE",N,"DIMENSIONAL SPACE X INTERVAL"), X UBOUND:READ("[",X[1],"=...,]"), X WNULL:VFUN(W,0), X SUB:MAPLIST("=",Y,AMP*CFUN+W), X DATABASE:MAP(LAMBDA([U],DIFF(U,AMP)=0),W), X ZWNULL:VFUN(ZWLIST,0), X NORM:INT(AFUN.CFUN), X TEMP1:EV(EQLAM,SUB,DIFF), X EQLAM X )$ X X X XG_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = 0),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:RATSIMP(TEMP4.AFUN), X BIFEQ:INT(TEMP5))$ X X X X XW_POLY(I,J):=BLOCK( X WMAX_DIFF:MAP(LAMBDA([U],'DIFF(U,AMP,I,LAM,J) = CONCAT(Z,U)),W), X TEMP2:DIFF(TEMP1,AMP,I,LAM,J), X DERIVSUBST:TRUE, X TEMP3:SUBST(WMAX_DIFF,TEMP2), X M:LENGTH(DATABASE), X FOR I THRU M DO TEMP3:EV(SUBST(DATABASE[M+1-I],TEMP3),DIFF), X DERIVSUBST:FALSE, X TEMP4:EXPAND(EV(TEMP3,AMP=0,LAM=0,WNULL)), X TEMP5:INT(EV(TEMP4,ZWNULL).AFUN), X TEMP6:TEMP4-TEMP5/NORM*CFUN, X FOR I THRU SPACE DO TEMP7:TRIGREDUCE(TEMP6,XVAR[I]), X W_DE:RATSIMP(TEMP7 ), X PRINT("NOW SOLVE THE EQUATIONS",W_DE,"=0! THEY ARE GIVEN IN W_DE"), X PRINT("CALL FEEDW() TO PROCEED"))$ X X X XFEEDW():=( X ADDBASE:[], X RESULT:MAKELIST((PRINT("ENTER",ZWLIST[K]),READ()),K,1,N), X F2:INT(RESULT.AFUN), X IF RATSIMP(F2)=0 X THEN X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X RESULT) X ELSE (ORTHO_RESULT:RATSIMP(RESULT- F2/NORM*CFUN), X/*PROJECTION Q*/ X ADDBASE:MAP("=",MAKELIST('DIFF(W[U],AMP,I,LAM,J),U,1,N), X ORTHO_RESULT)), X DATABASE:APPEND(DATABASE,ADDBASE), X ADDBASE)$ X X XINT(EXP):=(%INT:EXP,FOR I THRU SPACE DO %INT:INTEGRATE(TRIGREDUCE(%INT, X XVAR[I]),XVAR[I]), X RATSIMP(EV(%INT,UBOUND)-EV(%INT,LBOUND)))$ X XVFUN(LIST,VALUE):=MAP(LAMBDA([U],U=VALUE),LIST)$ END_OF_FILE if test 39932 -ne `wc -c <'program16'`; then echo shar: \"'program16'\" unpacked with wrong size! fi # end of 'program16' fi echo shar: End of shell archive. exit 0