eval_when([translate,batch,demo,load,loadfile], dv(a)::=buildq([a],define_variable(a,'a,any)) dvl(a)::=buildq([a],define_variable(a,[],list)),TTYOFF:TRUE)$ /* (c) Copyright 1981 Massachusetts Institute of Technology */ /* DYNAMALLOC:TRUE$ */ define_variable(IEQNPRINT,TRUE,boolean)$ define_variable(TECHLIST,'[FIRST,ALL,VLFRNK,TRANSFORM,COLLOCATE, FLFRNK2ND,TAILOR,FREDSERIES,NEUMANN, FLFRNK1ST,ABEL,FIRSTKINDSERIES],list)$ dv(tech)$ dv(px)$ dv(eqn)$ dv(fx)$ dv(uvar)$ dv(xvar)$ dv(kxu)$ dv(ax)$ dv(bx)$ dv(iesoln)$ dv(napprox)$ dv(pxlist)$ dv(eqnlist)$ dv(inteqn)$ dv(lowlim)$ dv(pofx)$ dv(iclist)$ dvl(unklist)$ dv(highlim)$ dv(%c)$ dv(eqno)$ dvl(ueqnlist)$ dv(x)$ dv(xhat)$ dv(opr)$ dv(guesslist)$ dv(guess)$ IEQN([ARGLIST]):= IF LENGTH(ARGLIST)<2 THEN ERROR("IEQN requires at least two arguments") ELSE IEQN1(ARGLIST[1], ARGLIST[2], IF LENGTH(ARGLIST)>2 THEN ARGLIST[3] ELSE [], IF LENGTH(ARGLIST)>3 THEN ARGLIST[4] ELSE [], IF LENGTH(ARGLIST)>4 THEN ARGLIST[5] ELSE [])$ OPR(E):=IF MAPATOM(E) THEN 'NONE ELSE INPART(E,0)$ ARG(E):=INPART(E,1)$ ALIAS(EXITFOR, RETURN)$ /* IDENTIFIES FUNCTIONS UNKNOWN TO MACSYMA */ UNFUN(E):=IS(?GETCHAR(E,1)='?\$)$ /* ***** MAIN PROGRAM *****/ IEQN1(EQNLIST,PXLIST,TECH,NAPPROX,GUESSLIST):= BLOCK([IESOLN, XVAR, UVAR, AX, BX, KXU, FX, EQN, PX, GUESS, SYSTEM, INFLAG, DISPFLAG, SINGSOLVE, SOLVERADCAN, KEEPFLOAT], INFLAG:TRUE, DISPFLAG:FALSE, SYSTEM:FALSE, IF NOT LISTP(EQNLIST) THEN EQNLIST:[EQNLIST], IF NOT LISTP(PXLIST) THEN PXLIST:[PXLIST], IF LENGTH(EQNLIST)#LENGTH(PXLIST) THEN ERROR("NUMBER OF EQUATIONS # NUMBER OF UNKNOWNS"), FOR UNK IN PXLIST DO IF MAPATOM(UNK) OR LENGTH(UNK)#1 OR NOT ATOM(XVAR:ARG(UNK)) OR NUMBERP(XVAR) THEN ERROR(UNK," IMPROPERLY SPECIFIED UNKNOWN"), IF TECH=[] THEN MYPRINT("DEFAULT 3RD ARG, TECHNIQUE: ",TECH:'FIRST) ELSE IF NOT MEMBER(TECH,TECHLIST) THEN ERROR(TECH," INVALID TECHNIQUE - SEE TECHLIST"), IF NAPPROX=[] THEN MYPRINT("DEFAULT 4TH ARG, NUMBER OF ITERATIONS OR COLL. PARMS.: ",NAPPROX:1) ELSE IF NOT INTEGERP(NAPPROX) OR NAPPROX<=0 THEN ERROR("NAPPROX MUST BE A POS. INTEGER"), IF GUESSLIST=[] THEN MYPRINT("DEFAULT 5TH ARG, INITIAL GUESS: ",GUESSLIST:'NONE), IF GUESSLIST='NONE THEN GUESS:'NONE ELSE (IF NOT LISTP(GUESSLIST) THEN GUESSLIST:[GUESS:GUESSLIST], IF LENGTH(GUESSLIST)#LENGTH(PXLIST) THEN ERROR("NUMBER OF GUESSES # NUMBER OF UNKNOWNS")), SINGSOLVE:SOLVERADCAN:KEEPFLOAT:TRUE, IESOLN:[], EQNLIST:MAPLIST(LAMBDA([EQN],NUM(RATSIMP(LHS(EQN)-RHS(EQN)))),EQNLIST), IF LENGTH(PXLIST)=1 THEN (PX:FIRST(PXLIST), EQN:FIRST(EQNLIST)) ELSE SYSTEM:TRUE, IF NOT SYSTEM THEN (TRY('VLFRNK), IF TECH='DONE THEN RETURN(REVERSE(IESOLN))), TRY('TRANSFORM), IF TECH='DONE THEN RETURN(REVERSE(IESOLN)), FOR EQNLIST IN SOLVESYS(EQNLIST,PXLIST) DO (TRY('FLFRNK2ND), TRY('TAILOR), IF NOT SYSTEM THEN TRY('FREDSERIES), TRY('NEUMANN)), IF NOT SYSTEM AND FIRSTKIND() THEN (TRY('FLFRNK1ST), TRY('ABEL), TRY('FIRSTKINDSERIES)), TRY('COLLOCATE), RETURN(REVERSE(IESOLN)))$ /* INVOKES SOLUTION METHOD CATCHING AND CONTAINING ERRORS. GLOBALS: TECH */ TRY(ROUTINE):=IF MEMBER(TECH,['ALL, 'FIRST, ROUTINE]) THEN BLOCK([ATTEMPT], ATTEMPT:ERRCATCH(CATCH(APPLY(ROUTINE,[]))), IF ATTEMPT=[] THEN (?PRINC('? IN\ ), ?PRINC(?STRIPDOLLAR(ROUTINE)), ?TERPRI()) ELSE IF FIRST(ATTEMPT) AND TECH='FIRST THEN TECH:'DONE)$ /* TEST FOR 1ST KIND AND EXTRACTS PARTS. GLOBALS: EQN, KXU, UVAR, XVAR, AX, BX */ FIRSTKIND():=BLOCK([INTEGRAL, XPOINTS, CONSISTENT], IF NOT FREEOF(PX,EQN) THEN RETURN(FALSE), INTEGRAL:CATCH(FINDINT(EQN)), IF INTEGRAL=FALSE THEN RETURN(FALSE), FX:MYSOLVE(EQN,INTEGRAL), IF FX=[] THEN RETURN(FALSE), FX:FIRST(FX), UVAR:INPART(INTEGRAL,2), KXU:LIN(ARG(INTEGRAL), SUBST(UVAR,XVAR,PX)), IF KXU=FALSE THEN RETURN(FALSE), AX:INPART(INTEGRAL,3), BX:INPART(INTEGRAL,4), IF KXU[2]#0 THEN FX:RATSIMP(FX-MYINT(KXU[2],UVAR,AX,BX)), XPOINTS:MYSOLVE(AX-BX,XVAR), KXU:KXU[1], IF XPOINTS=[] THEN RETURN(TRUE) ELSE CONSISTENT:FALSE, FOR XV IN XPOINTS DO IF RATSUBST(XV,XVAR,FX)=0 THEN EXITFOR(CONSISTENT:TRUE), RETURN(CONSISTENT))$ /* RETURNS FIRST INTEGRAL FOUND IN EXPR */ FINDINT(EXPR):=IF NOT MAPATOM(EXPR) THEN (IF OPR(EXPR)#'?%INTEGRATE THEN (FOR I IN EXPR DO FINDINT(I), FALSE) ELSE THROW(EXPR))$ /* GLOBALS: IEQNPRINT, IESOLN */ PRINTSOLN(SOLN,TECH,ORDER,KIND):=BLOCK([DISPFLAG, DEN], DISPFLAG:IEQNPRINT, IF KIND#[] THEN KIND:[KIND], IF ORDER#[] THEN KIND:CONS(ORDER,KIND), IF NOT LISTP(SOLN) THEN SOLN:[SOLN], SOLN:MAPLIST(LAMBDA([ELMT], block([elem,den], ELeM:RATSIMP(ELMT), DEN:DENOM(ELeM), IF NUMBERP(DEN) AND DEN#1 AND OPR(NUM(ELeM))="+" THEN MULTTHRU(ELeM) ELSE ELeM)),SOLN), IF REST(SOLN)=[] THEN SOLN:FIRST(SOLN), IESOLN:CONS(?DISPLINE(CONS(SOLN,CONS(TECH,KIND))),IESOLN), RETURN(TRUE))$ ALIAS(MYINT, INTEGRATE)$ MYPRINT(MSG1,MSG2):=IF IEQNPRINT THEN PRINT(MSG1,MSG2)$ /* IF EXPR = A*VAR+B RETURNS [A,B] ELSE FALSE */ LIN(EXPR,VAR):= BLOCK([LC], LC:RATSIMP(DIFF(EXPR,VAR)), IF FREEOF(VAR,LC) THEN [LC, RATSUBST(0,VAR,EXPR)])$ /* SOLVES EQN FOR UNK AND RETURNS A LIST OF ACTUAL SOLUTION VALUES */ MYSOLVE(EQN,UNK):=BLOCK([RESULT], RESULT:[], EQN:SOLVE(EQN,UNK), EQN:MAPLIST('RHS, apply('EV,[EQN])), FOR ANS IN EQN DO IF FREEOF(UNK, ANS) THEN RESULT:CONS(ANS, RESULT), RETURN(RESULT))$ /* SIMILAR TO MYSOLVE BUT FOR SYSTEMS. */ SOLVESYS(EQNS,UNKS):= (EQNS:SOLVE(EQNS,UNKS), EQNS:apply('EV,[EQNS]), MAPLIST(LAMBDA([EL], IF LISTP(EL) THEN MAPLIST('RHS,EL) ELSE [RHS(EL)]), EQNS))$ /* IF EXPR CAN BE EXPRESSED AS SUM(F[I](X)*G[I](U),I,1,N) THEN RESULT IS [[F1,G1],...[FN,GN]] ELSE []. GLOBALS: XVAR */ SUMFACTORS(EXPR,UVAR):=BLOCK([OTHER, REM], EXPR:CATCH(FRNK(EXPR)), IF EXPR=FALSE THEN RETURN([]), IF FREEOF(UVAR,EXPR) THEN RETURN([[EXPR,1]]), IF FREEOF(XVAR,EXPR) THEN RETURN([[1,EXPR]]), EXPR:EXPAND(EXPR), OTHER:MYPARTITION(EXPR,UVAR), EXPR:MYPARTITION(EXPR,XVAR), IF OTHER[1]0 THEN RETURN(FRNK(AR)^UP) ELSE THROW(FALSE), UP:PLUSSPLIT(EXPAND(UP)), IF UP#[] AND FREEOF(UVAR,XVAR,AR) THEN RETURN(BOX(AR^UP[1])*BOX(AR^UP[2])) ELSE THROW(FALSE)), IF MEMBER(OPR(AR),["*","+"]) THEN (E:PARTITION(AR,XVAR), IF NOT FREEOF(UVAR,E[2]) THEN THROW(FALSE)), IF OP='?%LOG AND OPR(AR)="*" THEN LOG(E[1])+LOG(E[2]) ELSE IF OPR(AR)="+" THEN IF OP='?%SIN THEN SIN(E[1])*COS(E[2])+COS(E[1])*SIN(E[2]) ELSE IF OP='?%COS THEN COS(E[1])*COS(E[2])-SIN(E[1])*SIN(E[2]) ELSE IF OP='?%SINH THEN SINH(E[1])*COSH(E[2])+COSH(E[1])*SINH(E[2]) ELSE IF OP='?%COSH THEN COSH(E[1])*COSH(E[2])+SINH(E[1])*SINH(E[2]) ELSE THROW(FALSE) ELSE THROW(FALSE))$ /* IF E=F(X)+G(U) RETURNS [F(X),G(U)] ELSE []. GLOBALS: XVAR, UVAR */ PLUSSPLIT(E):= IF OPR(E)="+" THEN (E:PARTITION(E,UVAR), IF FREEOF(XVAR,E[2]) THEN E ELSE []) ELSE IF FREEOF(UVAR,E) THEN [E,0] ELSE IF FREEOF(XVAR,E) THEN [0,E] ELSE []$ /* SOLVES EXS FOR UNKS AND SUBSTITUTES IN FORM */ SOLVEANDSUBST(EXS,UNKS,FORM,TECH):=BLOCK([SOLN,TRIAL], IF (SOLN:SOLVE(EXS,UNKS))=[] THEN (PRINT("FOR A ",TECH," SOLUTION SUBSTITUTE IN THE EXPRESSION:"), LDISP(FORM), APPLY('PRINT,CONS("THE VALUES OF ",UNKS)), PRINT("THAT MAKE THE FOLLOWING EXPRESSIONS SIMULTANEOUSLY ZERO"), APPLY('LDISP,EXS), RETURN(FALSE)), FOR SOL IN apply('EV,[SOLN]) DO (TRIAL:apply('EV,[FORM,SOL]), APPLY('PRINTSOLN,APPEND([TRIAL,TECH], IF TECH='COLLOCATE THEN [NAPPROX,TESTSOLN(TRIAL)] ELSE [[],[]]))), RETURN(TRUE))$ /* TESTS SOLUTION FOR EXACTNESS. GLOBALS: EQNLIST, PXLIST */ TESTSOLN(RESULTLIST):=BLOCK([FLAG], APPLY('LOCAL,MAPLIST('OPR,PXLIST)), FLAG:[], MAPLIST('DEFINE,PXLIST,RESULTLIST), RESULTLIST:apply('EV,[EQNLIST,'INTEGRATE,'RATSIMP]), FOR VAL IN RESULTLIST DO IF VAL#0 THEN EXITFOR(FLAG:'APPROXIMATE), RETURN(FLAG))$ /* FOUT=H(X,U)+Q(X)+R(U). RETURNS [N,H(X,U),Q(X),R(U)] WHERE N IS THE RANK OF FOUT. GLOBALS: XVAR, UVAR */ MYPARTITION(FOUT,VAR):=BLOCK([QX, RU, REM, CON], IF OPR(FOUT)#"+" OR OPR(FOUT:FACTOROUT(FOUT,VAR))#"+" THEN RETURN([1,[FOUT],0,0]), REM:CON:QX:RU:0, FOR FC IN FOUT DO IF FREEOF(UVAR,FC) THEN (IF FREEOF(XVAR,FC) THEN CON:CON+FC ELSE QX:QX+FC) ELSE IF FREEOF(XVAR,FC) THEN RU:RU+FC ELSE REM:REM+FC, FOUT:IF REM=0 THEN 0 ELSE IF OPR(REM)="+" THEN LENGTH(REM) ELSE (REM:[REM], 1), IF QX#0 THEN (FOUT:FOUT+1, QX:QX+CON, CON:0), IF RU#0 THEN (FOUT:FOUT+1, RU:RU+CON), RETURN([FOUT, REM, QX, RU]))$ /* REPLACES ALL INTEGRALS AND UNKNOWN FUNCTIONS IN EX WITH ZERO */ ZEROINT(EX):= IF MAPATOM(EX) THEN EX ELSE IF OPR(EX)='?%INTEGRATE OR UNFUN(OPR(EX)) THEN 0 ELSE MAP('ZEROINT,EX)$ /* ----------- THE FOLLOWING FUNCTIONS APPLY TO 1ST OR 2ND KIND EQNS */ /* VARIABLE-LIMIT FINITE-RANK 1ST OR 2ND KIND. REDUCTION TO ODE. GLOBALS: EQN, XVAR, PX */ VLFRNK():=BLOCK([INITIAL, LOWLIM, UNKLIST, INTEQN, POFX, KIND, FIRSTKIND, DIFFLIST, ICLIST], POFX:PX, INITIAL:TRUE, UNKLIST:DIFFLIST:[], KIND:'INCOMPLETE, FIRSTKIND:IS(FREEOF(PX,EQN)), INTEQN:VCONVERT(EQN), ICLIST:[XVAR=LOWLIM], FOR COUNT THRU LENGTH(UNKLIST) DO ( IF FIRSTKIND THEN FIRSTKIND:FALSE ELSE (IF INITIAL THEN GETICS(), POFX:DIFF(POFX,XVAR)), DIFFLIST:CONS(INTEQN,DIFFLIST), INTEQN:DIFF(INTEQN,XVAR), IF FREEOF('?%INTEGRATE, INTEQN) THEN EXITFOR(FALSE)), APPLY('REMOVE,[OPR(PX), 'ATVALUE]), IF NOT FREEOF('?%INTEGRATE, INTEQN) THEN ( DIFFLIST:SOLVE(DIFFLIST, UNKLIST), IF DIFFLIST=[] THEN THROW(FALSE), INTEQN:apply('EV,[INTEQN,DIFFLIST])), LOWLIM:DERIVDEGREE(INTEQN,PX,XVAR), ICLIST:REVERSE(ICLIST), IF LOWLIM=0 THEN (ICLIST:[], LOWLIM:3, IF (POFX:MYSOLVE(INTEQN,PX))#[] THEN (INTEQN:FIRST(POFX), KIND:[])), IF LOWLIM>2 OR (POFX:ODE2(INTEQN,PX,XVAR))=FALSE THEN RETURN(PRINTSOLN(INTEQN, 'VLFRNK, ICLIST, KIND)), POFX:RATSIMP(POFX), IF INITIAL AND LENGTH(ICLIST)-1=LOWLIM AND (INTEQN:ERRCATCH( APPLY(IF LOWLIM=1 THEN 'IC1 ELSE 'IC2, CONS(POFX,ICLIST))))#[] THEN (POFX:FIRST(INTEQN), ICLIST:[]), IF OPR(POFX)="=" AND LHS(POFX)=PX THEN (POFX:RHS(POFX), KIND:[]), IF ICLIST=[] AND KIND#[] THEN ICLIST:0, PRINTSOLN(POFX, 'VLFRNK, ICLIST, KIND))$ /* OBTAINS INITIAL CONDITIONS. GLOBALS: INTEQN, XVAR, LOWLIM, POFX, ICLIST, INITIAL */ GETICS():=BLOCK([VAL, INIT], VAL:AT(INTEQN, XVAR=LOWLIM), INIT:AT(POFX, XVAR=LOWLIM), INIT:MYSOLVE(VAL, INIT), IF INIT#[] THEN ( INIT:FIRST(INIT), ATVALUE(POFX, XVAR=LOWLIM, INIT), ICLIST:CONS(POFX=INIT, ICLIST)) ELSE INITIAL:FALSE)$ /* CONVERTS INTEGRANDS IN FUN TO FINITERANK FORM. UPPER LIMIT MUST BE XVAR AND LOWER LIMIT MUST BE A CONSTANT. GLOBALS: INITIAL, LOWLIM, UNKLIST, XVAR */ VCONVERT(FUN):= IF MAPATOM(FUN) THEN FUN ELSE IF OPR(FUN)#'?%INTEGRATE THEN MAP('VCONVERT, FUN) ELSE IF NOT FREEOF(XVAR,INPART(FUN,3)) OR INPART(FUN,4)#XVAR THEN THROW(FALSE) ELSE BLOCK([INTGR, NEWFUN, INT], IF LOWLIM='LOWLIM THEN LOWLIM:INPART(FUN,3) ELSE IF LOWLIM#INPART(FUN,3) THEN INITIAL:FALSE, INTGR:SUMFACTORS(ARG(FUN), INPART(FUN,2)), IF INTGR=[] THEN THROW(FALSE), NEWFUN:0, FOR TERM IN INTGR DO (INT:SUBSTINPART(TERM[2], FUN, 1), IF NOT MEMBER(INT,UNKLIST) THEN UNKLIST:CONS(INT, UNKLIST), NEWFUN:NEWFUN+TERM[1]*INT), RETURN(NEWFUN))$ /* LAPLACE TRANSFORM. GLOBALS: EQNLIST, PXLIST, XVAR */ TRANSFORM():=BLOCK([TEQNLIST, TRANSLIST, %S, PS, FLAG], PS:LAMBDA([FUN], LAPLACE(FUN, XVAR, %S)), FLAG:FALSE, TEQNLIST:MAPLIST('PS,EQNLIST), TRANSLIST:MAPLIST('PS,PXLIST), FOR SOLN IN SOLVESYS(TEQNLIST,TRANSLIST) DO IF FREEOF('?%INTEGRATE, '?%LAPLACE, SOLN) THEN (SOLN:MAPLIST(LAMBDA([FUN], ILT(FUN,%S,XVAR)), SOLN), PS:FREEOF('?%ILT, SOLN), PRINTSOLN(SOLN, 'TRANSFORM, IF PS THEN [] ELSE 0, IF PS THEN [] ELSE 'INCOMPLETE), FLAG:FLAG OR PS), RETURN(FLAG))$ /* COLLOCATION. GLOBALS: EQNLIST, PXLIST, NAPPROX, XVAR */ COLLOCATE():=BLOCK([LOWLIM, HIGHLIM, UNKLIST, %C, ELIST, FORM, INCR, POINT, NAME, LISTEQNS, SOLNLIST], APPLY('LOCAL, CONS('%C,MAPLIST('OPR,PXLIST))), LOWLIM:'MINF, HIGHLIM:'INF, MAP('GETLIMITS, EQNLIST), IF HIGHLIM<=LOWLIM THEN THROW(FALSE), SOLNLIST:LISTEQNS:UNKLIST:[], FOR UNK IN PXLIST DO ( NAME:OPR(UNK), FORM:APPROX(NAME,NAPPROX,XVAR), APPLY('DEFINE,[UNK,FORM]), SOLNLIST:CONS(FORM,SOLNLIST), FOR PARM:0 THRU NAPPROX-1 DO UNKLIST:CONS(%C[NAME,PARM],UNKLIST)), ELIST:apply('EV,[EQNLIST,'INTEGRATE,'EXPAND]), IF NOT FREEOF('?%INTEGRATE,ELIST) THEN THROW(FALSE), IF FREEOF(XVAR,ELIST) THEN LISTEQNS:ELIST ELSE (POINT:LOWLIM, INCR:IF NAPPROX>1 THEN (HIGHLIM-LOWLIM)/(NAPPROX-1) ELSE 0, FOR I THRU NAPPROX DO (LISTEQNS:APPEND(SUBST(POINT,XVAR,ELIST),LISTEQNS), POINT:POINT+INCR)), SOLVEANDSUBST(LISTEQNS,UNKLIST,REVERSE(SOLNLIST),'COLLOCATE))$ /* OBTAINS LARGEST LOWER LIMIT AND LEAST UPPER LIMIT FOR COLLOCATION POINTS GLOBALS: LOWLIM, HIGHLIM */ GETLIMITS(EXPR):= IF NOT MAPATOM(EXPR) THEN IF OPR(EXPR)#'?%INTEGRATE THEN FOR SUB IN EXPR DO GETLIMITS(SUB) ELSE BLOCK([LOW,HIGH], LOW:INPART(EXPR,3), HIGH:INPART(EXPR,4), IF NOT NUMBERP(LOW) OR NOT NUMBERP(HIGH) THEN THROW(FALSE), LOWLIM:MAX(LOW,LOWLIM), HIGHLIM:MIN(HIGH,HIGHLIM))$ /* APPROXIMATION FUNCTION FOR UNKNOWN FUNCTION FUN(VAR) */ APPROX(FUN,NPARMS,VAR):=SUM(%C[FUN,I-1]*VAR^(I-1),I,1,NPARMS)$ /* ---- THE FOLLOWING FUNCTIONS APPLY ONLY TO 2ND KIND EQNS */ /* FIXED-LIMIT, FINITE-RANK, 2ND KIND GLOBALS: EQNLIST, PXLIST */ FLFRNK2ND():=BLOCK([FRLIST, UNKLIST, UEQNLIST, EQNO, %C], APPLY('LOCAL, CONS('%C,MAPLIST('OPR, PXLIST))), UNKLIST:UEQNLIST:[], EQNO:0, FRLIST:MAPLIST('FCONVERT, EQNLIST), MAPLIST('DEFINE, PXLIST, FRLIST), UEQNLIST:apply('EV,[UEQNLIST,'INTEGRATE]), SOLVEANDSUBST(UEQNLIST, UNKLIST, FRLIST, 'FLFRNK2ND))$ /* RETURNS RESULT OF REPLACING ALL OCCURENCES OF 'INTEGRATE(F(X,U),U,A,B) IN FUN WITH SUM(Q[J](X)*%C[J],J,1,N) WHERE %C[J] IS 'INTEGRATE(R[J](U),U,A,B) AFTER EXPRESSING INTEGRANDS IN FINITE-RANK FORM (HERE F(X,U) BECOMES SUM(Q[J](X)*R[J](U),J,1,N). GLOBALS: EQNO (CURRENT NUMBER OF SUBSCRIPTED UNKNOWNS %C), UEQNLIST (LIST OF EQUATIONS RELATING %C'S) UNKLIST (LIST OF ALL %C UNKNOWNS TO BE SOLVED FOR) XVAR (INDEPEDENT VARIABLE) */ FCONVERT(FUN):= IF MAPATOM(FUN) THEN FUN ELSE IF OPR(FUN)#'?%INTEGRATE THEN MAP('FCONVERT, FUN) ELSE IF NOT FREEOF(XVAR,INPART(FUN,3)) OR NOT FREEOF(XVAR, INPART(FUN,4)) THEN THROW(FALSE) ELSE BLOCK([NEWFUN, INTGRND], INTGRND:SUMFACTORS(ARG(FUN), INPART(FUN,2)), IF INTGRND=[] THEN THROW(FALSE), NEWFUN:0, FOR TERM IN INTGRND DO (EQNO:EQNO+1, NEWFUN:NEWFUN+%C[EQNO]*TERM[1], UEQNLIST:CONS(%C[EQNO]-SUBSTINPART(TERM[2],FUN,1), UEQNLIST), UNKLIST:CONS(%C[EQNO], UNKLIST)), RETURN(NEWFUN))$ /* TAYLOR SERIES. GLOBALS: EQNLIST, PXLIST, XVAR, NAPPROX */ TAILOR():=BLOCK([XHAT, UFUN, EQTN, TFUN, NEQNS, VLIST, ORDER, VALUE, FACT], APPLY('LOCAL,APPEND([UFUN,EQTN,TFUN],MAPLIST('OPR,PXLIST))), MAP('GETXHAT, EQNLIST), NEQNS:0, VLIST:PXLIST, FOR EXPR IN EQNLIST DO ( NEQNS:NEQNS+1, EQTN[NEQNS]:EXPR, UFUN[NEQNS]:FIRST(VLIST), VLIST:REST(VLIST), TFUN[NEQNS]:VALUE:SUBST(XHAT,XVAR,EXPR), ATVALUE(UFUN[NEQNS], XVAR=XHAT, VALUE)), FACT:1, ORDER:NAPPROX, FOR I THRU NAPPROX DO ( FACT:FACT*I, FOR J THRU NEQNS DO ( VALUE:ERRCATCH(DIFF(EQTN[J],X)), IF VALUE=[] THEN EXITFOR(ORDER:I-1), EQTN[J]:FIRST(VALUE), VALUE:ERRCATCH(RATSIMP(AT(EQTN[J],XVAR=XHAT))), IF VALUE=[] THEN EXITFOR(ORDER:I-1), UFUN[J]:DIFF(UFUN[J],X), VALUE:FIRST(VALUE), ATVALUE(UFUN[J], XVAR=XHAT, VALUE), TFUN[J]:TFUN[J]+VALUE*(X-XHAT)^I/FACT), IF ORDER#NAPPROX THEN EXITFOR(FALSE)), FOR I:NEQNS STEP -1 THRU 1 DO VLIST:CONS(TFUN[I],VLIST), PRINTSOLN(VLIST, 'TAILOR, ORDER, TESTSOLN(VLIST)))$ /* OBTAINS EXPANSION POINT FOR TAYLOR SERIES BY SOLVING B(X)=A(X). GLOBALS: XVAR, XHAT */ GETXHAT(EXPR):= IF MAPATOM(EXPR) THEN FALSE ELSE IF OPR(EXPR)#'?%INTEGRATE THEN FOR SUB IN EXPR DO GETXHAT(SUB) ELSE IF XHAT='XHAT THEN (XHAT:MYSOLVE(INPART(EXPR,3)-INPART(EXPR,4),XVAR), IF XHAT=[] THEN THROW(FALSE) ELSE XHAT:FIRST(XHAT), IF RATSUBST(XHAT, XVAR, INPART(EXPR,3))#XHAT THEN THROW(FALSE)) ELSE IF SUBST(XHAT,XVAR,EXPR)#0 THEN THROW(FALSE)$ /* FREDHOLM-CARLEMAN SERIES. GLOBALS: NAPPROX, KXU, AX, BX, FX, XVAR, UVAR */ FREDSERIES():=BLOCK([ORDER, ALPHA, BTA, TOP, BOT, KIND, TVAR, EQTN, KXT], EQTN:FIRST(EQNLIST), ALPHA:CATCH(FINDINT(EQTN)), IF ALPHA=FALSE THEN RETURN(FALSE), UVAR:INPART(ALPHA,2), BTA:LIN(EQTN,ALPHA), IF BTA=FALSE THEN RETURN(FALSE), KXU:LIN(ARG(ALPHA),SUBST(UVAR,XVAR,PX)), IF KXU=FALSE THEN RETURN(FALSE), AX:INPART(ALPHA,3), BX:INPART(ALPHA,4), FX:BTA[2], IF KXU[2]#0 THEN FX:FX+MYINT(BTA[1]*KXU[2],UVAR,AX,BX), KXU:KXU[1]*BTA[1], ORDER:NAPPROX, KIND:'APPROXIMATE, KXT:SUBST(TVAR,UVAR,KXU), BOT:BTA:1, TOP:ALPHA:RAT(KXU,UVAR,XVAR), FOR I THRU NAPPROX DO ( BTA:-MYINT(RATSUBST(UVAR,XVAR,ALPHA),UVAR,AX,BX)/I, ALPHA:BTA*KXU+MYINT(KXT*RATSUBST(TVAR,XVAR,ALPHA),TVAR,AX,BX), BOT:BOT+BTA, IF ALPHA=0 THEN EXITFOR(KIND:[]), TOP:TOP+ALPHA), KXT:FX+MYINT(SUBST(UVAR,XVAR,FX)*TOP,UVAR,AX,BX)/RATDISREP(BOT), PRINTSOLN(KXT, 'FREDSERIES, ORDER, KIND))$ /* NEUMANN SERIES. GLOBALS: EQNLIST, PXLIST, GUESSLIST, NAPPROX */ NEUMANN():=BLOCK([ORDER, NEWGUESS], APPLY('LOCAL, MAPLIST('OPR,PXLIST)), ORDER:NAPPROX, IF GUESSLIST='NONE THEN GUESSLIST:MAPLIST('ZEROINT, EQNLIST), FOR COUNT THRU NAPPROX DO (MAPLIST('DEFINE, PXLIST, GUESSLIST), NEWGUESS:apply('EV,[EQNLIST,'INTEGRATE]), IF NOT FREEOF('?%INTEGRATE, NEWGUESS) THEN EXITFOR(ORDER:COUNT-1), GUESSLIST:MAPLIST('RATSIMP,NEWGUESS)), PRINTSOLN(GUESSLIST, 'NEUMANN, ORDER, TESTSOLN(GUESSLIST)))$ /* ---- THE FOLLOWING FUNCTIONS APPLY ONLY TO 1ST KIND EQNS */ /* FIXED-LIMIT FINITE-RANK 1ST KIND. GLOBALS: KXU, FX, UVAR, XVAR, AX, BX */ FLFRNK1ST():=BLOCK([SF, CNT, INTG, FORM, %C, UNKLIST, EQLIST,RES], LOCAL(%C), CNT:FORM:0, UNKLIST:[], IF NOT FREEOF(XVAR,[AX,BX]) THEN RETURN(FALSE), INTG:SUMFACTORS(KXU,UVAR), FOR TERM IN INTG DO (CNT:CNT+1, FORM:FORM+TERM[2]*%C[CNT], UNKLIST:CONS(%C[CNT],UNKLIST)), EQLIST:[], SF:NUM(RATSIMP(MYINT(KXU*FORM,UVAR,AX,BX)-FX,XVAR)), RES:0, IF OPR(SF)#"+" THEN SF:[SF], FOR TERM IN SF DO IF OPR(TERM)="*" THEN (INTG:PARTITION(TERM,XVAR), IF INTG[2]=1 THEN RES:RES+TERM ELSE EQLIST:CONS(INTG[1],EQLIST)) ELSE IF FREEOF(XVAR,TERM) THEN RES:RES+TERM ELSE ERROR("ERROR 1"), IF RES#0 THEN EQLIST:CONS(RES,EQLIST), IF LENGTH(EQLIST)#CNT THEN ERROR("ERROR 2"), SOLVEANDSUBST(EQLIST, UNKLIST, SUBST(XVAR,UVAR,FORM), 'FLFRNK1ST))$ /* SOME=ONE*PEACE+TWO, WHERE PEACE COULD BE A PRODUCT. */ MYDIVIDE(SOME,PEACE,VAR):=BLOCK([RES, ONE, TWO], ONE:TWO:0, IF OPR(SOME)#"+" THEN SOME:[SOME], FOR PRT IN SOME DO (RES:QUOTIENT(PRT,PEACE), IF RES#0 AND FREEOF(VAR,RES) THEN ONE:ONE+RES ELSE TWO:TWO+PRT), RETURN([ONE,TWO]))$ /* GENERALIZED ABEL METHOD. GLOBALS: KXU, AX, BX, FX, XVAR, UVAR */ ABEL():=BLOCK([POWER,DEN,FUN], IF NOT FREEOF(XVAR,AX) OR BX#XVAR OR OPR(KXU)#"^" THEN THROW(FALSE), POWER:-INPART(KXU,2), DEN:INPART(KXU,1), IF NOT FREEOF(UVAR, POWER) OR OPR(DEN)#"+" THEN THROW(FALSE), FUN:PARTITION(DEN,UVAR), IF NOT FREEOF(XVAR, FUN[2]) OR RATSIMP(SUBST(UVAR,XVAR,FUN[1])+FUN[2])#0 THEN THROW(FALSE), IF SIGN(POWER)#'POS OR SIGN(POWER-1)#'NEG THEN THROW(FALSE), FUN:MYINT(DIFF(-FUN[2],UVAR)*SUBST(UVAR,XVAR,FX)*DEN^(POWER-1),UVAR,AX,BX), FUN:SIN(%PI*POWER)/%PI*DIFF(FUN,XVAR), PRINTSOLN(FUN, 'ABEL, [], []))$ /* FIRSTKINDSERIES. GLOBALS: PX, NAPPROX, GUESS, EQN, PX */ FIRSTKINDSERIES():=BLOCK([KIND, ORDER, CORRECTION, TRIAL], APPLY('LOCAL,[OPR(PX)]), KIND:'APPROXIMATE, ORDER:NAPPROX, TRIAL:IF GUESS='NONE THEN ZEROINT(EQN) ELSE GUESS, FOR I THRU NAPPROX DO (APPLY('DEFINE,[PX, TRIAL]), CORRECTION:apply('EV,[EQN,'INTEGRATE]), IF NOT FREEOF('?%INTEGRATE, CORRECTION) THEN EXITFOR(ORDER:I-1), IF (CORRECTION:RATSIMP(CORRECTION))=0 THEN EXITFOR(KIND:[]), TRIAL:RATSIMP(TRIAL-CORRECTION)), PRINTSOLN(TRIAL, 'FIRSTKINDSERIES, ORDER, KIND))$ /* KILL(LABELS)$ */ TTYOFF:FALSE$