(* FINAL VERSION *) (* ************************************************************************ *) (* : Title : PDESpecialSolutions.m * : Authors : Douglas Baldwin, Unal Goktas, and Willy Hereman. * Department of Mathematical and Computer Sciences * Colorado School of Mines * Golden, Colorado 80401-1887, U.S.A * Contact email: whereman@mines.edu * : Last updated : 12-20-2001 at 21:00 by 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. \[NewLine]PDESpecialSolutions takes the option Form with the default value Tanh. Other options include Tanh, Sech, SechTanh, Cn, and Sn. \[NewLine]PDESpecialSolutions also takes the options NumericTest and SymbolicTest with the default values False. If NumericTest is set to True, PDESpecialSolutions tests the possible solutions based on 13 different sets of random numbers ranging from zero to one for all parameters and remaining constants. Solutions are accepted if they pass one or more of these tests. \[NewLine]If SymbolicTest is set to True the solutions are tested truly symbolically and the result of the test is shown in factored form. \[NewLine]In addition, PDESpecialSolutions has the option InputForm with the default value True. If InputForm is set to False, the output of PDESpecialSolutions will be in standard Mathematica output form." 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`." PDESpecialSolutionsmSolver::noSoln = "mSolver did not find any positive integer values for the dominate behavior, and will now test the values from 1 to 3 to see if the package missed a solution due to inequalities." PDESpecialSolutionsmSolver::fail = "The algorithm failed while attempting to find the positive integer values of m[n]." PDESpecialSolutions::solve = "The algorithm failed while attempting to find the coefficients for the polynomial solution." (* Options definition. *) Options[PDESpecialSolutions] = {Verbose -> False, Form -> tanh, NumericTest -> False, SymbolicTest -> False, InputForm -> False, DegreeOfThePolynomial -> {}}; sech=.; sechtanh=.; cn=.; sn=.; jacobisn=.; jacobicn=.; Begin["`Private`"] (* If called with non-lists, makes them lists. *) PDESpecialSolutions[eqns_?(!ListQ[#]&), funcs_, vars_, param_, opts___?OptionQ]:= PDESpecialSolutions[{eqns}, funcs, vars, param, opts] PDESpecialSolutions[eqns_, funcs_?(!ListQ[#]&), vars_, param_, opts___?OptionQ]:= PDESpecialSolutions[eqns, {funcs}, vars, param, opts] PDESpecialSolutions[eqns_, funcs_, vars_?(!ListQ[#]&), param_, opts___?OptionQ]:= PDESpecialSolutions[eqns, funcs, {vars}, param, opts] PDESpecialSolutions[eqns_, funcs_, vars_, param_?(!ListQ[#]&), opts___?OptionQ]:= PDESpecialSolutions[eqns, funcs, vars, {param}, opts] (* ************************************************************************ *) (* Start of function PDESpecialSolutions. *) PDESpecialSolutions[eqns_List, funcs_List, vars_List, param_List, opts___?OptionQ]:= Block[ {testarg = FreeQ[MapAll[Expand, eqns], (Power[b__, a__] /; (!FreeQ[a, _Symbol] | MemberQ[b, _Real | _Integer]))], theForm = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn -> sn, jacobicn -> cn}, VerboseTest = Verbose /. {opts} /. Options[PDESpecialSolutions], time = TimeUsed[], system, mSoln, newSystem, solution }, (* Protected Local Variables *) If[!testarg, Message[PDESpecialSolutions::poly]]; ( (* If verbose, prints system. *) If[VerboseTest, Print["The given system of differential equations is: "]; Print[ TableForm[eqns /. {Derivative[n__][F_][x__]:> SequenceForm[F, Subscript[SequenceForm @@ Flatten[ Table[#[[1]], {#[[2]]} ]& /@ Transpose[{{x}, {n}}] ] ] ], Sequence @@ ((#->Head[#])& /@ funcs) } ] ] ]; (* Step 1 in the paper. *) (* Transforms the PDE into a nonlinear ODE in $T=Tanh$ or $S=Sech$ *) If[VerboseTest, Print["Transform the differential equation(s) into a nonlinear ODE in "<> ToString[ (theForm /. {tanh -> "T=Tanh", sech -> "S=Sech", sechtanh -> "S=Sech", cn -> "CN=JacobiCN", sn -> "SN=JacobiSN"} ) ] ] ]; system = PDESpecialSolutionsVarChange[eqns /. (a__==b__):>(a-b), funcs, vars, opts]; If[VerboseTest, Print[PDESpecialSolutionsCleanUp[system /. myTrackingVariable[_]->1]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* Step 2 in the paper. *) (* Determines the maximal degree of the polynomial solutions. *) If[VerboseTest, Print["Determine the maximal degree of the polynomial solutions."] ]; mSoln = Internal`DeactivateMessages[ PDESpecialSolutionsmSolver[system, opts], Solve::svars ]; If[Length[mSoln] === 0, Message[PDESpecialSolutionsmSolver::fail];Abort[] ]; (* WILLY 12/19/2001 *) If[VerboseTest, Print[ PDESpecialSolutionsCleanUp[ DeleteCases[mSoln, max[_]->_, 2] ] ]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* Step 3 in the paper. *) (* Derives the algebraic system for the coefficients $a_{ij}, b_{ij}$. *) If[VerboseTest, Print["Derive the nonlinear algebraic system for the coefficients."] ]; newSystem = PDESpecialSolutionsBuildSystem[mSoln, system, vars, param, opts]; If[VerboseTest, Print["The nonlinear algebraic system is"]; Print[PDESpecialSolutionsCleanUp[newSystem /. myTrackingVariable[_]->1]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* Step 4 in the paper. *) (* Solve the nonlinear parameterized algebraic system. *) If[VerboseTest, Print["Solve the nonlinear algebraic system."] ]; solution = Algebra`AnalyzeAndSolve`AnalyzeAndSolve @@ #& /@ newSystem; If[Length[Flatten[solution]] === 0, Message[PDESpecialSolutions::solve];Abort[] ]; If[VerboseTest, Print[PDESpecialSolutionsCleanUp[solution]]; Print["Time Used:", TimeUsed[]-time] ]; time = TimeUsed[]; (* Step 5 in the paper. *) (* Build and test the solitary wave solutions. *) If[VerboseTest, If[( !(NumericTest /. {opts} /. Options[PDESpecialSolutions]) && !(SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), Print["Build the travelling wave solutions."], Print["Build and (numerically and/or symbolically) test the solutions."]; ] ]; solution = PDESpecialSolutionsBuildSolutions[solution, mSoln, eqns, funcs, vars, param, opts]; If[VerboseTest, Print["Time Used:", TimeUsed[]-time]; Print[]; If[( (NumericTest /. {opts} /. Options[PDESpecialSolutions]) || (SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), Print["The (numerically and/or symbolically) tested final solutions:"], Print["The list of possible solutions is:"] ]; ]; (* Returns solutions to user. *) Return[ If[InputForm /. {opts} /. Options[PDESpecialSolutions], InputForm[PDESpecialSolutionsCleanUp[solution]], PDESpecialSolutionsCleanUp[solution] ] ] ) /; testarg ] (* ************************************************************************ *) (* : Title : PDESpecialSolutionsVarChange *) (* : Author : Douglas Baldwin *) (* : Input : The equation list (F[i][x,y,z,t]) and the form tanh. *) (* : Output : The system in the form u[i][T] depending on whether it is tanh. *) (* : Summary : Converts from real space-time domain to u[i][T]. *) PDESpecialSolutionsVarChange[eqns_List, funcs_List, vars_List, opts___?OptionQ]:= Block[{Form, (* The form of the wanted solution. *) funcRules, (* conversion from user functions *) i=0, (* Iterator for myTrackingVariables *) eqList, (* The list of equations, in a usable form *) system, (* The system after transformation to ODE in T *) phi (* A temporary function. *) }, (* Protected Local Variables *) (* Sets the Form from user options *) Form = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn -> sn, jacobicn -> cn}; (* If form isn't correct, it sends a message and aborts. *) If[!(Or @@ ((Form === #)& /@ {tanh, sech, sechtanh, cn, sn})), Message[PDESpecialSolutions::form]; Abort[] ]; (* Creates a set of rules to converts the users functions *) (* (e.g., u[n,t]) to varChangeFunction[i][n,t]. *) funcRules = Table[ Head[funcs[[i]] ] -> varChangeFunction[Form][i], {i, Length[funcs]} ]; (* Adds tacking variables which will be used in the mSolve *) (* function to check for false balances coming form a single *) (* term. *) eqList = If[Head[#] === Plus, Plus @@ ((myTrackingVariable[++i]*#) & /@ List @@ #), myTrackingVariable[++i]*#] & /@ Expand[eqns]; (* Converts the user functions to varChangeFunction[i]. *) eqList = eqList //. funcRules; (* Changes the space of the functions. *) system = eqList /. { varChangeFunction[form_][n_][var__] :> phi[form][n][ Sum[c[i]*{var}[[i]], {i, Length[{var}]} ] ], Derivative[temp__][varChangeFunction[form_][n_]][var__] :> ( D[ phi[form][n][ Sum[c[i]*{var}[[i]], {i, Length[{var}]} ] ], Sequence @@ Transpose[{{var}, {temp}}] ] ) }; (* Parameterizes the sequence $c_1 x+c_2 y+\dots\to\eta$ *) system = system /. { phi[form_][n_][__] :> phi[form][n][eta], Derivative[i_][phi[form_][n_]][__] :> Derivative[i][phi[form][n]][eta] }; (* Repeatedly applies the chain rule. *) system = system /. {Derivative[a_][phi[tanh][n_]][eta] :> Expand[ Nest[(1-T^2)*D[#, T]&, u[n][T], a ] ], phi[tanh][n_][eta] :> u[n][T], Derivative[a_][phi[sech][n_]][eta] :> Expand[ Nest[-(S*Sqrt[1-S^2])*D[#, S]&, u[n][S], a ] ], phi[sech][n_][eta] :> u[n][S], Derivative[a_][phi[sechtanh][n_]][eta] :> Expand[ Nest[-(S*Sqrt[1-S^2])*D[#, S]&, u[n][S], a ] ], phi[sechtanh][n_][eta] :> u[n][S], Derivative[a_][phi[cn][n_]][eta] :> Expand[ Nest[-Sqrt[1-mod-CN^2+2*mod*CN^2-mod*CN^4]*D[#, CN]&, u[n][CN], a ] ], phi[cn][n_][eta] :> u[n][CN], Derivative[a_][phi[sn][n_]][eta] :> Expand[ Nest[Sqrt[1-SN^2-mod*SN^2+mod*SN^4]*D[#, SN]&, u[n][SN], a ] ], phi[sn][n_][eta] :> u[n][SN] }; system = MapAll[Expand, system]; (* Simplifies the system by removing extra (1-T^2)'s *) (* and T's from the system. *) system = If[(# /. T->-1) === 0 && (# /. T->1) === 0, Factor[#/(1-T^2)], Factor[#]]& /@ system; system = If[(# /. T -> 0) === 0, Factor[#/T], Factor[#]]& /@ system; system = If[(# /. S -> 0) === 0, Factor[#/S], Factor[#]]& /@ system; system = If[(# /. S -> -1) === 0 && (# /. S->1) === 0, Factor[#/Sqrt[1-S^2]], Factor[#]]& /@ system; system = If[(# /. CN -> 0) === 0, Factor[#/CN], Factor[#]]& /@ system; system = If[(((# /. CN -> -1) === 0 && (# /. CN -> 1) === 0 ) && (Factor[# /. CN -> Sqrt[-1+mod]/Sqrt[mod] ]) === 0), Factor[#/Sqrt[1-mod-CN^2+2*mod*CN^2-mod*CN^4]], Factor[#]]& /@ system; system = If[(# /. SN -> 0) === 0, Factor[#/SN], Factor[#]]& /@ system; system = If[(((# /. SN -> -1) === 0 && (# /. SN -> 1) === 0 ) && (Factor[# /. SN -> 1/Sqrt[mod] ]) === 0 ), Factor[#/Sqrt[1-SN^2-mod*SN^2+mod*SN^4]], Factor[#]]& /@ system; (* Returns the system back to the PDESpecialSolutionsr *) Return[MapAll[Factor, system]]; ]; (* ************************************************************************ *) (* : Title : PDESpecialSolutionsmSolver *) (* : Author : Douglas Baldwin *) (* : 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, opts___?OptionQ] := Block[{numberOfEquations, Form, myTrackingVariableMax, sqrtTerm, funcRules, mySystem, eqnList0, myMNList, replaceRules, eqnList, mList0, mList, mListStructure, perm, comboList, time = TimeUsed[], mRules, comboList2, mSoln2, freeVariable, freeVariableMax, mSolnFree, mSoln, mSoln0, checkList2}, (* Protects Local Variables *) (* Sets the form from the options. *) Form = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn->sn, jacobicn->cn}; (* Determines the number of myTrackingVariables *) myTrackingVariableMax = 0; system /. myTrackingVariable[a_Integer] :> myTrackingVariable[ If[a > myTrackingVariableMax, myTrackingVariableMax++]; a]; (* 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 }; (* Sets up a list of rules for replacing the form to F. *) funcRules = {tanh -> T, sech -> S, sechtanh -> S, cn -> CN, sn -> SN}; (* Breaks up system into most general form. *) mySystem = {Coefficient[#,Sqrt[Form /. sqrtRules],0]& /@ system, Coefficient[#,Sqrt[Form /. sqrtRules],1]& /@ system}; (* Simplifies the system by removing extra (1-T^2)'s *) (* and T's from the system. *) mySystem = If[(# /. T->-1) === 0 && (# /. T->1) === 0, Factor[#/(1-T^2)], Factor[#]]& /@ #& /@ mySystem; mySystem = If[(# /. T->0) === 0, Factor[#/T], Factor[#]]& /@ #& /@ mySystem; mySystem = If[(# /. S->0) === 0, Factor[#/S], Factor[#]]& /@ #& /@ mySystem; mySystem = If[(# /. S->-1) === 0 && (# /. S->1) === 0, Factor[#/Sqrt[1-S^2]], Factor[#]]& /@ #& /@ mySystem; mySystem = If[(# /. CN -> 0) === 0, Factor[#/CN], Factor[#]]& /@ mySystem; mySystem = If[(((# /. CN -> -1) === 0 && (# /. CN -> 1) === 0 ) && (Factor[# /. CN -> Sqrt[(-1+mod)/mod] ]) === 0 ), Factor[#/Sqrt[1-mod-CN^2+2*mod*CN^2-mod*CN^4]], Factor[#]]& /@ mySystem; mySystem = If[(# /. SN->0) === 0, Factor[#/SN], Factor[#]]& /@ mySystem; mySystem = 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[#]]& /@ mySystem; (* Replaces the function u[i_] with a function in *) (* F^m[i] or F^n[i]. *) eqnList0 = Expand[ DeleteCases[ Flatten[ {mySystem[[1]] /. u[i_] :> Function[{F}, F^m[i] ], mySystem[[2]] /. u[j_] :> Function[{F}, F^n[j] ]} ], 0, Infinity ] ]; (* For when there is a Sqrt[...] term, we use m and n *) myMNList = {}; eqnList0 /. {m[i_Integer] :> (myMNList = Append[myMNList, m[i]];m[i]), n[i_Integer] :> (myMNList = Append[myMNList, n[i]];n[i])}; myMNList = Union[myMNList]; (* Sets numberOfEquations to the length number of equations. *) numberOfEquations = Length[myMNList]; (* Breaks up expressions into lists of terms. *) eqnList = If[Head[#]===Plus, List @@ #, {#}]& /@ eqnList0; (* Pulls off the exponents of T and forms a list *) (* of expressions of the form {{{1+m[1]},{..}..}..} *) mList0 = Exponent[#, Form /. funcRules ]& /@ #& /@ eqnList; (* 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+\dotsc$ *) (* where $a,b,c,\dotsc,i,j,k,\dotsc\in\mathbb{R}$ into lists. *) mList = Map[ If[Head[#]===Plus, List @@ #, If[(Head[#]===m || Head[#]===n) || Head[#]===Times, {#}, # ] ]&, mList0, 2 ]; (* The following routine strips off the $a$ in the above *) (* expression and leaves only the underlying structure. *) mListStructure = Union[# /. {a_Integer, b__}:>{b}& /@ #]& /@ mList; (* Re-organizes the list of powers of T by the structure *) (* listed above. *) mList = Table[ Flatten[ Table[ Cases[#, Flatten[{_, mListStructure[[i,j]]}] | mListStructure[[i,j]] ]& /@ {mList[[i]]}, {j, Length[mListStructure[[i]] ]} ], 1 ], {i, Length[mList0]} ]; (* Determines the maximum $a$ in each power of T *) mList = {Max[# /. {a_, ___}:>If[IntegerQ[a], a, 0]& /@ #]}& /@ #& /@ 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+\dotsb+d\,m_i\,m_j$ *) mList = Plus @@ #& /@ #& /@ Table[ Table[ Join[mList[[i,j]], mListStructure[[i,j]] ], {j, Length[mList[[i]] ]} ], {i, Length[mList]} ]; (* For the purpose of the initial combinatorics, we must *) (* pad lists of only one element. *) mList = Table[ If[Length[ mList[[i]] ]===1, PadRight[mList[[i]], 2], mList[[i]] ], {i, Length[mList]} ]; (* A little delayed equal to simplify the expression below. *) perm := Permutations[Range[Length[mList[[j]] ] ] ]; (* Combines the expressions for the powers of F in order to *) (* find potential balances to the system. *) comboList = Flatten[ Table[ Table[ Partition[ Flatten[ Union[ Table[ Prepend[ Sort[ Drop[perm[[i]], -Length[ perm[[1]] ] + k] ], j ], {i, 1, Length[perm]} ] ] ], k+1 ], {k, 2, Length[mList[[j]] ]} ], {j, Length[mList]} ], 2 ]; (* Combines the polynomials in $m_i$ and solves for $m_i$ *) mRules = Table[ Union[ Flatten[ Table[ Solve[Equal @@ Flatten[ Table[ mList[[comboList[[i, 1]], comboList[[i, k]] ]], {k, 2, Length[comboList[[i]] ]} ] ], myMNList[[j]] ], {i, 1, Length[comboList]} ] ] ], {j, 1, numberOfEquations} ]; (* Prints out a warning if it's taking too long. *) If[(TimeUsed[]-time>3), Print["Still computing the maximal degree, please wait."] ]; (* Creates a set of rules for combining the *) (* above solutions of $m_i$ *) comboList2 = Partition[ Flatten[ Fold[Table, Array[beta, numberOfEquations], Table[{beta[i], Length[mRules[[i]] ]}, {i,numberOfEquations} ] ] ], numberOfEquations ]; (* Takes these solutions and uses them to find *) (* actual integer solutions to $m_i$ *) mSoln = Table[ Solve[ Table[Equal @@ mRules[[i, comboList2[[j, i]] ]], {i, numberOfEquations} ], myMNList ], {j, Length[comboList2]} ]; (* If there is a user case, it adds it to the computed *) (* polynomial solutions. *) If[(DegreeOfThePolynomial /. {opts} /. Options[PDESpecialSolutions]) =!= {}, mSoln = Append[mSoln, ToExpression[ StringReplace[ ToString[DegreeOfThePolynomial /. {opts}], {"m" -> "Calculus`PDESpecialSolutions`Private`m", "n" -> "Calculus`PDESpecialSolutions`Private`n" } ] ] ] ]; (* Cleans up the list by removing duplicates *) (* and extra braces *) mSoln = Union[ Table[ Flatten[mSoln[[i]] ], {i, Length[mSoln]} ] ]; (* Keeps only the solutions which are numeric. *) (* Id est, removes any complex solutions. *) mSoln2 = Complement[ mSoln, Flatten[ Union[ Table[ If[NumericQ[Evaluate[m[j] /. mSoln[[i]] ]], mSoln[[i]], Null ], {i, Length[mSoln]}, {j, numberOfEquations} ] ], 1 ] ]; (* Same as above, but removes all non-integer solutions. *) mSoln = Complement[ mSoln, Flatten[ Union[ Table[ If[IntegerQ[ Evaluate[m[j] /. mSoln[[i]] ]], Null, mSoln[[i]] ], {i, Length[mSoln]}, {j, numberOfEquations} ] ], 1 ] ]; (* Finally, it removes all solutions which are not positive. *) mSoln = Complement[mSoln, Flatten[Union[ Table[ If[Positive[ Evaluate[m[j] /. mSoln[[i]] ]], Null, mSoln[[i]] ], {i, Length[mSoln]}, {j, numberOfEquations} ] ], 1 ] ]; (* Sets up a local variable with these solutions. *) mSoln0 = mSoln; (* Pulls off the polynomials in m_i which are the *) (* dominant terms. *) checkList2= Table[ DeleteCases[ Union[ Flatten[ Table[ Table[ If[ Evaluate[mList0[[i,k]] /. mSoln[[j]] ]== Max[mList0[[i]]/.mSoln[[j]]], mList0[[i,k]] ], {k, Length[mList0[[i]]]} ], {j, Length[mSoln]} ] ] ], Null ], {i,numberOfEquations} ]; (* Checks to see if any of the highest powers are *) (* based on a term which is based on only one m_i *) (* If this is the case, it may self balance. *) freeVariable = DeleteCases[ Flatten[ Table[ If[ Length[ Cases[checkList2[[i]], ___Integer+__*m[j]|___Integer+m[j] ] ]>1 && numberOfEquations>1, m[j], Null ], {j, numberOfEquations}, {i, numberOfEquations} ] ], Null ]; (* The following code runs only if an inequality *) (* exists and solutions might have been missed. *) (* To find the missing solutions, it does a count- *) (* down scheme and then test them to see if any *) (* of them are valid choices. *) If[Length[freeVariable]=!=0, freeVariableMax = Max[ Table[freeVariable[[j]] /. mSoln[[i]], {i, Length[mSoln]}, {j, Length[freeVariable]} ] ]; mSolnFree = Partition[ Flatten[ Fold[Table, Table[m[i]->beta[i], {i, numberOfEquations}], Table[{beta[k],1,freeVariableMax},{k,numberOfEquations}] ] ], numberOfEquations ]; mSolnFree = Complement[ DeleteCases[ Table[ If[Plus @@ Flatten[ Table[ If[ Length[ Cases[mList0[[i]] /. mSolnFree[[j]], Max[mList0[[i]] /. mSolnFree[[j]]] ] ]>=2, 1, 0 ], {i, numberOfEquations} ] ]>=numberOfEquations, mSolnFree[[j]], Null ], {j, Length[mSolnFree]} ], Null ], mSoln ] ]; (* Displays message to user if freedom exists. *) Do[ If[ Length[ Cases[checkList2[[i]], ___Integer+__*m[j]|___Integer+m[j] ] ]>1 && numberOfEquations>1, Message[PDESpecialSolutionsmSolver::freedom, PDESpecialSolutionsCleanUp[Cases[checkList2[[i]], ___Integer+__*m[j]|___Integer+m[j] ] ], i ] ], {j, numberOfEquations}, {i, numberOfEquations} ]; If[Length[mSolnFree]=!=0, mSoln = Join[mSoln, mSolnFree] ]; (* The below function takes advantage of the tracking *) (* variables. Removes false solutions from mSoln, *) (* by checking balances of unique terms. *) mSoln = DeleteCases[ Table[ If[Plus @@ Flatten[ Table[ If[ Length[ Cases[#, Max[#] ]&[ Table[ Exponent[ Coefficient[ Expand[eqnList0[[i]] ], myTrackingVariable[k]] /. mSoln[[j]], T ], {k,1,myTrackingVariableMax} ] ] ]>=2, 1, 0 ], {i, numberOfEquations} ] ]>=numberOfEquations, mSoln[[j]], Null ], {j, Length[mSoln]} ], Null ]; (* If the form is of SechTanh, it uses the solution m_i=2, n_i=1 *) If[Form === sechtanh, mSoln = { Flatten[ Table[{m[i]->2, n[i]->1}, {i, Length[system]} ] ] } ]; (* If it doesn't find any find any solutions, it quits. *) If[Length[mSoln] == 0, Message[PDESpecialSolutionsmSolver::fail]; Abort[ ] ]; (* Appends the maximum value of m_i. *) mSoln = Table[ Join[mSoln[[j]], Table[max[i]->(Max[mList0[[i]] /. mSoln[[j]] ]), {i, Length[system]} ] ], {j, Length[mSoln]} ]; (* Prints out a warning if it's taking too long. *) If[ (TimeUsed[]-time>3) && !(Verbose /. {opts} /. Options[PDESpecialSolutions]), Print["Still computing the maximal degree, please wait."] ]; (* Returns the solutions. *) Return[mSoln]; ]; (* ************************************************************************ *) (* : Title : PDESpecialSolutionsBuildSystem. * : Author : Douglas Baldwin * : Date : 05-24-01 *) PDESpecialSolutionsBuildSystem[mSoln_List, system_List, vars_List, param_List, opts___?OptionQ]:= Block[{Form, numberOfEquations, sqrtRules, funcRules, uRules, tempuRules, newSystem, waveparameters, unknowns, nonzeros, maxF, 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, j}] + If[(Form /. sqrtRules) =!= 0, Sqrt[Form /. sqrtRules]* Sum[b[i,kk]*(Form /. funcRules)^kk, {kk, 0, n[i] /. Flatten[#]} ], 0 ] ] ) )& /@ mSoln; (* WILLY 12/19/2001 modification: introduced tempuRules *) If[Verbose /. {opts} /. Options[PDESpecialSolutions], Print["Seeking polynomial solutions of the form"]; tempuRules = DeleteCases[uRules, n[_]->_, 2]; Print[ PDESpecialSolutionsCleanUp[ DeleteCases[tempuRules, max[_]->_, 2 ] ] ] ]; (* Prints out a warning, if it's taking a long time. *) If[TimeUsed[]-time>3, Print["Still building the nonlinear algebraic system, please wait."] ]; (* Converts form of solutions given by solver to a pure *) (* function. *) toPure = (# /. (a__[var__]->temp__):> (a:>Function[{var}, temp]))&; (* Applies the expressions in T to the system and removes the *) (* tracking variables. *) newSystem = Expand[ system //. Append[toPure[#], myTrackingVariable[_]->1] ]& /@ uRules; (* Clears the denominator to simplify the system. *) newSystem = Map[Together[#*Denominator[Together[#]]]&, newSystem, 2]; (* Creates a list with the wave parameters (c[1], c[2], ...) to be *) (* passed to the solver. *) waveparameters = Join[ Table[c[i], {i, Length[vars]}], 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, 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. *) 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 = #[[2]]& /@ #& /@ (Sort[{# /. a_[b_,c_]:>{b,a,c},#}& /@ Flatten[#] ]& /@ unknowns); (* Expands the system. *) newSystem = MapAll[Expand,newSystem]; (* Determines the maximum value of F in each of the cases. *) maxF = Map[ Exponent[#, Form /. funcRules ]&, newSystem ]; (* Brakes the expressions up newSystem by powers of T. *) newSystem = Table[ Table[ (Sort[{Length[Expand[#]], #}& /@ Union[ Flatten[ DeleteCases[ Flatten[ Table[ Table[ Coefficient[ Expand[newSystem[[case,k]] ], Form /. funcRules, i ], {i, 0, maxF[[case,j]]}], {j, Length[maxF[[case]] ]}] ], 0 ] ] ] ] /. {_Integer, a__}:>a ), {k, Length[newSystem[[case]] ]} ], {case, Length[newSystem]}]; (* 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 }; (* Breaks up system into most general form. *) newSystem = Flatten[ {Coefficient[#,Sqrt[Form /. sqrtRules],0]& /@ #, Coefficient[#,Sqrt[Form /. sqrtRules],1]& /@ #} ]& /@ #& /@ newSystem; (* Converts from expressions to equations. *) newSystem = DeleteCases[(# == 0)& /@ Flatten[#]& /@ newSystem, False || True, Infinity]; (* 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."] ]; (* Combines into lists containing {newSystem, unknowns, *) (* waveparameters, param, nonzeros} for the solver. *) Return[ MapAll[Factor, Transpose[ {newSystem, unknowns, Table[waveparameters, {Length[newSystem]}], Table[param, {Length[newSystem]}], nonzeros } ] ] ] ]; (* ************************************************************************ *) (* : 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[{Form, (* Form defined by user. *) solutionRules, (* Local version of solution *) sqrtRules, (* Set by options. *) funcRules, (* Also set by options. *) uRules, (* Rules for creating the solutions. *) mySimplificationRules, (* Simplification rules *) numericTestSolutions, (* tests validity of solutions *) symbolicTestSolutions, (* Also, tests the validity of soln. *) solns = {}, (* confirmed solutions. *) warn = {}, (* Collects a list of warnings due to simplification. *) time = TimeUsed[] (* A timer for length warnings. *) }, (* Protected Local Variables. *) (* Sets the form from the options. *) Form = ToExpression[ ToLowerCase[ ToString[Form /. {opts} /. Options[PDESpecialSolutions]] ] ] /. {jacobisn -> sn, jacobicn -> cn}; (* Sets up a local version to modify *) solutionRules = MapAll[Factor, solution]; (* Rules for building the solutions. *) sqrtRules = {tanh -> 0, sech -> 0, sechtanh -> 1-Sech[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}] + phase]^2, cn -> 0, sn -> 0 }; funcRules = {tanh -> Tanh[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase], sech -> Sech[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase], sechtanh -> Sech[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase], cn -> CN[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase,mod], sn -> SN[Plus @@ Table[c[i]*vars[[i]], {i, Length[vars]}]+phase,mod] }; (* Converts the results of mSolve to an expression in T *) uRules = (m[i_]->j_) :> (funcs[[i]] -> Sum[a[i,kk]*(Form /. funcRules)^kk, {kk, 0, 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[ 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]}], max[_]->_, Infinity ], n[_]->_, Infinity], Rule[a_, b_] /; a == b && FreeQ[a, Alternatives @@ funcs], Infinity]; (* Simplification rules for reducing the solution *) (* even further . . . so as to find more general *) (* solutions. *) mySimplificationRules = { Sqrt[xxx__^2] :> (warn = Append[warn, "Sqrt[a^2]->a"]; xxx), Sqrt[Power[yyy__,2]] :> (warn = Append[warn, "Sqrt[a^2]->a"]; yyy), Sqrt[-zzz__^2] :> (warn = Append[warn, "Sqrt[-a^2]->I*a"]; I*zzz), Sqrt[-ttt_*zzz__^2] :> (warn = Append[warn, "Sqrt[-a*b^2]->I*b*Sqrt[a]"]; I*zzz*Sqrt[ttt]), (1-Sech[xx__]^2) -> Tanh[xx]^2, (1-Tanh[xx__]^2) -> Sech[xx]^2, JacobiDN[xx__, mm__]^2 :> 1-mm+mm*JacobiCN[xx,mm]^2 /; Form === cn, JacobiSN[xx__, mm__]^2 :> 1-JacobiCN[xx,mm]^2 /; Form === cn, JacobiDN[xx__, mm__] :> (warn = Append[warn,"JacobiDN[x,mod]->Sqrt[1-mod+mod*JacobiCN[x,mod]^2]"]; Sqrt[1-mm+mm*JacobiCN[xx,mm]^2] ) /; Form === cn, JacobiSN[xx__, mm__] :> (warn = Append[warn,"JacobiSN[x,mod]->Sqrt[1-JacobiCN[x,mod]^2]"]; Sqrt[1-JacobiCN[xx,mm]^2] ) /; Form === cn, JacobiDN[xx__, mm__]^2 :>1-mm*JacobiSN[xx,mm]^2 /; Form === sn, JacobiCN[xx__, mm__]^2 :> 1-JacobiSN[xx,mm]^2 /; Form === sn, JacobiDN[xx__, mm__] :> (warn = Append[warn, "JacobiDN[x,mod]->Sqrt[1-mod*JacobiSN[x,mod]^2]" ]; Sqrt[1-mm*JacobiSN[xx,mm]^2] ) /; Form === sn, JacobiCN[xx__, mm__] :> (warn = Append[warn, "JacobiSN[x,mod]->Sqrt[1-JacobiSN[x,mod]^2]]" ]; Sqrt[1-JacobiSN[xx,mm]^2] ) /; Form === sn }; (* Applies the above rules to the solutions. *) solutionRules = FixedPoint[ MapAll[Expand, MapAll[Factor, # //. mySimplificationRules ] //. mySimplificationRules ]&, solutionRules, 3 ]; (* Removes Trivial Solutions *) solutionRules = Select[#, !(And @@ (FreeQ[# /. (a_->b__):>b, Alternatives @@ vars]& /@ #))& ]& /@ solutionRules; If[Verbose /. {opts} /. Options[PDESpecialSolutionsr], Print["The possible non-trivial solutions (before any testing) are: "]; Print[ PDESpecialSolutionsCleanUp[solutionRules /. {CN -> JacobiCN, SN -> JacobiSN} ]] ]; (* Prints out a warning, if it's taking a long time. *) If[TimeUsed[]-time>3, Print["Still building the exact solutions, please wait."] ]; (* Depending on the NumericTest option, it either tests the *) (* solutions numerically, or it doesn't. *) If[ ( !(NumericTest /. {opts} /. Options[PDESpecialSolutions]) && !(SymbolicTest /. {opts} /. Options[PDESpecialSolutions]) ), CellPrint[ Cell[ "These solutions are not being tested numerically or symbolically. " <> "To test the solutions set either the NumericTest option " <> "to True, or set the SymbolicTest option to True, or both. ", "Message" ] ]; If[Length[warn]>0, CellPrint[ Cell[ "The following simplification rules are being used: " <> ToString[Union[warn]], "Message" ] ] ]; Return[ Map[Union, MapAll[Factor, MapAll[Expand, solutionRules /. {CN -> JacobiCN, SN -> JacobiSN} ] ], 2 ] ] ]; (* Converts form of solutions given by solver to a pure function. *) toPure = (# /. (a__[var__]->temp__):> (a:>Function[{var}, temp]))&; (* If statement for the numeric testing option. *) If[NumericTest /. {opts} /. Options[PDESpecialSolutions], (* Sends message to user. *) If[!(VerboseTest /. {opts} /. Options[DDESolve]), Print["Numerically testing the solutions."] ]; (* Tests numerically to make sure the above solutions are valid. *) numericalTestSolutions = {(eqns /. a__==b__:>a-b) /. toPure[MapAll[TrigToExp, #]], #}& /@ #& /@ solutionRules; (* Numerically tests the solutions by replacing all the symbols with *) (* random numbers 13 separate times. *) numericalTestSolutions = {Plus @@ (Table[ (* Prints out a warning, if it's taking a very long time. *) If[TimeUsed[]-time>3, time = TimeUsed[]+1; Print["Still numerically testing the solutions "<> "numerically, please wait." ] ]; And @@ ( ( Abs[# /. {a[_,_]->Random[], (* coefficients in polynomial *) b[_,_]->Random[], (* coefficients in polynomial *) Sequence @@ ((# -> Random[])& /@ If[Form === cn || Form === sn, Join[param, {mod}], param ] ) } /. {CN -> JacobiCN, SN-> JacobiSN}] < 10^(-10) )& /@ MapAll[ Expand, #[[1]] /. (Append[ (# -> Random[])& /@ Join[vars, Array[c, Length[vars]]], phase -> 0 ] ) ] ), {13} ] /. {True->1, False->0} ), # }& /@ #& /@ numericalTestSolutions; (* If it works for any of the 13 times, it keeps the solution. *) (* This works out to be about a 0.0001 chance of missing *) (* a correct solution *) solns = Cases[numericalTestSolutions, {a_Integer/;a>0, _List}, Infinity ] //. {_Integer, a_List}:>{a[[2]]}; ]; (* The second testing If statement. *) If[SymbolicTest /. {opts} /. Options[PDESpecialSolutions], (* Sends message to user. *) If[!(VerboseTest /. {opts} /. Options[DDESolve]), Print["Symbolically testing the solutions."] ]; (* This sets up the input for the next block of code. *) (* In that, it takes the question given by the user, *) (* and replaces the user defined functions (e.g. u[x,t]) *) (* with the functions found with the algorithm. To replace *) (* the partial derivatives correctly, we must first convert *) (* the solutions to pure functions. *) symbolicTestSolutions = {(eqns /. a__==b__:>a-b) /. toPure[MapAll[TrigToExp, #]], #}& /@ #& /@ (solutionRules /. {CN -> JacobiCN, SN-> JacobiSN}); (* Attempts to simplify the expressions as much as possible. *) symbolicTestSolutions = ExpToTrig[ FixedPoint[ MapAll[Factor, MapAll[Expand, #] /. mySimplificationRules ] /. mySimplificationRules &, #]& /@ symbolicTestSolutions]; (* Pulls of the solutions which simplify to zero, and adds *) (* them to solns (the final output applies the union, so *) (* duplicate solutions are not an issue at this point). *) solns = Join[solns, Cases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ] //. {{(0)..}, a_List}:>{a} ]; (* Removes all the good cases, so we may output the bad *) (* ones to the user. *) symbolicTestSolutions = DeleteCases[ DeleteCases[symbolicTestSolutions, {{(0)..}, _List}, Infinity ], {} ]; If[(Map[Union, symbolicTestSolutions,3] //. {{}}->{}) != {}, (* Sends out the left over (questionable) solutions to *) (* the user for hand inspection. *) CellPrint[ Cell[ "These solutions did not simplify to zero. " <> "Please test these equations by hand or interactively "<> "with Mathematica.", "Message" ] ]; Print[PDESpecialSolutionsCleanUp[#]]& /@ symbolicTestSolutions; ]; ]; If[Length[warn]>0, CellPrint[ Cell[ "The following simplification rules are being used " <> "to test the solutions: " <> ToString[Union[warn]], "Message" ] ] ]; (* Returns the solution rules *) Return[ Union[ MapAll[Factor, MapAll[Expand, solns /. {CN -> JacobiCN, SN-> JacobiSN} ] ] ] ]; ]; (* ************************************************************************ *) (* : Title : PDESpecialSolutionsCleanUp *) (* : Author : Douglas Baldwin *) (* : Summary : Remove "PDESpecialSolutionsPrivate`" from output. *) PDESpecialSolutionsCleanUp = ToExpression[ StringReplace[ToString[InputForm[#]], "Calculus`PDESpecialSolutions`Private`"->"" ] ]& End[] (* `Private` context *) Protect[PDESpecialSolutions] SetAttributes[PDESpecialSolutions, ReadProtected] EndPackage[] (* ************************************************************************ *) (* Nonlinear algebraic solver by Unal Goktas (WRI) and Willy Hereman *) (* written by Unal Goktas in October/November 2000 *) BeginPackage["Algebra`AnalyzeAndSolve`"] Unprotect[AnalyzeAndSolve] 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[] ] ] ] AnalyzeAndSolve[system: {HoldPattern[_ == _]..}, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_] := AnalyzeAndSolve[(#[[1]]-#[[2]]&) /@ system, primaryunknowns, secondaryunknowns, parameters, nonzeros] AnalyzeAndSolve[system_, primaryunknowns_, secondaryunknowns_,{},nonzeros_]:= Block[{constraints}, constraints = (And @@ ((# != 0 &) /@ nonzeros)); Internal`DeactivateMessages[ Solve[(And @@ ((# == 0 &) /@ system)) && constraints, Join[primaryunknowns, secondaryunknowns], Sort -> False], Solve::svars ] ] (* ************************************************************************ *) (* "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] ] (* ************************************************************************ *) (* RecursivelyKeepTrack[{}, __] := {} *) RecursivelyKeepTrack[{}, sol_, __] := (AppendTo[globalsol, sol]; {}) RecursivelyKeepTrack[system_, sol_, __] /; (!FreeQ[system, DirectedInfinity] || !FreeQ[system, Indeterminate]) := {} RecursivelyKeepTrack[system_, sol_, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_, constraints_] := Block[{a, b, c}, a = FactorListAndCleanUp[#, primaryunknowns, secondaryunknowns, parameters, nonzeros]& /@ system; a = Union[a]; a = Sort[a, (ComplexityLevel[#1, primaryunknowns, secondaryunknowns, parameters] <= ComplexityLevel[#2, primaryunknowns, secondaryunknowns, parameters]&)]; b = a[[1]]; b = SolveStepByStep[#, primaryunknowns, secondaryunknowns, parameters, constraints /. sol, sol]& /@ b; b = (Sequence @@ # &) /@ b; If[Length[b] == 0, Return[{}]]; b = Transpose[b]; c = (constraints && (And @@ #))& /@ b[[2]]; b = b[[1]]; a = Together[ Internal`DeactivateMessages[Rest[a] /. b, Power::infy, General::indet ] ]; a = Numerator[a]; a = DeleteCases[a, {___, 0, ___}, {2}]; MapThread[ RecursivelyKeepTrack[#1, #2, primaryunknowns, secondaryunknowns, parameters, nonzeros, #3]&, {a, b, c}] ] (* ************************************************************************ *) FactorListAndCleanUp[ system_, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_] := Block[{a}, a = Flatten[Map[First[#]&, FactorList /@ system, {2}]]; a = DeleteCases[a, _?(NumericQ[#] || MemberQ[nonzeros, #]&)]; Union[a] ] (* ************************************************************************ *) ComplexityLevel[expr_List, primaryunknowns_,secondaryunknowns_,parameters_]:= Max[ComplexityLevel[#, primaryunknowns, secondaryunknowns, parameters]& /@ expr] ComplexityLevel[expr_, primaryunknowns_, secondaryunknowns_, parameters_] := Block[{a, b}, a = Union[Cases[{expr}, q_ /; MemberQ[primaryunknowns, q], -1]]; b = Length[a]; ( SubComplexityLevel[expr, a, b, 1] ) /; b != 0 ] (* ************************************************************************ *) ComplexityLevel[expr_, primaryunknowns_, secondaryunknowns_, parameters_] := Block[{a, b}, a = Union[Cases[{expr}, q_ /; MemberQ[secondaryunknowns, q], -1]]; b = Length[a]; ( SubComplexityLevel[expr, a, b, 2] ) /; b != 0 ] (* ************************************************************************ *) ComplexityLevel[expr_, primaryunknowns_, secondaryunknowns_, parameters_] := Block[{a, b}, a = Union[Cases[{expr}, q_ /; MemberQ[parameters, q], -1]]; b = Length[a]; ( SubComplexityLevel[expr, a, b, 3] ) /; b != 0 ] (* ************************************************************************ *) SubComplexityLevel[expr_, a_, b_, level_] := Block[{c = Select[a, PolynomialQ[expr, #]&]}, Which[ Length[c] == 0, 100*b*level+LeafCount[expr], True, Min[Exponent[expr, c]] ] ] (* ************************************************************************ *) (* solves a single equation *) 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 ] (* ************************************************************************ *) 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 ] (* ************************************************************************ *) 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 ] (* ************************************************************************ *) SolveStepByStep[___] := {} 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 = Reduce[eqn == 0 && constraints, Flatten[{a, pars}], Sort -> False]; a = {ToRules[c]}; If[Length[a] == 0, Return[{}]]; c = Cases[{#}, p_ != q_, -1]& /@ If[Head[c] === Or, List @@ c, {c}]; b = Internal`DeactivateMessages[sol /. a, Power::infy, General::indet]; Union[Transpose[{Flatten /@ Thread[{b, a}, List, 2], c}]] ] (* ************************************************************************ *) End[] (* `Private` context *) SetAttributes[AnalyzeAndSolve, ReadProtected] EndPackage[] Print["Package PDESpecialSolutions.m was succesfully loaded."]; (* ************************************************************************ *)