/* PROGRAM NUMBER 6: TWOVAR(), A FUNCTION TO PERFORM A TWO VARIABLE
   EXPANSION TO ORDER EPSILON FOR N WEAKLY COUPLED PERTURBED NONLINEAR
   OSCILLATORS. SEE PAGE 94 IN "PERTURBATION METHODS, BIFURCATION
   THEORY AND COMPUTER ALGEBRA". */



TWOVAR():=BLOCK(

CHOICE:READ("DO YOU WANT TO ENTER NEW DATA (Y/N)"),

IF CHOICE = N THEN GO(JUMP1),

/* INPUT D.E.'S */
N:READ("NUMBER OF D.E.'S"),

PRINT("THE",N,"D.E.'S WILL BE IN THE FORM:"),
PRINT("X[I]'' + W[I]^2 X[I] = E F[I](X[1],...,X[",N,"],T)"),

FOR I:1 THRU N DO
    X[I]:READ("ENTER SYMBOL FOR X[",I,"]"),

FOR I:1 THRU N DO
    DEPENDS(X[I],T),

FOR I:1 THRU N DO
    W[I]:READ("ENTER W[",I,"]"),

FOR I:1 THRU N DO
    F[I]:READ("ENTER F[",I,"]"),

JUMP1,

/* UPDATE EQS FOR SUBSTITUTION OF RESONANT VALUES ON 2ND TIME THRU */
FOR I:1 THRU N DO(
   W[I]:EV(W[I]),
   F[I]:EV(F[I])),

/* ECHO BACK THE D.E.'S */
PRINT("THE D.E.'S ARE ENTERED AS:"),
FOR I:1 THRU N DO
    PRINT(X[I],"'' +",W[I]^2*X[I],"=",E*F[I]),

PRINT("THE METHOD ASSUMES A SOLUTION IN THE FORM:"),
PRINT("X[I] = X0[I] + E X1[I]"),
PRINT("WHERE X0[I] = A[I](ETA) COS W[I] XI + B[I](ETA) SIN W[I] XI"),
PRINT("WHERE XI = T AND ETA = E T"),

/* DON'T USE A OR B AS PARAMETERS IN THE GIVEN D.E.'S */
DEPENDS([A,B],ETA),

FOR I:1 THRU N DO
    X0[I]:A[I]*COS(W[I]*XI)+B[I]*SIN(W[I]*XI),

FOR I:1 THRU N DO
    G[I]:EV(F[I],T=XI),

FOR I:1 THRU N DO(
    FOR J:1 THRU N DO
        G[I]:EV(G[I],X[J]::X0[J])),
    
FOR I:1 THRU N DO(
    G[I]:G[I]-2*DIFF(X0[I],XI,1,ETA,1),
    G[I]:EV(G[I],DIFF),
    G[I]:EXPAND(TRIGREDUCE(EXPAND(G[I])))),    

/* COLLECT SECULAR TERMS */
FOR I:1 THRU N DO(
    S[I]:COEFF(G[I],SIN(W[I]*XI)),
    C[I]:COEFF(G[I],COS(W[I]*XI))),

/* DISPLAY SECULAR TERMS */
PRINT("REMOVAL OF SECULAR TERMS IN THE X1[I] EQS. GIVES:"),

FOR I:1 THRU N DO(
    PRINT(S[I],"= 0"),
    PRINT(C[I],"= 0")),

ABEQS:APPEND(MAKELIST(S[I],I,1,N),MAKELIST(C[I],I,1,N)),

CHOICE2:READ("DO YOU WANT TO TRANSFORM TO POLAR COORDINATES (Y/N)"),

IF CHOICE2=N THEN GO(JUMP2),

/* TRANSFORM TO POLAR COORDINATES */
DEPENDS([R,THETA],ETA),
TRANS:MAKELIST([A[I]=R[I]*COS(THETA[I]),B[I]=R[I]*SIN(THETA[I])],I,1,N),
INTEQS:EV(ABEQS,TRANS,DIFF),
RTHEQS:SOLVE(INTEQS,APPEND(MAKELIST(DIFF(R[I],ETA),I,1,N),
                     MAKELIST(DIFF(THETA[I],ETA),I,1,N))),
RTHEQS:TRIGSIMP(RTHEQS),
RTHEQS:EXPAND(TRIGREDUCE(RTHEQS)),
PRINT(RTHEQS),

JUMP2,

CHOICE3:READ("DO YOU WANT TO SEARCH FOR RESONANT PARAMETER VALUES (Y/N)"),

IF CHOICE3=N THEN GO(END),

/* OBTAIN FREQUENCIES WHICH APPEAR ON RHS'S OF D.E.'S */

/* DEFINE PATTERN MATCHING RULES TO ISOLATE FREQS */
MATCHDECLARE([DUMMY1,DUMMY2],TRUE),
DEFRULE(COSARG,DUMMY1*COS(DUMMY2),DUMMY2),
DEFRULE(SINARG,DUMMY1*SIN(DUMMY2),DUMMY2),

FOR I:1 THRU N DO(
/* CHANGE SUM TO A LIST */
    TEMP1:ARGS(G[I]),
/* REMOVE CONSTANT TERMS */
    TEMP2:MAP(TRIGIDENTIFY,TEMP1),
/* ISOLATE ARGUMENTS OF TRIG TERMS */
    TEMP3:APPLY1(TEMP2,COSARG,SINARG),
    TEMP4:EV(TEMP3,XI=1),
/* REMOVE DUPLICATE ARGUMENTS */
    FREQ[I]:SETIFY(TEMP4)),

/* DISPLAY FREQUENCIES */                   
FOR I:1 THRU N DO(
    PRINT(X[I],"EQ'S RESONANT FREQ =",W[I]),
    PRINT("FREQS ON RHS =",FREQ[I])),

JUMP3,

PAR:READ("WHICH PARAMETER TO SEARCH FOR ?"),

/* COMPUTE RESONANT VALUES OF DESIRED PARAMETER */
RESVALS:[],
FOR I:1 THRU N DO(
   FOR J:1 THRU LENGTH(FREQ[I]) DO(
      RES:SOLVE(PART(FREQ[I],J)=W[I],PAR),
      IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES),
      RES:SOLVE(PART(FREQ[I],J)=-W[I],PAR),
      IF (RES#[] AND RES#ALL) THEN RESVALS:APPEND(RESVALS,RES))),

/* ELIMINATE DUPLICATE PARAMETER VALUES */
RESVALUES:SETIFY(RESVALS),

/* DISPLAY RESONANT PARAMETER VALUES */
PRINT(RESVALUES),

CHOICE4:READ("DO YOU WANT TO SEARCH FOR ANOTHER PARAMETER (Y/N) ?"),

IF CHOICE4=Y THEN GO(JUMP3),

END," ")$

/* AUXILIARY FUNCTIONS */

TRIGIDENTIFY(EXP):=IF FREEOF(SIN,EXP) AND FREEOF(COS,EXP) THEN 0 ELSE EXP$

SETIFY(LIST):=(
    SET:[LIST[1]],
    FOR I:2 THRU LENGTH(LIST) DO(
        IF NOT MEMBER(LIST[I],SET) THEN SET:CONS(LIST[I],SET)),
    SET)$