(* ************************************************************************ *) (* : Title : PDESpecialSolutions.m * : Authors : Douglas Baldwin, Unal Goktas, and Willy Hereman. * Department of Mathematical and Computer Sciences * Colorado School of Mines * Golden, Colorado 80401-1887, U.S.A * Contact email: whereman@mines.edu or douglas@douglasbaldwin.com * : Last updated : Monday, February 28, 2005 at 12:38 AM by Douglas Baldwin * : Summary : Solves systems of nonlinear PDEs in terms of the hyperbolic * functions Tanh, Sech, and combinations of Sech and Tanh. * : Also computes exact solutions in terms of the Jacobi elliptic Cn and Sn functions. * : Warnings : This package uses the assumption that Sqrt[a^2] -> a, etc. * (see mySimplificationRules below) when simplifying. *) (* ************************************************************************ *) (* Algorithms and implementation by Douglas Baldwin, *) (* Unal Goktas (WRI) and Willy Hereman. *) (* Copyright 2001 *) BeginPackage["Calculus`PDESpecialSolutions`"] Unprotect[PDESpecialSolutions] PDESpecialSolutions::usage = "PDESpecialSolutions[eqns, funcs, vars, params, opts] solves a system of nonlinear partial differential equations (PDEs) for funcs, with independent variables vars and non-zero parameters params." PDESpecialSolutions::poly = "This system must be a polynomial of fixed order." PDESpecialSolutions::form = "The user entered form is not valid." PDESpecialSolutionsmSolver::freedom = "Freedom exists in this system, as `1` are both dominant powers in expression `2`." PDESpecialSolutions`dominantBehavior::freevalues = "The solution(s) `1` do not fix all the values of \!\(M\_i\). We will now generate all the possible values for the free \!\(M\_i\) upto the user option MaxDegreeOfThePolynomial which is set to 3 by default, but can be any non-negative integer. " PDESpecialSolutions`dominantBehavior::underdetermined = "There is too much freedom in choosing the values of \!\(M\_i\). We will now generate all the possible values for the free \!\(M\_i\) from zero to the user option MaxDegreeOfThePolynomial (with default 3), but can be any integer." PDESpecialSolutionsmSolver::DegreeOfThePolynomial = "Algorithm continuing with user given degree of the polynomial `1`." PDESpecialSolutionsmSolver::remove = "The potential solutions `1` are being removed because they are (i) negative, (ii) contain freedom, (iii) fail to balance highest exponent terms from two different terms in the origional system. If \!\(M\_i\) < 0, then the transform u -> 1/v might result in a system that PDESpecialSolutions may be able to solve." PDESpecialSolutionsBuildSystem::fail = "The systems of equations was inconsistant under the assumption that \!\(a\_\(i,M\_i\) \[NotEqual] 0\), please check the \!\(M\_i\) values by hand." PDESpecialSolutions::solve = "The algorithm failed while attempting to find the coefficients for the polynomial solution." (* Options definition. *) Options[PDESpecialSolutions] = { Verbose -> False, Form -> tanh, NumericTest -> True, SymbolicTest -> False, InputForm -> False, DegreeOfThePolynomial -> {}, MinDegreeOfThePolynomial -> 1, MaxDegreeOfThePolynomial -> 3, SolveAlgebraicSystem -> True, HighestOrderFirst -> False, mSolverVerbose -> False }; sech=.; sechtanh=.; cn=.; sn=.; jacobisn=.; jacobicn=.; Begin["`Private`"] (* If called with non-lists, makes them lists. *) PDESpecialSolutions[eqns_?(!ListQ[#]&), funcs_, vars_, param_, opts___?OptionQ]:= PDESpecialSolutions[{eqns}, funcs, vars, param, opts] PDESpecialSolutions[eqns_, funcs_?(!ListQ[#]&), vars_, param_, opts___?OptionQ]:= PDESpecialSolutions[eqns, {funcs}, vars, param, opts] PDESpecialSolutions[eqns_, funcs_, vars_?(!ListQ[#]&), param_, opts___?OptionQ]:= PDESpecialSolutions[eqns, funcs, {vars}, param, opts] PDESpecialSolutions[eqns_, funcs_, vars_, param_?(!ListQ[#]&), opts___?OptionQ]:= PDESpecialSolutions[eqns, funcs, vars, {param}, opts] (* ************************************************************************ *) (* Start of function PDESpecialSolutions. *) PDESpecialSolutions[eqns_List, funcs_List, vars_List, param_List, opts___?OptionQ]:= Block[ {testarg = FreeQ[MapAll[Expand, eqns], (Power[b__, a__] /; (!FreeQ[a, _Symbol] | MemberQ[b, _Real | _Integer]))], theForm = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn -> sn, jacobicn -> cn}, VerboseTest = Verbose /. {opts} /. Options[PDESpecialSolutions], time = TimeUsed[], system, mSoln, newSystem, solution }, (* Protected Local Variables *) If[!testarg, Message[PDESpecialSolutions::poly]]; ( (* If verbose, prints system. *) If[VerboseTest, Print["The given system of differential equations is: "]; Print[ TableForm[eqns /. {Derivative[n__][F_][x__]:> SequenceForm[F, Subscript[SequenceForm @@ Flatten[ Table[#[[1]], {#[[2]]} ]& /@ Transpose[{{x}, {n}}] ] ] ], Sequence @@ ((#->Head[#])& /@ funcs) } ] ] ]; (* Step 1 in the paper. *) (* Transforms the PDE into a nonlinear ODE in $T=Tanh$ or $S=Sech$ *) If[VerboseTest, Print["Transform the differential equation(s) into a nonlinear ODE in "<> ToString[ (theForm /. {tanh -> "T=Tanh", sech -> "S=Sech", sechtanh -> "S=Sech", cn -> "CN=JacobiCN", sn -> "SN=JacobiSN"} ) ] ] ]; system = PDESpecialSolutionsVarChange[ eqns /. (a__==b__):>(a-b), funcs, vars, opts ]; If[VerboseTest, Print[CleanUpFunction[system /. myTrackingVariable[_]->1]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* Step 2 in the paper. *) (* Determines the maximal degree of the polynomial solutions. *) If[VerboseTest, Print["Determine the maximal degree of the polynomial solutions."] ]; mSoln = Internal`DeactivateMessages[ PDESpecialSolutionsmSolver[system, theForm, opts], Solve::svars ]; If[VerboseTest, Print[CleanUpFunction[mSoln]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* Step 3 in the paper. *) (* Derives the algebraic system for the coefficients $a_{ij}, b_{ij}$. *) If[VerboseTest, Print["Derive the nonlinear algebraic system for the coefficients."] ]; {mSoln, newSystem} = PDESpecialSolutionsBuildSystem[mSoln, system, vars, param, opts]; If[VerboseTest, Print["The nonlinear algebraic system is"]; Print[CleanUpFunction[newSystem /. myTrackingVariable[_]->1]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* DB:3/21/2003 Output algebraic system. *) If[!(SolveAlgebraicSystem /. {opts} /. Options[PDESpecialSolutions]), Return[CleanUpFunction[newSystem]]; ]; (* Step 4 in the paper. *) (* Solve the nonlinear parameterized algebraic system. *) If[VerboseTest, Print["Solve the nonlinear algebraic system."] ]; solution = Algebra`AnalyzeAndSolve`AnalyzeAndSolve @@ #& /@ newSystem; If[Length[Flatten[solution]] === 0, Message[PDESpecialSolutions::solve];Return[{}]; ]; (* Reformat the output after Unal's change of AnalyzeAndSolve's output *) (* by Douglas Baldwin on 11/27/2002 *) solution = (# /. List[(a_List)..] :> a)& /@ #& /@ solution; If[VerboseTest, Print[CleanUpFunction[solution]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* Step 5 in the paper. *) (* Build and test the solitary wave solutions. *) If[VerboseTest, If[( !(NumericTest /. {opts} /. Options[PDESpecialSolutions]) && !(SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), Print["Build the travelling wave solutions."], Print["Build and (numerically and/or symbolically) test "<> "the solutions."] ] ]; solution = PDESpecialSolutionsBuildSolutions[solution, mSoln, eqns, funcs, vars, param, opts]; If[VerboseTest, Print["Time Used:", TimeUsed[]-time]; Print[]; If[( (NumericTest /. {opts} /. Options[PDESpecialSolutions]) || (SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), Print["The (numerically and/or symbolically) tested final solutions:"], Print["The list of possible solutions is:"] ]; ]; (* Returns solutions to user. *) Return[ If[InputForm /. {opts} /. Options[PDESpecialSolutions], InputForm[CleanUpFunction[solution]], CleanUpFunction[solution] ] ] ) /; testarg ] (* ************************************************************************ *) (* : Title : PDESpecialSolutionsVarChange *) (* : Author : Douglas Baldwin *) (* : Input : The equation list (F[i][x,y,z,t]) and the form tanh. *) (* : Output : The system in the form u[i][T] depending on whether it is tanh. *) (* : Summary : Converts from real space-time domain to u[i][T]. *) PDESpecialSolutionsVarChange[eqns_List, funcs_List, vars_List, opts___?OptionQ]:= Block[{Form, (* The form of the wanted solution. *) funcRules, (* conversion from user functions *) i=0, (* Iterator for myTrackingVariables *) eqList, (* The list of equations, in a usable form *) system, (* The system after transformation to ODE in T *) phi (* A temporary function. *) }, (* Protected Local Variables *) (* Sets the Form from user options *) Form = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn -> sn, jacobicn -> cn}; (* If form isn't correct, it sends a message and aborts. *) If[!(Or @@ ((Form === #)& /@ {tanh, sech, sechtanh, cn, sn})), Message[PDESpecialSolutions::form]; Abort[] ]; (* Creates a set of rules to converts the users functions *) (* (e.g., u[n,t]) to varChangeFunction[i][n,t]. *) funcRules = Table[ Head[funcs[[i]] ] -> varChangeFunction[Form][i], {i, Length[funcs]} ]; (* Adds tacking variables which will be used in the mSolve *) (* function to check for false balances coming form a single *) (* term. *) eqList = If[Head[#] === Plus, Plus @@ ((myTrackingVariable[++i]*#) & /@ List @@ #), myTrackingVariable[++i]*#] & /@ Expand[eqns]; (* Converts the user functions to varChangeFunction[i]. *) eqList = eqList //. funcRules; (* Changes the space of the functions. *) system = eqList /. { varChangeFunction[form_][n_][var__] :> phi[form][n][ Sum[c[i]*{var}[[i]], {i, Length[{var}]} ] ], Derivative[temp__][varChangeFunction[form_][n_]][var__] :> ( D[ phi[form][n][ Sum[c[i]*{var}[[i]], {i, Length[{var}]} ] ], Sequence @@ Transpose[{{var}, {temp}}] ] ) }; (* Parameterizes the sequence $c_1 x+c_2 y+...\to\eta$ *) system = system /. { phi[form_][n_][__] :> phi[form][n][eta], Derivative[i_][phi[form_][n_]][__] :> Derivative[i][phi[form][n]][eta] }; (* Repeatedly applies the chain rule. *) system = system /. {Derivative[a_][phi[tanh][n_]][eta] :> Expand[ Nest[(1-T^2)*D[#, T]&, u[n][T], a ] ], phi[tanh][n_][eta] :> u[n][T], Derivative[a_][phi[sech][n_]][eta] :> Expand[ Nest[-(S*Sqrt[1-S^2])*D[#, S]&, u[n][S], a ] ], phi[sech][n_][eta] :> u[n][S], Derivative[a_][phi[sechtanh][n_]][eta] :> Expand[ Nest[-(S*Sqrt[1-S^2])*D[#, S]&, u[n][S], a ] ], phi[sechtanh][n_][eta] :> u[n][S], Derivative[a_][phi[cn][n_]][eta] :> Expand[ Nest[-Sqrt[1-mod-CN^2+2*mod*CN^2-mod*CN^4]*D[#, CN]&, u[n][CN], a ] ], phi[cn][n_][eta] :> u[n][CN], Derivative[a_][phi[sn][n_]][eta] :> Expand[ Nest[Sqrt[1-SN^2-mod*SN^2+mod*SN^4]*D[#, SN]&, u[n][SN], a ] ], phi[sn][n_][eta] :> u[n][SN] }; system = MapAll[Expand, system]; (* Simplifies the system by removing extra (1-T^2)'s *) (* and T's from the system. *) system = If[(# /. T->-1) === 0 && (# /. T->1) === 0, Factor[#/(1-T^2)], Factor[#]]& /@ system; system = If[(# /. T -> 0) === 0, Factor[#/T], Factor[#]]& /@ system; system = If[(# /. S -> 0) === 0, Factor[#/S], Factor[#]]& /@ system; system = If[(# /. S -> -1) === 0 && (# /. S->1) === 0, Factor[#/Sqrt[1-S^2]], Factor[#]]& /@ system; system = If[(# /. CN -> 0) === 0, Factor[#/CN], Factor[#]]& /@ system; system = If[(((# /. CN -> -1) === 0 && (# /. CN -> 1) === 0 ) && (Factor[# /. CN -> Sqrt[-1+mod]/Sqrt[mod] ]) === 0), Factor[#/Sqrt[1-mod-CN^2+2*mod*CN^2-mod*CN^4]], Factor[#]]& /@ system; system = If[(# /. SN -> 0) === 0, Factor[#/SN], Factor[#]]& /@ system; system = If[(((# /. SN -> -1) === 0 && (# /. SN -> 1) === 0 ) && (Factor[# /. SN -> 1/Sqrt[mod] ]) === 0 ), Factor[#/Sqrt[1-SN^2-mod*SN^2+mod*SN^4]], Factor[#]]& /@ system; (* Returns the system back to the PDESpecialSolutions *) Return[MapAll[Factor, system]]; ]; (* ************************************************************************ *) (* : Title : PDESpecialSolutionsmSolver *) (* : Author : Douglas Baldwin *) (* : Date : Sunday, December 15, 2002 *) (* : Input : The system generated by varChange. *) (* : Output : The power of the polynomial solutions for u[i][T], id est u[i][T] = a[1,0]+a[1,1]*T+a[1, 2]*T^2+...+a[1,n]*T^n. *) PDESpecialSolutionsmSolver[system_List, theForm_, opts___?OptionQ] := Module[{mSoln}, mSoln = { ToExpression[ StringReplace[ ToString[DegreeOfThePolynomial /. {opts}], {"m" -> "Calculus`PDESpecialSolutions`Private`m", "n" -> "Calculus`PDESpecialSolutions`Private`n" } ] ] }; Message[PDESpecialSolutionsmSolver::DegreeOfThePolynomial, CleanUpFunction[mSoln] ]; Return[mSoln]; ] /; (DegreeOfThePolynomial /. {opts} /. Options[PDESpecialSolutions]) =!= {}; PDESpecialSolutionsmSolver[system_List, theForm_, opts___?OptionQ] := Module[{mSolverDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), myTrackingVariableMax, (* Max Tracking Variable. *) funcRules, (* Replaces F with T, S, CN, SN, etc. *) eqnList, (* The result of the system generation func. *) mList0, mList, (* Lists of the powers of F in m_i. *) myMList, (* List of m_i's: {m_1, m_2, m_3, ...} *) mRules, (* Rules for m_i, such as m_1 -> m_2. *) mSoln, mSoln0 (* List of explicit solutions for m_i. *) }, (* Determines max{i : myTrackingVariables[i]} by *) (* applying a rule to system which sets myTrackingVariable *) (* in the process. *) myTrackingVariableMax = 0; system /. myTrackingVariable[a_Integer] :> myTrackingVariable[ If[a > myTrackingVariableMax, myTrackingVariableMax = a ]; a]; (* Sets up a list of rules for replacing the form of F. *) funcRules = {tanh -> T, sech -> S, sechtanh -> S, cn -> CN, sn -> SN}; (* Breaks into parts P + Q Sqrt[...] and substitutes *) (* the ansatz chi F^{m_i} for u,v,...,w *) eqnList = mSolve`SystemGeneration[system, theForm, myTrackingVariableMax, opts]; If[mSolverDebug, Print["eqnList, System after ansatz substitution:"]; Print[CleanUpFunction[eqnList]]; ]; (* Pulls off the expressions for m_i *) mList0 = mSolve`ListFormation[eqnList, theForm, funcRules, myTrackingVariableMax, opts]; If[mSolverDebug, Print["mList0, Exponents of F before simplification:"]; Print[CleanUpFunction[mList0]]; ]; (* If the form is of SechTanh, it uses the solution m_i=2, n_i=1 *) If[theForm === sechtanh, (* Sets mSoln so m_i = 2 and n_i = 1 *) mSoln = { Flatten[ Table[{m[i]->2, n[i]->1}, {i, Length[system]} ] ] }; (* Returns the solution. *) Return[mSoln]; ]; (* Removes expressions which cannot be the highest power *) (* by trackingVariable. *) mList = Flatten[mSolve`Simplification[#, opts]& /@ #]& /@ mList0; If[mSolverDebug, Print["mList, mList0 after simplification:"]; Print[CleanUpFunction[mList]]; ]; (* Forms a list of all the m_i. *) myMList = {}; mList /. m[i_Integer] :> (myMList = Append[myMList, m[i]];m[i]); myMList = Union[myMList]; If[mSolverDebug, Print["The $m_i$s to be solved for are:"]; Print[CleanUpFunction[myMList]]; ]; (* Solves the expressions of m_i for m_i *) mRules = mSolve`RulesSolver[#, myMList, opts]& /@ mList; If[mSolverDebug, Print["mRules, the solution iteration yields:"]; Print[CleanUpFunction[mRules]]; ]; (* Uses the previous results to determine *) (* explicit solutions for m_i. *) mSoln0 = mSolve`PowerSolver[mRules, myMList, opts]; If[mSolverDebug, Print["mSoln, the possible solutions before pruning:"]; Print[CleanUpFunction[mSoln0]]; ]; (* DB:5/15/2004 New version from Painleve code. *) mSoln = Join[ mSoln0, mSolve`FixFreeM[mSoln0, myMList, opts]]; (* Remove bad solutions. *) mSoln = mSolve`SystemCleanUp[eqnList, mList, mSoln, myMList, theForm /. funcRules, opts ]; (* DB:3/21/2003 Warn the user when potential solutions are removed. *) If[Length[Complement[Union[Sort /@ mSoln0], mSoln] ] > 0, StylePrint[ "The potential solutions "<> ToString[ InputForm[ CleanUpFunction[ Complement[Union[Sort /@ mSoln0], mSoln] ] ] ] <> " are being removed because they are (i) negative, (ii) contain "<> "freedom, (iii) fail to balance highest exponent terms from two "<> "different terms in the origional system. If \!\(M\_i\) < 0, "<> "then the transform u -> 1/v might result in a system that "<> "PDESpecialSolutions may be able to solve.", "Message" ]; ]; If[mSolverDebug, Print["mSoln, after being pruned:"]; Print[CleanUpFunction[mSoln]]; ]; (* If it doesn't find any find any solutions, it quits. *) If[Length[mSoln] === 0, StylePrint[ "The algorithm failed while attempting to find a positive values "<> "of \!\(M\_i\). The list of rules that constrain the system are " <> ToString[ InputForm[ CleanUpFunction[mRules] ] ] <> ". The origional powers in \!\(M\_i\) are " <> ToString[ InputForm[ CleanUpFunction[mList] ] ] <> ". The package will now test all the possible solutions in the " <> "range MinDegreeOfThePolynomial to MaxDegreeOfThePolynomial.", "Message" ]; (* Generate all the possible solutions from -3 to -1. *) mSoln = mSolve`SystemCleanUp[eqnList, mList, mSolve`GenerateAlternativeSolutions[ myMList, opts ], myMList, theForm /. funcRules, opts ] ]; (* If it doesn't find any find any solutions, it quits. *) If[Length[mSoln] == 0, Abort[ ] ]; If[mSolverDebug, Print["mSoln, the final version before being returned:"]; Print[CleanUpFunction[mSoln]]; ]; (* Returns the solutions. *) Return[mSoln]; ] /; (DegreeOfThePolynomial /. {opts} /. Options[PDESpecialSolutions]) === {}; (* Breaks up the system into the correct form. *) mSolve`SystemGeneration[system_, theForm_, myTrackingVariableMax_, opts___?OptionQ] := Block[{SystemGenerationDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), sqrtRules, (* A set of rules for replacing what's *) (* under the radical. *) mySystem, (* The system passed to the function after *) (* transformation. *) eqnList (* The result which will be returned. *) }, (* Sets up the way in which the system will be divided up *) (* in the next step, based on the Options set by the user. *) sqrtRules = {tanh -> Null, sech -> 1-S^2, sechtanh -> 1-S^2, cn -> 1-mod-CN^2+2*mod*CN^2-mod*CN^4, sn -> 1-SN^2-mod*SN^2+mod*SN^4 }; (* DB:2/28/2005 Added new Coefficient function, for 5.0 compatibility. *) myCoefficient[p_Plus, q_, 0] := Plus @@ Select[p, FreeQ[#,q]&]; myCoefficient[p_, q_, 0] := Plus @@ Select[{p}, FreeQ[#,q]&]; myCoefficient[p_Plus, q_, 1] := Plus @@ Cases[p, z_. q -> z]; myCoefficient[p_, q_, 1] := Plus @@ Cases[{p}, z_. q -> z]; (* Breaks up system into P and Q, where system = P + Sqrt[...]*Q. *) mySystem = {myCoefficient[#,Sqrt[theForm /. sqrtRules], 0]& /@ MapAll[Expand, system], myCoefficient[#,Sqrt[theForm /. sqrtRules], 1]& /@ MapAll[Expand, system]}; If[SystemGenerationDebug, Print["The system divided into P and Q, is: "]; Print[CleanUpFunction[mySystem]]; ]; mySystem = mSolve`SimplifySystem[mySystem]; If[SystemGenerationDebug, Print["The system after simplification is: "]; Print[CleanUpFunction[mySystem]]; ]; (* Replaces the function u[i_] with a function in *) (* F^m[i] or F^n[i]. *) eqnList = Expand[ DeleteCases[ Flatten[ mySystem /. u[i_] :> Function[{F}, \[Chi][i]*F^m[i] ] ], 0 ] ]; If[SystemGenerationDebug, Print["The system after substitution of $chi F^{m_i}$"]; Print[CleanUpFunction[eqnList]]; ]; (* Returns the list after substitution of the ansatz. *) Return[eqnList]; ]; (* Simplifies the system by removing extra (1-T^2), (1-S^2), *) (* S, T, etc. from the system. *) mSolve`SimplifySystem[system_]:= Module[{mySystem = system}, mySystem = If[(# /. T->-1) === 0 && (# /. T->1) === 0, Expand[#/(1-T^2)], Expand[#]]& /@ #& /@ mySystem; mySystem = If[(# /. T->0) === 0, Expand[#/T], Expand[#]]& /@ #& /@ mySystem; mySystem = If[(# /. S->0) === 0, Expand[#/S], Expand[#]]& /@ #& /@ mySystem; mySystem = If[(# /. S->-1) === 0 && (# /. S->1) === 0, Expand[#/Sqrt[1-S^2]], Expand[#]]& /@ #& /@ mySystem; mySystem = If[(# /. CN -> 0) === 0, Expand[#/CN], Expand[#]]& /@ mySystem; mySystem = If[(((# /. CN -> -1) === 0 && (# /. CN -> 1) === 0 ) && (Expand[# /. CN -> Sqrt[(-1+mod)/mod] ]) === 0 ), Expand[#/Sqrt[1-mod-CN^2+2*mod*CN^2-mod*CN^4]], Expand[#]]& /@ mySystem; mySystem = If[(# /. SN->0) === 0, Expand[#/SN], Expand[#]]& /@ mySystem; mySystem = If[(((# /. SN->-1) === 0 && (# /. SN->1) === 0 ) && (Expand[# /. SN -> 1/Sqrt[mod] ]) === 0 ), Expand[#/Sqrt[1-SN^2-mod*SN^2+mod*SN^4]], Expand[#]]& /@ mySystem; (* Returns the simplified system. *) Return[mySystem]; ]; mSolve`ListFormation[eqnList0_, Form_, funcRules_, myTrackingVariableMax_, opts___?OptionQ]:= Block[{ mPowerListFormationDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), eqnList = eqnList0, mList }, eqnList = ( Table[Coefficient[#, myTrackingVariable[i] ], {i, 1, myTrackingVariableMax} ]& /@ eqnList0) /. 0:>Sequence[]; If[mPowerListFormationDebug, Print["The system further divided, is: "]; Print[CleanUpFunction[eqnList]]; ]; (* Breaks up expressions into lists of terms. *) eqnList = If[Head[#]===Plus, List @@ #, {#}]& /@ #& /@ MapAll[Expand, eqnList]; If[mPowerListFormationDebug, Print["The system where + -> {} is:"]; Print[CleanUpFunction[eqnList]]; ]; (* Pulls off the exponents of T and forms a list *) (* of expressions of the form {{{1+m[1]},{..}..}..} *) (* DB:5/15/2004 moved Union. *) mList = ( Union[ Exponent[#, Form /. funcRules ] ]& /@ #)& /@ eqnList; If[mPowerListFormationDebug, Print["The list of powers of F is:"]; Print[CleanUpFunction[mList]]; ]; (* Returns the number of equations and the list of powers. *) Return[mList]; ]; mSolve`Simplification[mList0_, opts___?OptionQ]:= Module[{mPowerListSimplificationDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), mList, (* The list of powers. *) mListStructure (* Structure of powers. *) }, (* Protected Local Variables. *) (* The following lines breaks up the list of exponents *) (* of T and then removes invalid cases. *) (* Breaks up the expressions of $a+b \, m_i+c\, m_j+...$ *) (* where $a,b,c,...,i,j,k,...\in\mathbb{R}$ into lists. *) mList = If[Head[#]===Plus, List @@ #, If[(Head[#]===m || Head[#]===n) || Head[#]===Times, {#}, # ] ]& /@ mList0; If[mPowerListSimplificationDebug, Print["Breaks up $a+b m_i+c m_j+...$ into lists:"]; Print[CleanUpFunction[mList]]; ]; (* The following routine strips off the $a$ in the above *) (* expression and leaves only the underlying structure. *) mListStructure = Union[mList /. {a_Integer, b__}:>{b}]; If[mPowerListSimplificationDebug, Print["The above sans $a$:"]; Print[CleanUpFunction[mListStructure ]]; ]; (* Re-organizes the list of powers of T by the structure *) (* listed above. *) mList = Cases[mList, {_, Sequence @@ #} | #]& /@ mListStructure; If[mPowerListSimplificationDebug, Print["The list of powers reorganized by structure:"]; Print[CleanUpFunction[mList ]]; ]; (* Determines the maximum $a$ in each power of T *) mList = {Max[# /. {a_, ___}:>If[IntegerQ[a], a, 0]& /@ #]}& /@ mList; If[mPowerListSimplificationDebug, Print["The list of powers after powers which cannot be "<> "the highest power are removed:"]; Print[CleanUpFunction[mList ]]; ]; (* Creates a list of the maximum powers of T, such *) (* that all the members of the list are of the form *) (* $a_{\rm max}+b\,m_i+c\,m_j+...+d\,m_i\,m_j$ *) mList = (Plus @@ Flatten[#])& /@ Transpose[{mList, mListStructure}]; If[mPowerListSimplificationDebug, Print["Formulates the powers back into the correct form:"]; Print[CleanUpFunction[mList ]]; ]; (* Returns the simplified list. *) Return[mList]; ]; mSolve`RulesSolver[mList0_, myMList_, opts___?OptionQ]:= Module[{rulesSolverDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), mList, eqnList, mRules (* List of rules from first solve. *) }, (* Protected local variables. *) (* Makes sure it is the simplest list possible *) (* for combinatorial purposes. *) mList = Union[mSolve`Simplification[mList0] ]; If[rulesSolverDebug, Print["The mList used in Power Solver is"]; Print[CleanUpFunction[mList]]; ]; (* Forms a list of equations to be solved for m_i. *) eqnList = Flatten[ Map[ Thread[ Table[#[[1]], {Length[#] -1} ] == Drop[#, 1] ]&, Table[Drop[mList, i], {i, 0, Length[mList] - 2} ] ], 1 ]; If[rulesSolverDebug, Print["The Equations to be solved first"]; Print[CleanUpFunction[eqnList]]; ]; (* Does the first run of solving. *) mRules = Union[ Flatten[ Solve[#, myMList]& /@ eqnList ] ]; If[rulesSolverDebug, Print["The first set of solutions are"]; Print[CleanUpFunction[mRules]]; ]; Return[mRules]; ]; mSolve`PowerSolver[mRules0_, myMList_, opts___?OptionQ]:= Module[{powerSolverDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), mRules = Sort[mRules0, (Length[#1] < Length[#2])&], numberOfEquations = Length[myMList] }, (* Protected local variables. *) If[powerSolverDebug, Print["The mRules used in Power Solver (sorted from " <> "shortest to longest) is:"]; Print[CleanUpFunction[mRules]]; ]; (* A little warning message if there is a combinatorial explosion. *) If[Binomial[Length[Flatten[ mRules ] ], numberOfEquations] > 100, CellPrint[ Cell["Formulating and solving up to " <> ToString[ Binomial[Length[Flatten[ mRules ] ], numberOfEquations] ] <> " systems of equations. Please be patient.", "Message" ] ]; ]; (* Forms a list of equations to be solved for m_i. *) eqnList = Outer[List, Sequence @@ (mRules /. Rule->Equal)]; (* Since Outer creates a list numberOfEquations deep, we *) (* drop it down into the correct form. DB:5/20/2003 *) eqnList = Partition[ Flatten[ eqnList ], numberOfEquations ]; If[powerSolverDebug, Print["The next set of equations to be solved are:"]; Print[CleanUpFunction[eqnList]]; ]; (* Takes these solutions and uses them to find *) (* actual integer solutions to $m_i$ *) mSoln = Union[ Sequence @@ Solve[#, myMList]& /@ eqnList ]; If[powerSolverDebug, Print["The solutions are:"]; Print[CleanUpFunction[mSoln]]; ]; (* Returns the solutions *) Return[mSoln]; ]; mSolve`SystemCleanUp[eqnList0_, mList_, mSoln0_, myMList_, F_, opts___?OptionQ]:= Module[{systemCleanUpDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), mSoln = Union[Sort /@ mSoln0], (* A local version of mSoln0. *) chiList = {}, (* The list of chi constraints. *) eqnList, (* A local copy of the equations. *) numberOfEquations = Length[eqnList0] (* The number of equations. *) }, (* Protected local variables. *) (* Applies the solutions to the expressions of m_i. *) mSoln = Transpose[{(mList //. #)& /@ mSoln, mSoln}]; If[systemCleanUpDebug, Print["The expressions of m_i after possible " <> "solution substitution:"]; Print[CleanUpFunction[mSoln]]; ]; mSoln = mSoln /. {a_List, {(_Rule)..}} :> Sequence[] /; Not[ And @@ (FreeQ[ a, #] & /@ myMList ) ]; If[systemCleanUpDebug, Print["mSoln, after the underdetermined systems are removed:"]; Print[CleanUpFunction[mSoln]]; ]; (* Removes solutions that do not balance at least two independent *) (* Terms in the original system. *) mSoln = mSoln /. { {a_List, b : {(_Rule) ..}} :> b /; And @@ ((Length[ Cases[#, Max[#]] ] >= 2)& /@ a), {a_List, b : {(_Rule) ..}} :> Sequence[] }; If[systemCleanUpDebug, Print["mSoln, after power balances are checked:"]; Print[CleanUpFunction[mSoln]]; ]; (* Removes solutions with any of the $M_i = 0$. DB:5/25/2003 *) mSoln = mSoln /. a : {(_Rule) ..} :> Sequence[] /; Or @@ ( #[[2]] === 0 & /@ a ); If[systemCleanUpDebug, Print["mSoln, after zero solutions are removed:"]; Print[CleanUpFunction[ mSoln ]]; ]; (* If we are dealing with a system of equations, than we remove negative *) (* and rational solutions for the dominate powers. DB:5/25/2003 *) mSoln = mSoln //. { a : {(_Rule) ..} :> Sequence[] /; Or @@ (! IntegerQ[ #[[2]] ]& /@ a ), a : {(_Rule) ..} :> Sequence[] /; Or @@ ( #[[2]] < 0 & /@ a ) }; If[systemCleanUpDebug, Print[ "If it is a system, we remove negative and non-integer solutions:" ]; Print[CleanUpFunction[ mSoln ]]; ]; (* Removes Tracking Variables. *) eqnList = eqnList0 /. myTrackingVariable[_]:>1; If[systemCleanUpDebug, Print["eqnList, sans tracking variables is:"]; Print[CleanUpFunction[eqnList]]; ]; (* Forms chiList. *) eqnList /. \[Chi][i_] :> ( chiList = Append[chiList, \[Chi][i] ]; 1 ); chiList = Union[chiList]; If[systemCleanUpDebug, Print["chiList, the list of chi s are:"]; Print[CleanUpFunction[chiList]]; ]; (* Finds constraints on chi. *) eqnList = MapAll[Factor, Plus @@ Coefficient[(eqnList /. #), F^Max[mList /. mSoln] ]& /@ mSoln ]; If[systemCleanUpDebug, Print["mSoln, Checking in eqnList:"]; Print[CleanUpFunction[eqnList]]; ]; (* Returns the good solutions. *) Return[mSoln]; ]; (* Taken from Painleve Test code on 5/15/2004. *) mSolve`FixFreeM[ mSoln_List, myMList_List, opts___?OptionQ ] := Module[{ fixFreeMDebug = (* Debugging bool flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), mSoln0, mFree, mFreeValues, mFixedValues, mValues }, (* Protected Local Variables *) (* Remove symbolic solutions, based on the code in SystemCleanUp. *) (* DB:7/29/2004 added removal of rational solutions. *) mSoln0 = ( {Rest /@ (# /. Rule -> List), #}& /@ mSoln ) /. {a_List, {(_Rule)..}} :> Sequence[] /; Not[ And @@ (FreeQ[ a, #] & /@ myMList ) ]; mSoln0 = mSoln0 /. {a_List, {(_Rule)..}} :> Sequence[] /; Not[ And @@ FreeQ[ a, Rational] ]; mSoln0 = #[[2]]& /@ mSoln0; (* Warn the user when potential solutions are removed. *) If[Length[Complement[Union[Sort /@ mSoln], mSoln0] ] > 0, StylePrint[ "The potential solutions "<> ToString[ InputForm[ CleanUpFunction[ Complement[Union[Sort /@ mSoln], mSoln0] ] ] ] <> " are being removed because they are under determined. ", "Message" ]; ]; If[fixFreeMDebug, Print["mSoln, after the underdetermined systems are removed:"]; Print[CleanUpFunction[mSoln0]]; ]; (* Pick out the solutions with freedom. *) mFree = Select[mSoln0, Length[#] < Length[myMList]&]; If[fixFreeMDebug, Print["The dominant behaviors with one or more free \!\(M\_i\):"]; Print[CleanUpFunction[ mFree ]]; ]; If[Length[mFree] >= 1, Message[PDESpecialSolutions`dominantBehavior::freevalues, CleanUpFunction[ mFree ] ]; (* Substitutes in the values that are fixed, leaving the free values. *) mFreeValues = (myMList /. #)& /@ mFree; mFixedValues = (myMList /. #)& /@ Complement[mSoln0, mFree]; (* DB:2/8/2004 *) If[Length[mFixedValues] === 0, Message[PDESpecialSolutions`dominantBehavior::underdetermined]; mFixedValues = mFreeValues /. m[_] :> 0 ]; If[fixFreeMDebug, Print["The free m values are:"]; Print[CleanUpFunction[ mFreeValues ]]; Print["The fixed m values are:"]; Print[CleanUpFunction[ mFixedValues ]]; ]; mValues = Transpose[ { #, Sequence @@ Cases[ mFixedValues, # /. (m[i_] :> _) ] } ]& /@ mFreeValues; If[fixFreeMDebug, Print["The fixed values divided by free values:"]; Print[CleanUpFunction[ mValues ]]; ]; mValues = Sequence @@ Thread[ If[Head[First[#]] === Integer, First[#], Range[ MinDegreeOfThePolynomial /. {opts} /. Options[PDESpecialSolutions], Max[Rest[#]] ] ]& /@ # ]& /@ mValues; If[fixFreeMDebug, Print["Fixing the free values from the " <> "MinDegreeOfThePolynomial upto max:"]; Print[CleanUpFunction[ mValues ]]; ]; mValues = ( Rule @@ #& /@ Transpose[ { myMList, # } ] )& /@ mValues; If[fixFreeMDebug, Print["The fixed values reformated:"]; Print[CleanUpFunction[ mValues ]]; ]; Return[ mValues ]; ]; Return[mSoln0]; ]; (* Finds free solutions that result from inequalities. *) (* Added direct from old code on 11/20/2002 by DB. *) mSolve`GenerateAlternativeSolutions[ mList_List, opts___? OptionQ ]:= Module[{ n (* The number of alpha[i] *) }, (* Protected Local Variables *) n = Length[ mList ]; Return[ Thread[ mList -> # ] & /@ Flatten[ Outer[ List, Sequence @@ Table[ Range[ ( MinDegreeOfThePolynomial /. {opts} /. Options[PDESpecialSolutions] ), ( MaxDegreeOfThePolynomial /. {opts} /. Options[PDESpecialSolutions] ) ], {n} ] ], n - 1 ] ] ]; (* ************************************************************************ *) (* : Title : PDESpecialSolutionsBuildSystem. * : Author : Douglas Baldwin * : Date : 05-24-01 *) (* : Update : 1/12/2003 by Douglas Baldwin * Added debug statements and highest order equation solving. * Changes marked with DB:Date *) (* : Update : 1-21-2003 by Douglas Baldwin * Added StripNonZero to the equations before they are sent to the solver. * This is mainly for printing the simplest system to the user, because * AnalyzeAndSolve already strips the non-zeros. *) PDESpecialSolutionsBuildSystem[mSoln_List, system_List, vars_List, parameters_List, opts___?OptionQ]:= Block[{buildSystemDebug = False, Form, numberOfEquations, sqrtRules, funcRules, uRules, tempuRules, newSystem, waveparameters, param, unknowns, nonzeros, maxF, highestOrderSystem, highestOrderSolutions, myTemp, time = TimeUsed[]}, (* Protects Local Variables. *) (* Sets the form from the options. *) Form = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn->sn, jacobicn->cn}; (* Sets a local variable to the number of equations given by user. *) numberOfEquations = Length[system]; (* Sets up the way in which the system will be divided up *) (* in the next step, based on the Options set by the user. *) sqrtRules = {tanh -> 0, sech -> 0, sechtanh -> 1-S^2, cn -> 0, sn -> 0 }; (* Sets up a list of rules for replacing the form to F. *) funcRules = {tanh -> T, sech -> S, sechtanh -> S, cn -> CN, sn -> SN}; (* Converts the results of mSolve to an expression in T *) uRules = (# /. (m[i_]->j_) :> (u[i][Form /. funcRules] -> Evaluate[ Sum[a[i,kk]*(Form /. funcRules)^kk, {kk, 0, Numerator[j]} ]^(Sign[j] / Denominator[j]) + If[(Form /. sqrtRules) =!= 0, Sqrt[Form /. sqrtRules]* Sum[b[i,kk]*(Form /. funcRules)^kk, {kk, 0, n[i] /. Flatten[#]} ], 0 ] ] ) )& /@ mSoln; (* WILLY 12/19/2001 modification: introduced tempuRules *) If[Verbose /. {opts} /. Options[PDESpecialSolutions], Print["Seeking polynomial solutions of the form"]; tempuRules = DeleteCases[uRules, n[_]->_, 2]; Print[CleanUpFunction[tempuRules]]; ]; (* Prints out a warning, if it's taking a long time. *) If[TimeUsed[]-time>3, Print["Still building the nonlinear algebraic system, please wait."] ]; (* Converts form of solutions given by solver to a pure *) (* function. *) toPure = (# /. (a__[var__]->temp__):> (a:>Function[{var}, temp]))&; (* Applies the expressions in T to the system and removes the *) (* tracking variables. *) newSystem = Expand[ system //. Append[toPure[#], myTrackingVariable[_]->1] ]& /@ uRules; If[buildSystemDebug, Print["Apply the expressions in T to the sytem gives: "]; Print[CleanUpFunction[ newSystem ]]; ]; (* Clears the denominator to simplify the system. *) newSystem = Map[Together[#*Denominator[Together[#]]]&, newSystem, 2]; If[buildSystemDebug, Print["Clearing the denominator, we get: "]; Print[CleanUpFunction[ newSystem ]]; ]; (* Creates a list with the wave parameters (c[1], c[2], ...) to be *) (* passed to the solver. DB:2-4-2003 Put mod with param instead of c_is.*) waveparameters = Table[c[i], {i, Length[vars]}]; (* DB:2-4-2003 Including mod with parameters. *) param = Join[parameters, If[Form === cn || Form === sn, {mod}, {}] ]; (* Creates a list of all the a[i,j]'s to be passed to the solver. *) unknowns = Table[ {Table[a[i,j], {j, 0, Numerator[m[i] /. #]}], If[(n[i] /. # /. n[_]->0) > 0, Table[b[i,j], {j, 0, n[i] /. # /. n[_]->0}], Null ]}, {i, numberOfEquations}]& /@ mSoln; unknowns = DeleteCases[Sequence @@ #& /@ #& /@ unknowns, Null, Infinity]; (* Creates a list of those variables which must be nonzero, so as *) (* to simplify the work of the solver. *) (* DB:1/12/2003 Commented out the Form =!= cn and sn. *) nonzeros = Join[param, If[Form =!= sechtanh (* && Form =!= cn && Form =!= sn *), Last /@ #, {} ], waveparameters ]& /@ unknowns; (* Flattens the sublists and reorders them to speed up the solver. *) unknowns = Map[ Reverse, #[[2]]& /@ #& /@ (Sort[{# /. a_[b_,c_]:>{c,a,b},#}& /@ Flatten[#] ]& /@ unknowns) ]; If[buildSystemDebug, Print["The unknowns, nonzeros and wave parameters are: "]; Print[CleanUpFunction[{unknowns, nonzeros, waveparameters} ]]; ]; (* Expands the system. *) newSystem = MapAll[Expand,newSystem]; (* Determines the maximum value of F in each of the cases. *) maxF = Map[ Exponent[#, Form /. funcRules ]&, newSystem ] /. -Infinity -> 0 (* DB:3/21/2003 Infinity rule added. *); If[buildSystemDebug, Print["The dominate power in each expression of new system is: "]; Print[CleanUpFunction[ maxF ]]; ]; (* UnsortedUnion from Mathematica Book -- added by DB on 1/11/2003. *) (* Removed on Tuesday, January 14, 2003 by DB UnsortedUnion[x_] := Module[{f}, f[y_] := (f[y] = Sequence[]; y); f /@ x]; *) (* If the HigestOrderFirst option is set, then: *) If[HighestOrderFirst /. {opts} /. Options[PDESpecialSolutions], (* Creates a list of the highest order system. *) (* Rewritten by DB:1/14/2003 *) highestOrderSystem = (# == 0)& /@ MapThread[Coefficient[#1, Form /. funcRules, #2] &, {#, Exponent[#, Form /. funcRules]}]& /@ newSystem; If[buildSystemDebug, Print["The highest order equations are: "]; Print[CleanUpFunction[ highestOrderSystem ]]; ]; (* DB:1-21-2003, Strips the non-zeros at this step. *) highestOrderSystem = MapThread[StripNonZero, {highestOrderSystem, nonzeros} ]; If[buildSystemDebug, Print["Removes non-zeros from the system."]; Print[CleanUpFunction[ highestOrderSystem ]]; ]; (* DB:1/12/2003 highestOrderSytem is then formated to *) (* input into AnalyzeAndSolve *) highestOrderSystem = MapAll[Factor, Transpose[ {highestOrderSystem, unknowns, Table[waveparameters, {Length[newSystem]}], Table[param, {Length[newSystem]}], nonzeros } ] ]; If[buildSystemDebug, Print["The highest order system formated for AnalyzeAndSolve: "]; Print[CleanUpFunction[ highestOrderSystem ]]; ]; (* DB:1/12/2003 Find the solutions using AnalyzeAndSolve *) highestOrderSolutions = Algebra`AnalyzeAndSolve`AnalyzeAndSolve @@ #& /@ highestOrderSystem; If[buildSystemDebug, Print["The highest order equations' solutions are: "]; Print[CleanUpFunction[ highestOrderSolutions ]]; ]; ]; (* Sets up the way in which the system will be divided up *) (* in the next step, based on the Options set by the user. *) sqrtRules = {tanh -> Null, sech -> 1-S^2, sechtanh -> 1-S^2, cn -> 1-mod-CN^2+2*mod*CN^2-mod*CN^4, sn -> 1-SN^2-mod*SN^2+mod*SN^4 }; (* DB:2/28/2005 Added new Coefficient function, for 5.0 compatibility. *) myCoefficient[p_Plus, q_, 0] := Plus @@ Select[p, FreeQ[#,q]&]; myCoefficient[p_, q_, 0] := Plus @@ Select[{p}, FreeQ[#,q]&]; myCoefficient[p_Plus, q_, 1] := Plus @@ Cases[p, z_. q -> z]; myCoefficient[p_, q_, 1] := Plus @@ Cases[{p}, z_. q -> z]; (* Breaks up system into most general form. *) newSystem = Flatten[ {myCoefficient[#,Sqrt[Form /. sqrtRules],0]& /@ #, myCoefficient[#,Sqrt[Form /. sqrtRules],1]& /@ #} ]& /@ #& /@ newSystem; If[buildSystemDebug, Print["After breaking the system up into {P,Q} (where P + Sqrt[]*Q): "]; Print[CleanUpFunction[ newSystem ]]; ]; (* Brakes the expressions up newSystem by powers of T. *) (* Modified version by DB:1/11/2003 *) newSystem = Table[ Table[ Union[ DeleteCases[ Flatten[ Table[ Table[ Coefficient[ Expand[newSystem[[case,k]] ], Form /. funcRules, i ], {i, 0, maxF[[case,j]]}], {j, Length[maxF[[case]] ]} ] ], 0 ] ], {k, Length[newSystem[[case]] ]} ], {case, Length[newSystem]} ]; If[buildSystemDebug, Print["The system after being broken up is: "]; Print[CleanUpFunction[ newSystem ]]; ]; (* Converts from expressions to equations. *) newSystem = DeleteCases[(# == 0)& /@ Flatten[#]& /@ newSystem, False || True, Infinity]; If[buildSystemDebug, Print["After flattening and removing True/Falses:"]; Print[CleanUpFunction[ newSystem ]]; ]; (* DB:1-21-2003, Strips the non-zeros at this step. *) newSystem = MapThread[StripNonZero, {newSystem, nonzeros} ]; If[buildSystemDebug, Print["Removes non-zeros from the system."]; Print[CleanUpFunction[ newSystem ]]; ]; (* Prints out a warning if it's taking too long. *) If[ (TimeUsed[]-time>3) && !(Verbose /. {opts} /. Options[PDESpecialSolutions]), Print["Still building the nonlinear algebraic system, please wait."] ]; If[HighestOrderFirst /. {opts} /. Options[PDESpecialSolutions], (* DB:11/13/2002 Version to work with highest order solutions *) (* Combines into lists containing {newSystem, unknowns, *) (* waveparameters, param, nonzeros} for the solver. *) newSystem = MapAll[Factor, (* Expands the cases over the solutions for the highest terms *) (Sequence @@ ( { (* Added Flattens on 11/13/02 to compensate *) (* for the output of highestOrderSolutions *) (* StripNonZero added on 2/2/2003 by DB. *) StripNonZero[ Flatten[(First[#] /. Last[#]) /. True -> Sequence[] ], #[[-2]] ], Sequence @@ Drop[Rest[#],-1], Flatten[Last[#]] }& /@ Thread[#, List, -1] ) ) & /@ (* Groups the information in the same form as the input *) (* to the solver. *) Transpose[ {newSystem, unknowns, Table[waveparameters, {Length[newSystem]}], Table[param, {Length[newSystem]}], nonzeros, highestOrderSolutions } ] ], (* Combines into lists containing {newSystem, unknowns, *) (* waveparameters, param, nonzeros} for the solver. *) newSystem = MapAll[Factor, Transpose[ {newSystem, unknowns, Table[waveparameters, {Length[newSystem]}], Table[param, {Length[newSystem]}], nonzeros } ] ] ]; If[buildSystemDebug, Print["The system formatted for AnalyzeAndSolve: "]; Print[CleanUpFunction[ newSystem ]]; ]; (* DB:3/21/2003 Allows mSoln to be updated in main. *) newSystem = Thread[{mSoln, newSystem}]; If[buildSystemDebug, Print["Threading mSoln"]; Print[CleanUpFunction[ newSystem]]; ]; (* DB:3/21/2003 Remove inconsistant systems. *) newSystem = newSystem /. { {_,{{___, False, ___}, __}} :> Sequence[], {_,{{}, __}} :> Sequence[] }; If[buildSystemDebug, Print["Removes inconsistant systems: "]; Print[CleanUpFunction[ newSystem ]]; ]; (* If there are no systems to be solved, then abort. *) (* DB:3/21/2003 *) If[Length[newSystem] === 0, Message[ PDESpecialSolutionsBuildSystem::fail ]; Abort[ ] ]; (* DB:3/21/2003 Transpose added after threading mSoln above. *) Return[ Transpose[newSystem] ]; ]; (* ************************************************************************ *) (* : Title : PDESpecialSolutionsBuildSolutions. *) (* : Author : Douglas Baldwin *) (* : Summary : Builds (includes massive simplification) *) (* and test the final solutions. *) PDESpecialSolutionsBuildSolutions[solution_, mSoln_, eqns_, funcs_, vars_, param_, opts___?OptionQ] := Block[{NumericDebug = False, (* Debug for numeric tester. *) Form, (* Form defined by user. *) solutionRules, (* Local version of solution *) sqrtRules, (* Set by options. *) funcRules, (* Also set by options. *) uRules, (* Rules for creating the solutions. *) mySimplificationRules, (* Simplification rules *) numericTestSolutions, (* tests validity of solutions *) symbolicTestSolutions, (* Also, tests the validity of soln. *) solns = {}, (* confirmed solutions. *) warn = {}, (* Collects a list of warnings due to simplification. *) time = TimeUsed[] (* A timer for length warnings. *) }, (* Protected Local Variables. *) (* Sets the form from the options. *) Form = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn -> sn, jacobicn -> cn}; (* Sets up a local version to modify *) solutionRules = MapAll[Factor, solution]; (* Rules for building the solutions. *) sqrtRules = {tanh -> 0, sech -> 0, sechtanh -> 1-Sech[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}] + phase]^2, cn -> 0, sn -> 0 }; funcRules = {tanh -> Tanh[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase], sech -> Sech[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase], sechtanh -> Sech[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase], cn -> CN[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase,mod], sn -> SN[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase,mod] }; (* Converts the results of mSolve to an expression in T *) uRules = (m[i_]->j_) :> (funcs[[i]] -> Sum[a[i,kk]*(Form /. funcRules)^kk, {kk, 0, Numerator[j]} ]^(1/Denominator[j]) + Sqrt[Form /. sqrtRules]* Sum[b[i,kk]*(Form /. funcRules)^kk, {kk, 0, n[i] /. Flatten[#]} ] )&; (* Builds the solutions from the above rules and *) (* the powers of F listed in the passed mSoln. *) solutionRules = DeleteCases[ DeleteCases[ Table[ Table[ Join[(mSoln[[jj]] /. uRules[mSoln[[jj]] ]) //. solutionRules[[jj, ii]], (# -> (# /. solutionRules[[jj,ii]]))& /@ If[Form === cn || Form === sn, Append[param, mod], param ] ], {ii, Length[solutionRules[[jj]]]} ], {jj, Length[mSoln]} ], n[_]->_, Infinity ], Rule[a_, b_] /; a == b && FreeQ[a, Alternatives @@ funcs], Infinity ]; (* Simplification rules for reducing the solution *) (* even further . . . so as to find more general *) (* solutions. *) mySimplificationRules = { Sqrt[xxx__^2] :> (warn = Append[warn, "Sqrt[a^2]->a"]; xxx), Sqrt[Power[yyy__,2]] :> (warn = Append[warn, "Sqrt[a^2]->a"]; yyy), Sqrt[-zzz__^2] :> (warn = Append[warn, "Sqrt[-a^2]->I*a"]; I*zzz), Sqrt[-ttt_*zzz__^2] :> (warn = Append[warn, "Sqrt[-a*b^2]->I*b*Sqrt[a]"]; I*zzz*Sqrt[ttt]), (1-Sech[xx__]^2) -> Tanh[xx]^2, (1-Tanh[xx__]^2) -> Sech[xx]^2, JacobiDN[xx__, mm__]^2 :> 1-mm+mm*JacobiCN[xx,mm]^2 /; Form === cn, JacobiSN[xx__, mm__]^2 :> 1-JacobiCN[xx,mm]^2 /; Form === cn, JacobiDN[xx__, mm__] :> (warn = Append[warn,"JacobiDN[x,mod]->Sqrt[1-mod+mod*JacobiCN[x,mod]^2]"]; Sqrt[1-mm+mm*JacobiCN[xx,mm]^2] ) /; Form === cn, JacobiSN[xx__, mm__] :> (warn = Append[warn,"JacobiSN[x,mod]->Sqrt[1-JacobiCN[x,mod]^2]"]; Sqrt[1-JacobiCN[xx,mm]^2] ) /; Form === cn, JacobiDN[xx__, mm__]^2 :>1-mm*JacobiSN[xx,mm]^2 /; Form === sn, JacobiCN[xx__, mm__]^2 :> 1-JacobiSN[xx,mm]^2 /; Form === sn, JacobiDN[xx__, mm__] :> (warn = Append[warn, "JacobiDN[x,mod]->Sqrt[1-mod*JacobiSN[x,mod]^2]" ]; Sqrt[1-mm*JacobiSN[xx,mm]^2] ) /; Form === sn, JacobiCN[xx__, mm__] :> (warn = Append[warn, "JacobiSN[x,mod]->Sqrt[1-JacobiSN[x,mod]^2]]" ]; Sqrt[1-JacobiSN[xx,mm]^2] ) /; Form === sn }; (* Applies the above rules to the solutions. *) solutionRules = FixedPoint[ MapAll[Expand, MapAll[Factor, # //. mySimplificationRules ] //. mySimplificationRules ]&, solutionRules, 3 ]; (* DB:3/21/2003 Removes zero solutions. *) solutionRules = solutionRules /. {___,_-> 0,___} :> Sequence[]; (* Removes constant solutions *) solutionRules = Select[#, !(And @@ (FreeQ[# /. (a_->b__):>b, Alternatives @@ vars]& /@ #))& ]& /@ solutionRules; If[Verbose /. {opts} /. Options[PDESpecialSolutions], Print["The possible non-trivial solutions (before any testing) are: "]; Print[ CleanUpFunction[ MapAll[Factor, solutionRules /. {CN -> JacobiCN, SN -> JacobiSN} ] ] ]; ]; (* Prints out a warning, if it's taking a long time. *) If[TimeUsed[]-time>3, Print["Still building the exact solutions, please wait."] ]; (* Depending on the NumericTest option, it either tests the *) (* solutions numerically, or it doesn't. *) If[ ( !(NumericTest /. {opts} /. Options[PDESpecialSolutions]) && !(SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), CellPrint[ Cell[ "These solutions are not being tested numerically or symbolically. " <> "To test the solutions set either the NumericTest option " <> "to True, or set the SymbolicTest option to True, or both. ", "Message" ] ]; If[Length[warn]>0, CellPrint[ Cell[ "The following simplification rules are being used: " <> ToString[Union[warn]], "Message" ] ] ]; Return[ Map[Union, MapAll[Factor, MapAll[Expand, solutionRules /. {CN -> JacobiCN, SN -> JacobiSN} ] ], 2 ] ] ]; (* Converts form of solutions given by solver to a pure function. *) toPure = (# /. (a__[var__]->temp__):> (a:>Function[{var}, temp]))&; (* If statement for the numeric testing option. *) If[NumericTest /. {opts} /. Options[PDESpecialSolutions], (* Sends message to user. *) If[!(VerboseTest /. {opts} /. Options[DDESolve]), Print["Numerically testing the solutions."] ]; (* Tests numerically to make sure the above solutions are valid. *) numericTestSolutions = {(eqns /. a__==b__:>a-b) /. toPure[MapAll[TrigToExp, #]], #}& /@ #& /@ solutionRules; randomVarRules = Append[ (# -> Random[Real, {-1, 1}])& /@ Join[vars, Array[c, Length[vars]]], phase -> 0 ]; randomAijRules = { a[_,_]->Random[Real, {-1, 1}], b[_,_]->Random[Real, {-1, 1}], Sequence @@ ((# -> Random[Real, {-1, 1}])& /@ If[Form === cn || Form === sn, Join[param, {mod}], param ] ) }; (* Numerically tests the solutions by replacing all the symbols with *) (* random numbers. *) numericTestSolutions = { Or[ And @@ ( Abs[ # ] < 10^(-$MachinePrecision/2)& /@ ( MapAll[ Expand, #[[1]] /. randomVarRules ] /. randomAijRules /. {CN -> JacobiCN, SN-> JacobiSN} ) ), And @@ ( Abs[ # ] < 10^(-$MachinePrecision/2)& /@ ( MapAll[ Expand, #[[1]] /. ( randomVarRules /. (a_ -> b_) :> (a -> -b) ) ] /. ( randomAijRules /. (a_ -> b_) :> (a -> -b) ) /. {CN -> JacobiCN, SN-> JacobiSN} ) ) ], # }& /@ #& /@ numericTestSolutions; (* If it works for any of the trials, it keeps the solution. *) solns = numericTestSolutions /. { {True, a_} :> a[[2]], {False, _} :> Sequence[] }; numericTestSolutions = Complement[solutionRules, solns]; If[(Map[Union, numericTestSolutions, 3] //. {{}}->{}) != {}, CellPrint[ Cell[ "These solutions did not simplify to be less than 10^(" <> ToString[-$MachinePrecision/2] <> ") when random numbers in [-1,1] were substituted for the " <> "unknowns. " <> "This could be caused from Mathematica's assumption that " <> "Sqrt[a^2] -> a when a is numeric. " <> "Please test these equations by hand or interactively "<> "with Mathematica.", "Message" ] ]; Print[CleanUpFunction[#]]& /@ numericTestSolutions; ]; ]; (* End of numeric test if statement. *) (* The second testing If statement. *) If[SymbolicTest /. {opts} /. Options[PDESpecialSolutions], (* Sends message to user. *) If[!(VerboseTest /. {opts} /. Options[DDESolve]), Print["Symbolically testing the solutions."] ]; (* This sets up the input for the next block of code. *) (* In that, it takes the equestion given by the user, *) (* and replaces the user defined functions (e.g. u[x,t]) *) (* with the functions found with the algorithm. To replace *) (* the partial derivatives correctly, we must first convert *) (* the solutions to pure functions. *) symbolicTestSolutions = {(eqns /. a__==b__:>a-b) /. toPure[MapAll[TrigToExp, #]], #}& /@ #& /@ (solutionRules /. {CN -> JacobiCN, SN-> JacobiSN}); (* Attempts to simplify the expressions as much as possible. *) symbolicTestSolutions = ExpToTrig[ FixedPoint[ MapAll[Factor, MapAll[Expand, #] /. mySimplificationRules ] /. mySimplificationRules &, #, 10]& /@ symbolicTestSolutions ]; (* Pulls of the solutions which simplify to zero, and adds *) (* them to solns (the final output applies the union, so *) (* duplicate solutions are not an issue at this point). *) solns = Join[solns, Cases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ] //. {{(0)..}, a_List}:>{a} ]; (* Removes all the good cases, so we may output the bad *) (* ones to the user. *) symbolicTestSolutions = DeleteCases[ DeleteCases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ], {} ]; If[(Map[Union, symbolicTestSolutions,3] //. {{}}->{}) != {}, (* Sends out the left over (questionable) solutions to *) (* the user for hand inspection. *) CellPrint[ Cell[ "These solutions did not simplify to zero. " <> "Please test these equations by hand or interactively "<> "with Mathematica.", "Message" ] ]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; ]; ]; If[Length[warn]>0, CellPrint[ Cell[ "The following simplification rules are being used " <> "to test the solutions: " <> ToString[Union[warn]], "Message" ] ] ]; (* Returns the solution rules *) Return[ Union[ MapAll[Factor, MapAll[Expand, solns /. {CN -> JacobiCN, SN-> JacobiSN} ] ] ] ]; ]; (* ************************************************************************ *) (* : Title : CleanUpFunction *) (* : Author : Douglas Baldwin *) (* : Summary : Remove "PDESpecialSolutionsPrivate`" from output. *) CleanUpFunction = ToExpression[ StringReplace[ToString[InputForm[#]], "Calculus`PDESpecialSolutions`Private`"->"" ] ]& (* ************************************************************************ *) (* : Title : StripNonZero *) (* : Author : Douglas Baldwin *) (* : Summary : Strips parameters from expressions. *) StripNonZero[{a : (_List) ..}, param_List] := StripNonZero[#, param] & /@ {a}; StripNonZero[theList_List, param_List] := Module[ {result, stripDebug = False}, (* Local Variables *) If[stripDebug, Print["The origional expression is: "]; Print[CleanUpFunction[ MapAll[Factor, theList] ]]; ]; (* Maps factor to every term, so as to have a constant base. *) result = FactorList /@ (theList /. Equal[a_,b_]:>a-b); If[stripDebug, Print["The results of FactorList is: "]; Print[CleanUpFunction[ result ]]; ]; If[stripDebug, Print["The rules are: "]; Print[ CleanUpFunction[ {({#^_:1, _} :> Sequence[] & /@ param), {_?NumericQ, _} :> Sequence[] } ] ]; ]; (* Remove terms that are non-zero. *) result = result /. ({#^_:1, _} :> Sequence[] & /@ param); (* Remove terms that are numeric. *) result = result /. {_?NumericQ, _} :> Sequence[]; If[stripDebug, Print["The result after apply the non-zero rules: "]; Print[CleanUpFunction[ result ]]; ]; (* Puts it back into standard form, {a, b} to a*b *) result = Times @@ (First[#]^Last[#]&) /@ #& /@ result; If[stripDebug, Print["Converting back to the standard form: "]; Print[CleanUpFunction[ result ]]; ]; (* Converts expressions back into equations. *) If[! FreeQ[theList, Equal], result = Equal[#,0]& /@ result; ]; If[stripDebug, Print["The final result of stripping the non-zeros is: "]; Print[CleanUpFunction[ result ]]; ]; Return[ result ]; ]; End[] (* `Private` context *) Protect[PDESpecialSolutions] SetAttributes[PDESpecialSolutions, ReadProtected] EndPackage[] (* ************************************************************************ *) (* Nonlinear algebraic solver by Unal Goktas (WRI) and Willy Hereman *) (* written by Unal Goktas in October/November 2000 *) (* Comments added by Douglas Baldwin on Monday, January 13, 2003 *) BeginPackage["Algebra`AnalyzeAndSolve`"] Unprotect[AnalyzeAndSolve] AnalyzeAndSolve::soln = "The solultion `1` is being removed because it is inconsistent." Begin["`Private`"] If[$VersionNumber < 4, Attributes[Internal`DeactivateMessages] = {HoldAll}; Internal`DeactivateMessages[body_, msgs___MessageName] := Module[{wasOn = Select[Hold[msgs], Head[#] =!= $Off &], result}, CheckAbort[ Scan[Off, wasOn]; result = body; Scan[On, wasOn]; result, Scan[On, wasOn]; Abort[] ] ] ] If[$VersionNumber >= 5, Developer`SetSystemOptions[ "ReduceOptions" -> {"ReorderVariables" -> False} ]; myReduce[eqn_, vars_] := LogicalExpand[ Reduce[eqn, Reverse[vars]] /. Equal -> myEqual ] /. myEqual -> Equal, myReduce[eqn_, vars_] := Reduce[eqn, vars, Sort -> False] ] (* : Title : CleanUpFunction *) (* : Author : Douglas Baldwin *) (* : Summary : Remove "Algebra`AnalyzeAndSolve`Private`" from output. *) (* : Added for debugging for AnalyzeAndSolve by DB:1/14/2003 *) CleanUpFunction = ToExpression[ StringReplace[ToString[InputForm[#]], { "Calculus`PDESpecialSolutions`Private`"->"", "Algebra`AnalyzeAndSolve`Private`"->""} ] ]& (* This is the first call of AnalyzeAndSolve. In this call, *) (* The equations are converted into expressions equal to zero. *) AnalyzeAndSolve[system: {HoldPattern[_ == _]..}, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_] := AnalyzeAndSolve[(#[[1]]-#[[2]]&) /@ system, primaryunknowns, secondaryunknowns, parameters, nonzeros] (* ************************************************************************ *) (* Takes the newly recast version of AnalyzeAndSolve and converts its inputs *) (* into a form readable by RecursivelyKeepTrack. *) (* Warning: It is assumed that "system" is polynomial in "primaryunknowns", *) (* "secondaryunknowns", and "parameters". *) AnalyzeAndSolve[system_, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_] := Block[{a, globalsol = {}, constraints}, a = Together[({#}& /@ system) /. {{}}]; constraints = (And @@ ((# != 0 &) /@ nonzeros)); MapThread[ RecursivelyKeepTrack[#1, #2, primaryunknowns, secondaryunknowns, parameters, nonzeros, constraints]&, {a, {{}}}]; Union[globalsol] ] (* ************************************************************************ *) (* Collects all the solutions at the end into globalsol. *) (* Origional code . . . replaced by the two following functions. *) RecursivelyKeepTrack[{}, sol_, __] := (AppendTo[globalsol, sol]; {}) (* Terminates that recursive branch if the solution expodes or indeterminate.*) RecursivelyKeepTrack[system_, sol_, __] /; (!FreeQ[system, DirectedInfinity] || !FreeQ[system, Indeterminate]) := {} (* Main RecursivelyKeepTrack call. It takes the system, cleans it up, *) (* sorts it by complexity level (heuristic based on the degree of the *) (* polynomial and size of the expression), takes first equation solves *) (* it, simplifies the solution and recursively calls it self on each *) (* unique branch of the solution. *) RecursivelyKeepTrack[system_, sol_, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_, constraints_] := Block[{a, b, c, recursiveDebug = False}, (* Breaks the expression into a list and removes nonzeros and numbers. *) a = FactorListAndCleanUp[#, primaryunknowns, secondaryunknowns, parameters, nonzeros]& /@ system; (* Removes duplicates. Message added DB:5/23/2003 *) a = Union[a] /. {} :> (Message[AnalyzeAndSolve::soln, CleanUpFunction[sol]];{}); (* Sorts the system by the "ComplexityLevel" heuristic which is the *) (* degree of the polynomial*100 + LeafCount size of the expression. *) a = Sort[a, (ComplexityLevel[#1, primaryunknowns, secondaryunknowns, parameters] <= ComplexityLevel[#2, primaryunknowns, secondaryunknowns, parameters]&)]; (* DB:1/14/2003 Print statement added for debugging. *) If[recursiveDebug, Print["The order of the equations in AnalyzeAndSolve is: "]; Print[ CleanUpFunction[ {ComplexityLevel[#1, primaryunknowns, secondaryunknowns, parameters],#1}& /@ a ] ]; ]; (* Takes the first equation of the system of equations (which is *) (* the simplest according to the complexity heuristic) *) b = a[[1]]; (* Then solves the first equation for only one variable. *) b = SolveStepByStep[#, primaryunknowns, secondaryunknowns, parameters, constraints, sol]& /@ b; (* Flattens from the inside out: {{{a,b},{c,d}}} -> {{a,b},{c,d}}. *) b = (Sequence @@ # &) /@ b; (* DB:2/4/2003 Added expand. *) b = MapAll[Factor, b]; If[recursiveDebug, Print["The solution from the first equation is:"]; Print[CleanUpFunction[ b ]]; ]; (* If there is no solution for the `simplest' equation, it terminates *) (* this recursive branch. *) If[Length[b] == 0, Return[{}]]; (* Transposes the solution: {{a,b},{c,d}} -> {{a,c}, {b,d}} . *) b = Transpose[b]; (* Adds constraints to the list, if there are any. .. && b && d *) c = (constraints && (And @@ #))& /@ b[[2]]; (* Drops constraints. *) b = b[[1]]; (* Applies the solution to the rest of the system and simplifies via *) (* a together and a numerator. *) a = Together[ Internal`DeactivateMessages[Rest[a] /. b, Power::infy, General::indet ] ]; a = Numerator[a]; (* Removes any cases which it eliminates one of the other equations. *) a = DeleteCases[a, {___, 0, ___}, {2}]; (* Takes the new solution, the one equation shorter system and the *) (* constraint and continues the recursive process. *) MapThread[ RecursivelyKeepTrack[#1, #2, primaryunknowns, secondaryunknowns, parameters, nonzeros, #3]&, {a, b, c}] ] (* ************************************************************************ *) (* Performs a factor list and then removes from it any nonzeros or numbers. *) (* Ex: -4 + 8a - 4b + 8ab - b^2 + 2ab^2 *) (* -> {{1, 1}, {-1 + 2*a, 1}, {2 + b, 2}} *) (* -> {{-1 + 2*a}, {2 + b}} *) FactorListAndCleanUp[ system_, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_] := Block[{a = FactorList /@ system}, (* DB:2/4/2003 Remove terms that are non-zero. *) a = a /. ({#^_:1, _} :> Sequence[] & /@ nonzeros); (* DB:2/4/2003 Remove terms that are numeric. *) a = a /. {_?NumericQ, _} :> Sequence[]; a = Flatten[Map[First[#]&, a, {2}]]; Union[a] ] (* ************************************************************************ *) (* New complexity analyzer by Douglas Baldwin on DB:2/3/2003 *) ComplexityLevel[expr_List, primaryunknowns_, secondaryunknowns_, parameters_]:= Max[ ComplexityLevel[#,primaryunknowns,secondaryunknowns,parameters] & /@ expr ] ComplexityLevel[expr_,primaryUnknowns_,secondaryUnknowns_,parameters_]:= Module[ { primaryComplexity = Select[primaryUnknowns, !FreeQ[expr,#]&], secondaryComplexity= Select[secondaryUnknowns, !FreeQ[expr,#]&], parameterComplexity= Select[parameters, !FreeQ[expr,#]&], exprLength = If[Head[Expand[expr]] === Plus, Length[Expand[expr]], 1], complexityDebug = False }, primaryComplexity = Exponent[expr, primaryComplexity] ; secondaryComplexity = Exponent[expr, secondaryComplexity]; parameterComplexity = Exponent[expr, parameterComplexity]; If[complexityDebug, Print["The precept sequence for the complexity analyzer is: "]; Print[ CleanUpFunction[{primaryComplexity, secondaryComplexity, parameterComplexity, exprLength} ] ]; ]; Return[ Min @@ Flatten[{1*primaryComplexity, 2*secondaryComplexity, 3*parameterComplexity}] + exprLength ]; ]; (* ************************************************************************ *) (* This computes the leafs of the complexity analysis tree. *) (* The complexity heuristic (smaller is better) is either *) (* a) the exponent of the expression if it is polynomial, or *) (* b) 100*b*level + size of expression, where b is the number *) (* of expressions which are of the level type (primary, secondary, *) (* parameter), and the level: primary = 1, secondary = 2, param = 3 *) SubComplexityLevel[expr_, a_, b_, level_] := Block[{c = Select[a, PolynomialQ[expr, #]&]}, Which[ Length[c] == 0, 100*b*level+LeafCount[expr], True, Min[Exponent[expr, c]] ] ] (* ************************************************************************ *) (* solves a single equation *) (* Takes the factor terms which have the primary unknowns and pass it *) (* on to the SubSolveStepByStep to be solved along with the secondary *) (* unknowns, parameters, constraints and solutions. *) SolveStepByStep[eqn_, primaryunknowns_, secondaryunknowns_, parameters_, constraints_, sol_] := Block[{a}, a = Union[Cases[{eqn}, q_ /; MemberQ[primaryunknowns, q], -1]]; ( SubSolveStepByStep[eqn, a, Flatten[{secondaryunknowns, parameters}], constraints, sol ] ) /; Length[a] != 0 ] (* ************************************************************************ *) (* Takes the factor terms which have the secondaryunknowns and pass it *) (* on to the SubSolveStepByStep to be solved along with parameters, *) (* constraints and solutions. *) SolveStepByStep[eqn_, primaryunknowns_, secondaryunknowns_, parameters_, constraints_, sol_] := Block[{a}, a = Union[Cases[{eqn}, q_ /; MemberQ[secondaryunknowns, q], -1]]; ( SubSolveStepByStep[eqn, a, parameters, constraints, sol] ) /; Length[a] != 0 ] (* ************************************************************************ *) (* Takes the factor terms which have the paramaters and pass it on to the *) (* SubSolveStepByStep to be solved along with cosntraints and solutions. *) SolveStepByStep[eqn_, primaryunknowns_, secondaryunknowns_, parameters_, constraints_, sol_] := Block[{a}, a = Union[Cases[{eqn}, q_ /; MemberQ[parameters, q], -1]]; ( SubSolveStepByStep[eqn, a, {}, constraints, sol] ) /; Length[a] != 0 ] (* ************************************************************************ *) (* Takes anything that isn't composed of primary unknowns, secondary *) (* unknowns, or parameters and returns an empty set. *) SolveStepByStep[___] := {} (* Takes the equation and primary unknowns (which are either THE primary *) (* unknowns, THE secondary unknowns, or THE parameters from above). *) (* 1. If the equation is polynomial in any of the unknowns, it selects *) (* them and sorts them by power/exponent and then includes the unknowns *) (* that the equation is not polynomial in. If the equation is not *) (* polynomial in any of the unknowns, it takes the order that they *) (* were given in. *) (* 2. Uses Reduce to solve the equation (along with find the constraints). *) (* 3. Then it converts the results to rules (same form as Solve) and *) (* returns an empty list if it is empty. *) (* 4. If there is a solution, it removes a != b type results from the *) (* form of the solution outputed by Reduce. *) (* 5. It then applies the results of the solution to the rest of the *) (* solution. *) (* 6. Formats the solution so that each sublist contains the latest solution *) (* branches along with the previous solutions in a flattened and unique *) (* list. *) SubSolveStepByStep[eqn_, primaryunknowns_, pars_, constraints_, sol_] := Block[{a, b, c}, a = Select[primaryunknowns, PolynomialQ[eqn, #]&]; If[Length[a] != 0, a = Sort[a, Exponent[eqn, #1] <= Exponent[eqn, #2]&]; a = Flatten[{a, Complement[primaryunknowns, a]}], a = primaryunknowns ]; c = myReduce[eqn == 0 && constraints, Flatten[{a, pars}]]; a = {ToRules[c]}; If[Length[a] == 0, Return[{}]]; c = Cases[{#}, p_ != q_, -1]& /@ If[Head[c] === Or, List @@ c, {c}]; b = Internal`DeactivateMessages[sol /. a, Power::infy, General::indet]; Union[Transpose[{Flatten /@ Thread[{b, a}, List, 2], c}]] ] (* ************************************************************************ *) End[] (* `Private` context *) SetAttributes[AnalyzeAndSolve, ReadProtected] EndPackage[] Print["Package PDESpecialSolutions.m was successfully loaded."]; Print["Version 2, last updated by Baldwin on February 28, 2005."]; (* ************************************************************************ *)