(* Last updated: March 21, 2010 at 15:40 by Hereman in Boulder *) (* Purpose: replace the old AnalyzeAndSolve -- dated January 8, 2010 *) (* with Doug Baldwin new AnalyzeAndSolve -- dated March 13, 2010 *) (* Fix various bugs. Corrections are marked with WH:03/17/2010 *) (* Previously updated: January 8, 2010 by Douglas Baldwin and Willy Hereman *) (* in Boulder at 16:50 *) (* See DB:01/08/2010 Added Version 7.0+ check. *) (* Previously updated: January 5, 2010 by Hereman in Boulder at 17:18 *) (* Avoid the duplication problem that comes from the symbolic test; *) (* changes marked with WH:01/05/2010 *) (* ************************************************************************ *) (* : Title : PDESpecialSolutionsV2.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 : Friday, January 8, 2009 at 16:50 by Hereman * : Previous update : Sunday, December 13, 2009 at 23:30 by Goktas and on January 23, 2009 by Hereman * : Based on code February 28, 2005 by Baldwin, Goktas, and Hereman * : 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\) up to 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 original 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 inconsistent 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. *) (* WH:01/23/2009 NumericTest is set to True by default, *) (* SymbolicTest is also set to True by default *) Options[PDESpecialSolutions] = { Verbose -> False, Form -> tanh, NumericTest -> True, SymbolicTest -> True, 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. *) (* Transforming the PDE into a nonlinear ODE in $T=Tanh$ or $S=Sech$ *) If[VerboseTest, Print["Transforming 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. *) (* Determining the maximal degree of the polynomial solutions. *) If[VerboseTest, Print["Determining 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. *) (* Deriving the algebraic system for the coefficients $a_{ij}, b_{ij}$. *) (* WH:01/05/2010 moved the print statement lower *) (* If[VerboseTest, Print["Deriving the nonlinear algebraic system for the coefficients."] ]; *) {mSoln, newSystem} = PDESpecialSolutionsBuildSystem[mSoln, system, vars, param, opts]; If[VerboseTest, Print["Deriving the nonlinear algebraic system for the coefficients."] ]; If[VerboseTest, Print["The nonlinear algebraic system is"]; Print[CleanUpFunction[newSystem /. myTrackingVariable[_]->1]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* DB:03/21/2003 Output algebraic system. *) If[!(SolveAlgebraicSystem /. {opts} /. Options[PDESpecialSolutions]), Return[CleanUpFunction[newSystem]]; ]; (* Step 4 in the paper. *) (* Solving the nonlinear parameterized algebraic system. *) If[VerboseTest, Print["Solving 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 *) (* DB: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. *) (* Building and testing the solitary wave solutions. *) If[VerboseTest, If[( !(NumericTest /. {opts} /. Options[PDESpecialSolutions]) && !(SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), Print["Building the travelling wave solutions (without testing)."], Print["Building and (numerically and/or symbolically) testing "<> "the solutions."] ] ]; solution = PDESpecialSolutionsBuildSolutions[solution, mSoln, eqns, funcs, vars, param, opts]; (* start if 10 *) If[VerboseTest, Print["Time Used:", TimeUsed[]-time]; Print[]; (* WH:01/05/2010 adjusted the print statements *) Which[ ( ( NumericTest /. {opts} /. Options[PDESpecialSolutions] ) && !( SymbolicTest /. {opts} /. Options[PDESpecialSolutions] ) ), Print["The final solutions were only tested numerically."], ( (SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) && !( NumericTest /. {opts} /. Options[PDESpecialSolutions] ) ), Print["The final solutions were only tested symbolically."], ( ( NumericTest /. {opts} /. Options[PDESpecialSolutions] ) && ( SymbolicTest /. {opts} /. Options[PDESpecialSolutions] ) ), Print["The final solutions were tested both numerically "<> "and symbolically."], ( !( NumericTest /. {opts} /. Options[PDESpecialSolutions] ) && !( SymbolicTest /. {opts} /. Options[PDESpecialSolutions] ) ), Print["The list of possible final solutions (without "<> "testing them numerically or symbolically):"] ] (* end of which *) ]; (* end if 10 *) (* Returns solutions to user. *) Return[ If[InputForm /. {opts} /. Options[PDESpecialSolutions], InputForm[CleanUpFunction[solution]], CleanUpFunction[solution] ] ] ) /; testarg ]; (* end of block PDESpecialSolutions *) (* ************************************************************************ *) (* : 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 = (system /. (p_ -> q_) :> (p -> Expand[q])); (* 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[(system /. (p_ -> q_) :> (p -> Factor[q]))]; ]; (* end of block PDESpecialSolutionsVarChange *) (* ************************************************************************ *) (* : 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]) =!= {}; (* end of module PDESpecialSolutionsmSolver *) PDESpecialSolutionsmSolver[system_List, theForm_, opts___?OptionQ] := Module[{mSolverDebug = (* Debugging Boolean 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 function. *) 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:05/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:03/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 original 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 original 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]) === {}; (* end of module PDESpecialSolutionsmSolver *) (* Breaks up the system into the correct form. *) mSolve`SystemGeneration[system_, theForm_, myTrackingVariableMax_, opts___?OptionQ] := Block[{SystemGenerationDebug = (* Debugging Boolean 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:02/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]; ]; (* end of block mSolve`SystemGeneration *) (* 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]; ]; (* end of module mSolve`SimplifySystem *) mSolve`ListFormation[eqnList0_, Form_, funcRules_, myTrackingVariableMax_, opts___?OptionQ]:= Block[{ mPowerListFormationDebug = (* Debugging Boolean 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 @@ #, {#}]& /@ #& /@ (eqnList /. (p_ -> q_) :> (p -> Expand[q])); 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:05/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]; ]; (* end of block mSolve`ListFormation *) mSolve`Simplification[mList0_, opts___?OptionQ]:= Module[{mPowerListSimplificationDebug = (* Debugging Boolean 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]; ]; (* end of module mSolve`Simplification *) mSolve`RulesSolver[mList0_, myMList_, opts___?OptionQ]:= Module[{rulesSolverDebug = (* Debugging Boolean 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]; ]; (* end of module mSolve`RulesSolver *) mSolve`PowerSolver[mRules0_, myMList_, opts___?OptionQ]:= Module[{powerSolverDebug = (* Debugging Boolean 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:05/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]; ]; (* end of module mSolve`PowerSolver *) mSolve`SystemCleanUp[eqnList0_, mList_, mSoln0_, myMList_, F_, opts___?OptionQ]:= Module[{systemCleanUpDebug = (* Debugging Boolean 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:05/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:05/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 = ( Plus @@ Coefficient[(eqnList /. #), F^Max[mList /. mSoln] ]& /@ mSoln ) /. (p_ -> q_) :> (p -> Factor[q]); If[systemCleanUpDebug, Print["mSoln, Checking in eqnList:"]; Print[CleanUpFunction[eqnList]]; ]; (* Returns the good solutions. *) Return[mSoln]; ]; (* end of module mSolve`SystemCleanUp *) (* Taken from Painleve Test code DB:05/15/2004. *) mSolve`FixFreeM[ mSoln_List, myMList_List, opts___?OptionQ ] := Module[{ fixFreeMDebug = (* Debugging Boolean flag *) (mSolverVerbose /. {opts} /. Options[PDESpecialSolutions]), mSoln0, mFree, mFreeValues, mFixedValues, mValues }, (* Protected Local Variables *) (* Remove symbolic solutions, based on the code in SystemCleanUp. *) (* DB:07/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:02/08/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 up to max:"]; Print[CleanUpFunction[ mValues ]]; ]; mValues = ( Rule @@ #& /@ Transpose[ { myMList, # } ] )& /@ mValues; If[fixFreeMDebug, Print["The fixed values reformatted:"]; Print[CleanUpFunction[ mValues ]]; ]; Return[ mValues ]; ]; Return[mSoln0]; ]; (* end of module mSolve`FixFreeM *) (* Finds free solutions that result from inequalities. *) (* Added direct from old code DB:11/20/2002. *) 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 ] ] ]; (* end module mSolve`GenerateAlternativeSolutions *) (* ************************************************************************ *) (* : Title : PDESpecialSolutionsBuildSystem. * : Author : Douglas Baldwin * : Date : 05/24/2001 *) (* : Update : 01/12/2003 by Douglas Baldwin * Added debug statements and highest order equation solving. * Changes marked with DB:Date *) (* : Update : 01/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, pureuRules, 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 ] ] (* closes evaluate *) ) )& /@ mSoln; If[buildSystemDebug, Print["Before conversion in pure functions, these are the uRules : "]; Print[CleanUpFunction[ uRules ]]; ]; (* WH: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]))&; (* WH:01/23/2009 added the following line to see the format of pureuRules *) pureuRules = toPure[#]& /@ uRules; If[buildSystemDebug, Print["After conversion into pure function rules, the pureuRules are: "]; Print[CleanUpFunction[ pureuRules ]]; ]; (* Applies the expressions in T to the system and removes the *) (* tracking variables by setting these to 1 *) newSystem = Expand[ system //. Append[toPure[#], myTrackingVariable[_]->1] ]& /@ uRules; If[buildSystemDebug, Print["Apply the expressions in T to the system 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:02/04/2003 Put mod with param instead *) (* of c_is.*) waveparameters = Table[c[i], {i, Length[vars]}]; (* DB:02/04/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:01/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 = newSystem /. (p_ -> q_) :> (p -> Expand[q]); (* Determines the maximum value of F in each of the cases. *) maxF = Map[ Exponent[#, Form /. funcRules ]&, newSystem ] /. -Infinity -> 0 (* DB:03/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 DB:01/11/2003. *) (* Removed on Tuesday, January 14, 2003 DB:01/23/2003 *) (* 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:01/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:01/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:01/12/2003 highestOrderSytem is then formatted to *) (* input into AnalyzeAndSolve *) highestOrderSystem = Transpose[ {highestOrderSystem, unknowns, Table[waveparameters, {Length[newSystem]}], Table[param, {Length[newSystem]}], nonzeros } ] /. (p_ -> q_) :> (p -> Factor[q]); If[buildSystemDebug, Print["The highest order system formatted for AnalyzeAndSolve: "]; Print[CleanUpFunction[ highestOrderSystem ]]; ]; (* DB:01/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:02/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:01/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 all True or False items:"]; Print[CleanUpFunction[ newSystem ]]; ]; (* DB:01/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 = ( (* 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 DB:02/02/2003. *) 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 } ] ) /. (p_ -> q_) :> (p -> Factor[q]), (* Combines into lists containing {newSystem, unknowns, *) (* waveparameters, param, nonzeros} for the solver. *) newSystem = ( Transpose[ {newSystem, unknowns, Table[waveparameters, {Length[newSystem]}], Table[param, {Length[newSystem]}], nonzeros } ] ) /. (p_ -> q_) :> (p -> Factor[q]); ]; If[buildSystemDebug, Print["The system formatted for AnalyzeAndSolve: "]; Print[CleanUpFunction[ newSystem ]]; ]; (* DB:03/21/2003 Allows mSoln to be updated in main. *) newSystem = Thread[{mSoln, newSystem}]; If[buildSystemDebug, Print["Threading mSoln"]; Print[CleanUpFunction[ newSystem]]; ]; (* DB:03/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:03/21/2003 *) If[Length[newSystem] === 0, Message[ PDESpecialSolutionsBuildSystem::fail ]; Abort[ ] ]; (* DB:03/21/2003 Transpose added after threading mSoln above. *) Return[ Transpose[newSystem] ]; ]; (* end of block PDESpecialSolutionsBuildSystem *) (* ************************************************************************ *) (* : 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 in *) (* numerical test. *) symbolicTestSolutions, (* tests the validity of solutions. *) (* in symbolic test. *) newSymbolicTestSolutions, (* used in evaluation of expresssions *) (* with roots of cubics, quartics, etc. *) evalNewSymbolicTestSolutions, (* ibid. *) mySymbol, (* ibid. *) variousRoots, (* ibid. *) ruleForMySymbol, (* ibid. *) rulesForMySymbols, (* ibid. *) tempSymbolicTestSolutions, (* used to control simplifications with *) (* simplify function *) 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 = solution /. (p_ -> q_) :> (p -> Factor[q]); (* 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]), (* WH:01/05/2010 the patterns on left may not be recognized *) (1-Sech[xx__]^2) -> Tanh[xx]^2, (1-Tanh[xx__]^2) -> Sech[xx]^2, (* begin replacement, both below created an infinite loop *) (* Tanh[xx__]^2 -> 1-Sech[xx]^2, *) (* Sech[xx__]^2 -> 1-Tanh[xx]^2, *) (* end replacement *) 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[ (( (( # //. mySimplificationRules ) /. (p_ -> q_) :> (p -> Factor[q])) //. mySimplificationRules ) /. (p_ -> q_) :> (p -> Expand[q]))&, solutionRules, 3 ]; (* DB:03/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[ ( solutionRules /. {CN -> JacobiCN, SN -> JacobiSN} ) /. (p_ -> q_) :> (p -> Factor[q]) ] ]; ]; (* Prints out a warning, if it's taking a long time. *) If[TimeUsed[]-time>3, Print["Still building the exact solutions, please wait."] ]; (* WH:03/17/2010 formatted the cellprint lines were too long, *) (* or print too small in Mathematica v. 7.0 *) (* Depending on the NumericTest option, it either tests the *) (* solutions numerically, or it doesn't. *) (* start if 0 *) If[ ( !(NumericTest /. {opts} /. Options[PDESpecialSolutions]) && !(SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), CellPrint[ Cell[ "These solutions are not being tested numerically or symbolically", "Message" ] ]; CellPrint[ Cell[ "To test the solutions set either the NumericTest option to True ", "Message" ] ]; CellPrint[ Cell[ "or set the SymbolicTest option to True, or both. ", "Message" ] ]; (* WH:03/17/2010 formatted the cellprint lines were too long, *) (* or print too small in Mathematica v. 7.0 *) If[Length[warn]>0, CellPrint[ Cell[ "The following simplification rules are being used:", "Message" ] ]; CellPrint[ Cell[ ToString[Union[warn]], "Message" ] ]; ]; (* end if *) Return[ Map[Union, ( ( solutionRules /. {CN -> JacobiCN, SN -> JacobiSN} ) /. (p_ -> q_) :> (p -> Expand[q]) ) /. (p_ -> q_) :> (p -> Factor[q]), 2 ] ] ]; (* end if 0 *) (* Converts form of solutions given by solver to a pure function. *) toPure = (# /. (a__[var__]->temp__):> (a:>Function[{var}, temp]))&; (* WH:01/23/2009 these tests can be removed later *) (* Print["At xxxxxxx0, this is NumericTest :"]; Print[NumericTest]; Print["At xxxxxxx1, this is opts :"]; Print[opts]; Print["At xxxxxxx2, this is Options[PDESpecialSolutions] :"]; Print[Options[PDESpecialSolutions]]; Print["At xxxxxxx3, this is Verbose :"]; Print[Verbose]; Print["At xxxxxxx4, this is "<> "NumericTest /. {opts} /. Options[PDESpecialSolutions] :"]; Print[NumericTest /. {opts} /. Options[PDESpecialSolutions]]; Print["At xxxxxxx5, this is "<> "!(Verbose /. {opts} /. Options[PDESpecialSolutions]): "]; Print[!(Verbose /. {opts} /. Options[PDESpecialSolutions])]; *) (* WH:01/23/2009 always print the message about numerically testing *) (* If statement for the numerically test option. *) If[NumericTest /. {opts} /. Options[PDESpecialSolutions], (* Sends message to user. *) If[(Verbose /. {opts} /. Options[PDESpecialSolutions]), Print["Numerically testing the solutions."] ]; (* was: If[!(VerboseTest /. {opts} /. Options[PDESolve]), Print["Numerically test the solutions."] ]; *) (* WH:01/04/2010 Added a print statement *) (* WH:03/17/2010 solutionRules enters numerical test *) If[VerboseTest, Print["The following solutions will be numerically tested:"]; Print[CleanUpFunction[solutionRules]] ]; (* WH:03/17/2010 does this change format to a list of double lists? *) (* code Baldwin: *) (* Tests numerically to make sure the above solutions are valid. *) numericTestSolutions = {(eqns /. a__ == b__:>a-b) /. toPure[(# /. (p_ -> q_) :> (p -> TrigToExp[q]))], #}& /@ #& /@ solutionRules; (* WH:03/17/2010 *) (* Print["At Pt. R1, numericTestSolutions: "]; Print[CleanUpFunction[numericTestSolutions]]; *) 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)& /@ ( (( #[[1]] /. randomVarRules ) /. (p_ -> q_) :> (p -> Expand[q])) /. randomAijRules /. {CN -> JacobiCN, SN-> JacobiSN} ) ), And @@ ( Abs[ # ] < 10^(-$MachinePrecision/2)& /@ ( (( #[[1]] /. ( randomVarRules /. (a_ -> b_) :> (a -> -b) ) ) /. (p_ -> q_) :> (p -> Expand[q])) /. ( randomAijRules /. (a_ -> b_) :> (a -> -b) ) /. {CN -> JacobiCN, SN-> JacobiSN} ) ) ], # }& /@ #& /@ numericTestSolutions; (* WH:03/17/2010 *) (* Print["At Pt. R2, numericTestSolutions: "]; Print[CleanUpFunction[numericTestSolutions]]; *) (* If it works for any of the trials, it keeps the solution. *) solns = numericTestSolutions /. { {True, a_} :> a[[2]], {False, _} :> Sequence[] }; (* WH:03/17/2010 *) (* Print["At Pt. R3, solns: "]; Print[CleanUpFunction[solns]]; Print["The above soln are the solutions that PASSED the NUMERICAL test"]; *) (* WH:01/04/2010 Added a print statement *) If[VerboseTest, Print["The following solutions passed the numerical test:"]; Print[CleanUpFunction[solns]] ]; (* WH:03/17/2010 *) (* Print["At Pt. R3-bis, solutionRules: "]; Print[CleanUpFunction[solutionRules]]; *) (* WH:03/17/2010 change the lists solns and solutionRules artifically *) (* to test if complement works *) (* solns = solns /. { 15 -> N[Pi,4] }; *) (* Print["At Pt. R3-tris, solns in raw form: "]; Print[solns]; *) (* Print["At Pt. R3-quatro (after modification), solns: "]; Print[CleanUpFunction[solns]]; Print["length of solns (after modification):", Length[solns[[1]]]]; *) (* solutionRules = solutionRules /. { 15 -> N[Pi,6] }; *) (* Print["At Pt. R3-5, solutionRules in raw form: "]; Print[solutionRules]; *) (* Print["At Pt. R3-6 (after modification), solutionRules: "]; Print[CleanUpFunction[solutionRules]]; Print["length of solutionRules (after modification):", Length[solutionRules[[1]]]]; *) (* Print["Now take complement of solutionRules and solns"]; *) (* WH:03/17/2010 the complement below malfunctioned *) (* Baldwin's code was: *) (* numericTestSolutions = Complement[solutionRules, solns]; *) (* WH:03/17/2010 replaced by *) numericTestSolutions = { Complement[solutionRules[[1]], solns[[1]]] }; (* WH:03/17/2010 *) (* Print["At Pt. R4, coming from taking the complement, numericTestSolutions: "]; Print[CleanUpFunction[numericTestSolutions]]; Print["length of numericTestSolutions:", Length[numericTestSolutions[[1]]]]; *) (* Print["The above numericTestSolution are the solutions that DID NOT PASS "<> " the NUMERICAL test"]; *) (* WH:03/17/2010 formatted the cellprint lines were too long, *) (* or print too small in Mathematica v. 7.0 *) If[(Map[Union, numericTestSolutions, 3] //. {{}}->{}) != {}, CellPrint[ Cell[ "After subsitution of the solutions into the PDE, the resulting "<> "expressions did not simplify to be less than 10^(" <> ToString[-$MachinePrecision/2] <> ") ", "Message" ] ]; CellPrint[ Cell[ "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. ", "Message" ] ]; CellPrint[ Cell[ "Please test the corresponding solutions by hand or interactively "<> "with Mathematica.", "Message" ] ]; Print[CleanUpFunction[#]]& /@ numericTestSolutions; ]; ]; (* End of numeric test if statement. *) (* WH:01/23/2009 these tests can be removed later *) (* Print["At yyyyyyy0, this is SymbolicTest :"]; Print[NumericTest]; Print["At yyyyyyy1, this is opts :"]; Print[opts]; Print["At yyyyyyy2, this is Options[PDESpecialSolutions] :"]; Print[Options[PDESpecialSolutions]]; Print["At yyyyyyy3, this is Verbose :"]; Print[Verbose]; Print["At yyyyyyy4, this is "<> "NumericTest /. {opts} /. Options[PDESpecialSolutions] :"]; Print[NumericTest /. {opts} /. Options[PDESpecialSolutions]]; Print["At yyyyyyy5, this is "<> "!(Verbose /. {opts} /. Options[PDESpecialSolutions]): "]; Print[!(Verbose /. {opts} /. Options[PDESpecialSolutions])]; *) (* WH:01/23/2009 always print the message about symbolic testing *) (* The second testing If statement. *) If[SymbolicTest /. {opts} /. Options[PDESpecialSolutions], (* Sends message to user. *) If[(Verbose /. {opts} /. Options[PDESpecialSolutions]), Print["Symbolically testing the solutions."] ]; (* was: If[!(VerboseTest /. {opts} /. Options[PDESolve]), Print["Symbolically test the solutions."] ]; *) (* WH:01/04/2010 Added a line of code and a print statement. *) (* Under the default setting, solns would come out of the numeric test. *) (* therefore, the solns would be {} if numeric test is set to False. *) (* So, if numeric test is false, set solns = solutionRules, but simplify *) (* solutionRules a bit *) If[!(NumericTest /. {opts} /. Options[PDESpecialSolutions]), solns = ( ( solutionRules /. {CN -> JacobiCN, SN -> JacobiSN} ) /. (p_ -> q_) :> (p -> Expand[q]) ) /. (p_ -> q_) :> (p -> Factor[q]) ]; (* WH:03/17/2010 only test solutions (soln) that passed the numeric test. *) (* Added a print statement to specify which solutions are being tested. *) If[VerboseTest, If[(NumericTest /. {opts} /. Options[DDESpecialSolutions]), Print["Only the solutions that passed the numeric test will "<> "be tested symbolically."]; Print["To test all solutions, set the option NumericTest to False."] ] ]; If[VerboseTest, Print["The following solutions will be tested symbolically:"]; Print[CleanUpFunction[solns]] ]; (* WH:03/17/2007 artificially mess up one of the solutions to test if the *) (* code handles that correctly *) (* solutionRules = solutionRules /. { 15 -> N[Pi,6] }; *) (* Print["At Pt. U1 (after modification), solutionRules: "]; Print[CleanUpFunction[solutionRules]]; Print["length of solutionRules (after modification):", Length[solutionRules[[1]]]]; *) (* WH:03/17/2010 only test the solutions solns that passed the numeric test *) (* if the numeric test was skipped then solns = solutionRules *) (* This sets up the input for the next block of code. *) (* In that, it takes the equation 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. *) (* WH:03/17/2010 replaced solutionRules by solns *) symbolicTestSolutions = {(eqns /. a__ == b__:>a-b) /. toPure[(# /. (p_ -> q_) :> (p -> TrigToExp[q]))], #}& /@ #& /@ ( (* solutionRules *) solns /. {CN -> JacobiCN, SN-> JacobiSN}); (* WH:03/17/2010, check for format changes! *) (* Print["At Pt. U2, solutionRules: "]; Print[CleanUpFunction[solutionRules]]; *) (* Attempts to simplify the expressions as much as possible. *) (* WH:01/05/2010 works fine on cases without roots of cubics, etc. *) symbolicTestSolutions = ExpToTrig[ FixedPoint[ (( (# /. (p_ -> q_) :> (p -> Expand[q])) /. mySimplificationRules ) /. (p_ -> q_) :> (p -> Factor[q])) /. mySimplificationRules &, #, 10]& /@ symbolicTestSolutions ]; (* WH:03/17/2010 *) (* Print["At Pt. U3, after applying FixedPoint, symbolicTestSolutions: "]; Print[CleanUpFunction[symbolicTestSolutions]]; Print["Check the format in comparison with solutionRules"]; *) (* Print["At Pt. B1, after fixedpoint algorithm, symbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; *) (* WH:01/05/2010 added new piece of code to deal with expressions *) (* involving unevaluated roots of cubics, quartics, etc. *) (* start if 2 *) If[( Map[Union, symbolicTestSolutions,3] //. {{}}->{} ) != {} && ( FreeQ[symbolicTestSolutions, Root] === False ), Print["Simplifying the remaining expressions (involving roots "<> "of cubics, quartics, etc.) "]; (* Print["At Pt. B2, symbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; *) newSymbolicTestSolutions = symbolicTestSolutions /. { Root[a__,i_] -> mySymbol[i] }; (* Print["At Pt. B3, newSymbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ newSymbolicTestSolutions; *) variousRoots = Union[Cases[symbolicTestSolutions, Root[a__,i_],Infinity]]; (* Print["At Pt. B4, variousRoots = "]; Print[CleanUpFunction[#]]& /@ variousRoots; *) ruleForMySymbol[i_] := { mySymbol[i] -> Part[variousRoots,i] }; (* Print["At Pt. B5, ruleForMySymbol[i] = "]; Print[CleanUpFunction[ruleForMySymbol[i]]]; *) rulesForMySymbols = Flatten[Table[ ruleForMySymbol[i],{i,1,Length[variousRoots]}]]; (* Print["At Pt. B6, rulesForMySymbols = "]; Print[CleanUpFunction[#]]& /@ rulesForMySymbols; *) evalNewSymbolicTestSolutions = Map[Simplify,newSymbolicTestSolutions /. rulesForMySymbols]; (* Print["At Pt. B7, evalNewSymbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ evalNewSymbolicTestSolutions; *) symbolicTestSolutions= evalNewSymbolicTestSolutions; (* Print["At Pt. B9, symbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; *) ]; (* end if 2 *) (* WH:01/05/2010 reluctantly added a simplify because expressions *) (* involving sech did not fully simplify *) tempSymbolicTestSolutions = DeleteCases[ DeleteCases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ], {} ]; (* Print["At Pt. C1, tempSymbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ tempSymbolicTestSolutions; *) If[tempSymbolicTestSolutions != {}, symbolicTestSolutions = Map[Simplify, symbolicTestSolutions]; (* Print["At Pt. C4, after mapping simplify, symbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions *) ]; (* for debugging only *) (* If[( Map[Union, symbolicTestSolutions,3] //. {{}}->{} ) != {}, result = Map[Union, symbolicTestSolutions,3] //. {{}}->{}; Print["At Pt. C2, test for not equal to {}, result ="]; Print[CleanUpFunction[#]]& /@ result; symbolicTestSolutions = Map[Simplify, symbolicTestSolutions]; Print["At Pt. C3, after mapping simplify, symbolicTestSolutions ="]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; ]; *) (* for old code January 23, 2009 WH:01/05/2010 *) (* Print["At Pt. 1, before joining the symbolically tested solutions, solns = "]; Print[CleanUpFunction[solns]]; *) (* WH:01/05/2010 do no longer join the original (not-tested) and the tested *) (* solutions. *) (* (* Pulls off 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} ]; *) (* by the following *) (* solns = Cases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ] //. {{(0)..}, a_List} :> {a} ; *) (* WH:03/17/2010 if a_List -> {a} *) (* Print["At Pt. U4, if a_List -> {a}, then solns: "]; Print[CleanUpFunction[solns]]; *) (* WH:03/17/2010 change a_List -> {a}, into a_List -> a *) solns = Cases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ] //. {{(0)..}, a_List} :> a ; (* WH:03/17/2010 *) (* Print["At Pt. U5, if a_List -> a, then solns: "]; Print[CleanUpFunction[solns]]; Print["At Pt. U6, continue with this format for the rest of the compuations"]; *) (* in former code of January 23, 2009 WH:01/05/2010 *) (* Print["At Pt. 2, after joining the symbolically tested solutions, solns = "]; Print[CleanUpFunction[solns]]; *) If[VerboseTest, Print["The following solutions (in factored form) passed "<> "the symbolic test:"]; Print[CleanUpFunction[solns]] ]; (* WH:01/05/2010 *) (* in former code of January 23, 2009, applying a Union here did not work *) (* Print["At Pt. 3, after applying a Union, solns = "]; Print[CleanUpFunction[solns]]; *) (* Removes all the good cases, so we may output the bad *) (* ones to the user. *) (* Print["At Pt. A1, before deleting good cases, symbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; *) symbolicTestSolutions = DeleteCases[ DeleteCases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ], {} ]; (* Print["At Pt. A2, after deleting good cases, symbolicTestSolutions = "]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; *) (* WH:03/17/2010 formatted the cellprint lines were too long, *) (* or print too small in Mathematica v. 7.0 *) (* start if 1 *) If[ ( Map[Union, symbolicTestSolutions,3] //. {{}}->{} ) != {}, (* Sends out the left over (questionable) solutions to *) (* the user for hand inspection. *) CellPrint[ Cell[ "After substitution of the solution into the PDE,"<> " the resulting expressions did not simplify to zero. ", "Message" ] ]; CellPrint[ Cell[ "Please test the corresponding solutions by hand or interactively "<> "with Mathematica.", "Message" ] ]; Print[CleanUpFunction[#]]& /@ symbolicTestSolutions; ]; (* end if 1 *) ]; (* end if symbolic test *) (* WH:03/17/2010 formatted the cellprint lines were too long, *) (* or print too small in Mathematica v. 7.0 *) If[Length[warn]>0, CellPrint[ Cell[ "The following simplification rules are being used ", "Message" ] ]; CellPrint[ Cell[ ToString[Union[warn]], "Message" ] ] ]; (* end if *) (* Returns the solution rules *) (* Print["At Pt. 5, before applying the simplificatons, solns = "]; Print[CleanUpFunction[solns]]; *) (* for debugging *) (* seeSolns = Union[ ( ( solns /. {CN -> JacobiCN, SN -> JacobiSN} ) /. (p_ -> q_) :> (p -> Expand[q]) ) /. (p_ -> q_) :> (p -> Factor[q]) ]; *) (* Print["At Pt. 6, after applying the simplificatons, returned are seeSolns = "]; Print[CleanUpFunction[seeSolns]]; *) Return[ Union[ ( ( solns /. {CN -> JacobiCN, SN-> JacobiSN} ) /. (p_ -> q_) :> (p -> Expand[q]) ) /. (p_ -> q_) :> (p -> Factor[q]) ] ]; ]; (* end of block PDESpecialSolutionsBuildSolutions *) (* ************************************************************************ *) (* : 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 original expression is: "]; Print[CleanUpFunction[ theList /. (p_ -> q_) :> (p -> Factor[q]) ]]; ]; (* 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 of module StripNonZero *) End[] (* `Private` context *) Protect[PDESpecialSolutions] SetAttributes[PDESpecialSolutions, ReadProtected] EndPackage[] (* START OLD ANALYZEANDSOLVE (* ************************************************************************ *) (* 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 *) (* DB:01/13/2003 *) BeginPackage["Algebra`AnalyzeAndSolve`"] Unprotect[AnalyzeAndSolve] AnalyzeAndSolve::soln = "The solution `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[] ] ] (* end of module Internal`DeactivateMessages *) ]; (* end of if versionnumber < 4 *) If[$VersionNumber >= 5, (* DB:01/08/2010 Added Version 7.0+ check. *) If[$VersionNumer >= 7, SetSystemOptions[ "ReduceOptions" -> {"ReorderVariables" -> False} ], Developer`SetSystemOptions[ "ReduceOptions" -> {"ReorderVariables" -> False} ] ]; (* end if versionnumber >= 7 *) myReduce[eqn_, vars_] := LogicalExpand[ Reduce[eqn, Reverse[vars]] /. Equal -> myEqual ] /. myEqual -> Equal, myReduce[eqn_, vars_] := Reduce[eqn, vars, Sort -> False] ]; (* end if versionnumber >= 5 *) (* : Title : CleanUpFunction *) (* : Author : Douglas Baldwin *) (* : Summary : Remove "Algebra`AnalyzeAndSolve`Private`" from output. *) (* : Added for debugging for AnalyzeAndSolve DB:01/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] ]; (* end of block AnalyzeAndSolve *) (* ************************************************************************ *) (* Collects all the solutions at the end into globalsol. *) (* Original code . . . replaced by the two following functions. *) RecursivelyKeepTrack[{}, sol_, __] := (AppendTo[globalsol, sol]; {}) (* Terminates that recursive branch if the solution explodes 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:05/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:01/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:02/04/2003 Added expand. *) b = b /. (p_ -> q_) :> (p -> Factor[q]); 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}] ]; (* end of block RecursivelyKeepTrack *) (* ************************************************************************ *) (* 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:02/04/2003 Remove terms that are non-zero. *) a = a /. ({#^_:1, _} :> Sequence[] & /@ nonzeros); (* DB:02/04/2003 Remove terms that are numeric. *) a = a /. {_?NumericQ, _} :> Sequence[]; a = Flatten[Map[First[#]&, a, {2}]]; Union[a] ]; (* end of block FactorListAndCleanUp *) (* ************************************************************************ *) (* New complexity analyzer by Douglas Baldwin DB:02/03/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 ]; ]; (* end of module ComplexityLevel *) (* WH: March 13, 2010 put sets of parenthesis around (a) and (b) *) (* ************************************************************************ *) (* This computes the leaves 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]] ] ]; (* end of block SubComplexityLevel *) (* ************************************************************************ *) (* 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 ]; (* end of block SolveStepByStep 1 *) (* ************************************************************************ *) (* 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 ]; (* end of block SolveStepByStep 2 *) (* ************************************************************************ *) (* Takes the factor terms which have the parameters and pass it on to the *) (* SubSolveStepByStep to be solved along with constraints 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 ]; (* end of block SolveStepByStep 3 *) (* ************************************************************************ *) (* 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 outputted 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 of block SubSolveStepByStep *) (* ************************************************************************ *) End[] (* `Private` context *) SetAttributes[AnalyzeAndSolve, ReadProtected] EndPackage[] END OLD ANALYZEANDSOLVE *) (* START NEW ANALYZEANDSOLVE *) (* ************************************************************************ *) (* : Title : AnalyzeAndSolveV2.m * : Author : Douglas Baldwin, Unal Goktas, and Willy Hereman * : Email : solve@douglasbaldwin.com * : Date Modified : Thursday, March 4, 2010 by Douglas Baldwin * : History : First developed by Unal Goktas (WRI) and Willy Hereman and * written by Unal Goktas in October/November 2000; * comments added and complexity analyzer redesigned by Douglas Baldwin * in February 2003; * redesigned by DB in February 2007 to use depth-first search rather than * a recursive breadth-first search; * DB finished redesign in March 2010. *) (* ************************************************************************ *) 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[] ] ] (* end of module Internal`DeactivateMessages *) ]; (* end of if versionnumber < 4 *) If[$VersionNumber >= 5, (* DB:1/8/2010 Added Version 7.0+ check. *) If[$VersionNumer >= 7, SetSystemOptions[ "ReduceOptions" -> {"ReorderVariables" -> False} ], Developer`SetSystemOptions[ "ReduceOptions" -> {"ReorderVariables" -> False} ]]; (* end if versionnumber >= 7 *) myReduce[eqn_, vars_] := LogicalExpand[ Reduce[eqn, Reverse[vars]] /. Equal -> myEqual ] /. myEqual -> Equal, myReduce[eqn_, vars_] := Reduce[eqn, vars, Sort -> False] ]; (* end if versionnumber >= 5 *) (* ************************************************************************ *) (* : 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`"->"", "Calculus`DDESpecialSolutions`Private`"->"", "Algebra`AnalyzeAndSolve`Private`"->""} ] ]& (* ************************************************************************ *) (* : Title : AnalyzeAndSolve *) (* : Author : Douglas Baldwin *) (* : Summary : Takes a list of equations which depend on primary unknowns, * secondary unknowns, and parameters and solves for the unknowns and * parameters in an depth-first search which takes advantage of the given * nonzero unknowns to simplify the expressions. *) AnalyzeAndSolve[system_, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_] := Block[{ sys = system /. Equal[a_,b_]:>a-b, (* A local copy of system. *) constraints, (* Constraints on the possible solutions. *) theList, (* The search tree. *) soln = {}, (* The solution. *) currentSystem, (* The system of equations to be solved. *) currentSolutions, (* Expressions for known variables. *) currentSoln, (* The new solution from the simplest equation. *) currentConstraints, (* The constraints on future solutions. *) theNextNode, (* The next node to be expanded. *) depthFirstSolveDebug = False, recursiveDebug = False }, (* Protected local variables. *) sys = Numerator[Together[({#}& /@ sys)]]; (* Takes the nonzeros given by the user and puts them in the form *) (* needed by Reduce. *) constraints = (And @@ ((# != 0 &) /@ nonzeros)); (* The first node of the search tree. *) theList = {{{sys,{},constraints}}}; While[True, (* Begin While A1 *) If[( !FreeQ[theList[[-1,1,1]], DirectedInfinity] || !FreeQ[theList[[-1,1,1]], Indeterminate] || !FreeQ[theList[[-1,1,1]], False] ), (* Goes back to the next branch. *) While[ Length[theList[[-1]]] === 1, theList = Drop[theList,-1]; (* If all branches have been searched, return soln. *) If[ Length[theList] === 0, Return[soln] ]; ]; (* Removes the branch just solved for. *) theList = Append[Drop[theList,-1], Rest[theList[[-1]]]]; If[depthFirstSolveDebug, Print["The pruned search tree:"]; Print[CleanUpFunction[ theList ]]; ]; ]; If[theList[[-1,1,1]] === {}, (* Begin If A2 *) (* * * * * * * * * * * * * * If TRUE of If A2 * * * * * * * * * * * * * * *) If[depthFirstSolveDebug, Print["The current search tree is:"]; Print[CleanUpFunction[ theList ]]; ]; (* Adds new solution to soln. *) AppendTo[soln, theList[[-1,1,2]] ]; If[depthFirstSolveDebug, Print["The complete set of solutions is:"]; Print[CleanUpFunction[ soln ]]; ]; (* Goes back to the next branch. *) While[ Length[Last[theList]] === 1, theList = Drop[theList,-1]; (* If all branches have been searched, return soln. *) If[ Length[theList] === 0, Return[soln] ]; ]; (* Removes the branch just solved for. *) theList = Append[Drop[theList,-1], Rest[theList[[-1]]]]; If[depthFirstSolveDebug, Print["The pruned search tree:"]; Print[CleanUpFunction[ theList ]]; ], (* * * * * * * * * * * * * * If FALSE of If A2 * * * * * * * * * * * * * * *) (* Applies the current solutions to the system. *) currentSystem = FactorListAndCleanUp[#, primaryunknowns, secondaryunknowns, parameters, nonzeros]& /@ theList[[-1,1,1]]; (* Pulls off the first solution of each list. *) currentSolutions = theList[[-1,1,2]]; (* Pulls off the corresponding contraints on *) (* the above solutions. *) currentConstraints = theList[[-1,1,3]]; If[depthFirstSolveDebug, Print["The last entry in theList is:"]; Print[CleanUpFunction[ Last[theList] ]]; Print["The current system is:"]; Print[CleanUpFunction[ currentSystem ]]; Print["The current list of solutions is:"]; Print[CleanUpFunction[ currentSolutions ]]; Print["The current list of contraints is:"]; Print[CleanUpFunction[ currentConstraints ]]; ]; (* Removes duplicates. Message added DB:5/23/2003 *) currentSystem = Union[currentSystem] /. {} :> (Message[AnalyzeAndSolve::soln, CleanUpFunction[currentSolutions]];{}); If[depthFirstSolveDebug, Print["The system with duplicates removed is:"]; Print[CleanUpFunction[ currentSystem ]]; ]; (* Sorts the system by the "ComplexityLevel" heuristic. *) (* DB:3/2/2010 Changed back to Sort[List, Test]. *) currentSystem = Sort[currentSystem, (ComplexityLevel[#1, primaryunknowns, secondaryunknowns, parameters] <= ComplexityLevel[#2, primaryunknowns, secondaryunknowns, parameters]&)]; (* DB:1/14/2003 Print statement added for debugging. *) If[depthFirstSolveDebug || recursiveDebug, Print["The order of the equations in AnalyzeAndSolve is: "]; Print[ CleanUpFunction[ {ComplexityLevel[#1, primaryunknowns, secondaryunknowns, parameters],#1}& /@ currentSystem ] ]; ]; (* Takes the first equation of the system of equations (which *) (* is the simplest according to the complexity heuristic) *) (* and then solves it for only one variable. *) currentSoln = SolveStepByStep[#, primaryunknowns, secondaryunknowns, parameters, currentConstraints, currentSolutions]& /@ currentSystem[[1]]; (* Flattens from the inside out: {{{a,b},{c,d}}} -> {{a,b},{c,d}}. *) currentSoln = (Sequence @@ # &) /@ currentSoln; currentSoln = currentSoln /. (p_ -> q_) :> (p -> Factor[q]); If[depthFirstSolveDebug || recursiveDebug, Print["The solution from the first equation is:"]; Print[CleanUpFunction[ currentSoln ]]; ]; (* If there is no solution for the `simplest' equation, *) (* it terminates this branch of the search tree. *) If[Length[currentSoln] == 0, (* Begin If A3 *) (* * * * * * * * * * * * * * If TRUE of If A3 * * * * * * * * * * * * * * *) (* Goes back to the next branch. *) While[ Length[Last[theList]] === 1, theList = Drop[theList,-1]; (* If all branches have been searched, return soln. *) If[ Length[theList] === 0, Return[soln] ]; ]; theList = Append[Drop[theList,-1], Rest[theList[[-1]]]]; If[depthFirstSolveDebug, Print["The pruned search tree:"]; Print[CleanUpFunction[ theList ]]; ], (* * * * * * * * * * * * * * If FALSE of If A3 * * * * * * * * * * * * * * *) (* Transposes the solution: {{a,b},{c,d}} -> {{a,c}, {b,d}} *) currentSoln = Transpose[currentSoln]; (* Adds constraints to the list, if there are any. .. && b && d *) currentConstraints = (currentConstraints && (And @@ #))& /@ currentSoln[[2]]; (* Drops constraints. *) currentSoln = currentSoln[[1]]; (* Applies the solution to the rest of the system *) (* and simplifies via a together and a numerator. *) currentSystem = Together[ Internal`DeactivateMessages[ Rest[currentSystem] /. currentSoln, Power::infy, General::indet ] ]; currentSystem = Numerator[currentSystem]; If[depthFirstSolveDebug, Print["The current system is:"]; Print[CleanUpFunction[ currentSystem ]]; ]; (* Breaks the expression into a list and removes nonzeros and numbers. *) currentSystem = Numerator[currentSystem]; (* If any of the equations are satisfied, remove them. *) currentSystem = DeleteCases[currentSystem, {___, 0, ___}, {2}]; (* Reformat for the system and solutions to be passed as the next node. *) theNextNode = MapThread[ {#1, #2, #3}&, {currentSystem, currentSoln, currentConstraints}]; (* Remove any invalid solutions. *) theNextNode = DeleteCases[theNextNode, {a_List,_List,_And} /; (!FreeQ[a, DirectedInfinity] || !FreeQ[a, Indeterminate]), {2}]; If[depthFirstSolveDebug, Print["The current (cleaned up) system is:"]; Print[CleanUpFunction[ currentSystem ]]; Print["The theList argument being passed to the next depthFirstSolve is:"]; Print[CleanUpFunction[ Append[theList,theNextNode ] ]]; ]; (* Takes the new solution, the one equation shorter system and the *) (* constraint and continues to the next step lower node. *) AppendTo[theList, theNextNode]; ]; (* End If A3 *) ]; (* End If A2 *) ]; (* End While A1 *) ] (* ************************************************************************ *) (* 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] ]; (* end of block FactorListAndCleanUp *) (* ************************************************************************ *) (* 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 ]; ]; (* end of module ComplexityLevel *) (* ************************************************************************ *) (* This computes the leaves 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]] ] ]; (* end of block SubComplexityLevel *) (* ************************************************************************ *) (* 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 ]; (* end of block SolveStepByStep 1 *) (* ************************************************************************ *) (* 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 ]; (* end of block SolveStepByStep 2 *) (* ************************************************************************ *) (* Takes the factor terms which have the parameters and pass it on to the *) (* SubSolveStepByStep to be solved along with constraints 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 ]; (* end of block SolveStepByStep 3 *) (* ************************************************************************ *) (* 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 outputted 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 of block SubSolveStepByStep *) (* ************************************************************************ *) End[] (* `Private` context *) SetAttributes[AnalyzeAndSolve, ReadProtected] EndPackage[] Print["Package AnalyzeAndSolveV2.m of March 13, 2010 was successfully loaded."]; (* ************************************************************************ *) (* END NEW ANALYZEANDSOLVE *) Print["Package PDESpecialSolutionsV2.m was successfully loaded."]; Print["Version 2, last updated by Hereman on March 21, 2010."]; Print["Previously updated by Goktas on December 13, 2009"]; Print["and by Hereman on January 8 and January 23, 2009."]; Print["Based on code written by Baldwin et al, February 28, 2005."]; (* ************************************************************************ *)