(* : Title : Painleve Integrabilty Test *) (* : Author : Douglas Baldwin *) (* : Summary : This package test a system of PDEs for the Painleve property. *) (* : Package Version : 1.0 *) (* : Data : August 8, 2000 *) (* : Mathematica Version : 4.0 *) (* : Keywords : Painleve, integrability *) (* : Warnings : None *) (* : Limitations : Only deals with cases with positive resonances, as all negative resonances would require the Taylor series rather than the Laurent series. *) (* : Example : PainleveTest[{D[u[x, t], t] + 6*u[x, t]*D[u[x, t], x] + D[u[x, t], {x, 3}] == 0}, {u[x, t]}, {x, t}, Kruskal->x] *) BeginPackage["Painleve`"] PainleveTest::usage = "PainleveTest[{eqn1, eqn2, ...}, {funcs}, {vars}, opt___] test for the Painleve property of a system of nonlinear ODEs and/or PDEs. \[NewLine]\[NewLine] PainleveTest[] brings up a menu of test cases. \[NewLine]\[NewLine] The Painleve Test is made up of three steps: 1) the dominant behavior is determined, 2) the resonances are computed, and 3) the constants of integration are determined. For the system of ODEs/PDEs to be of the Painleve Property all three steps must passed unambiguously. If freedom or compatibility conditions exist, a message will be sent to the user. \[NewLine]\[NewLine] The optional argument are, i) Kruskal->var, which will use the Kruskal Simplification. For example Kruskal->x, for varaibles {x,y,z,t} would assumes g[x,y,z,t]->x-h[y,z,t] and the partial derivative with respect to x of u[i,j] are zero, ii) WorkLog->\"FileName\", will write to a file in the current directory the equations, the alphas, the resonances, the constants, and the net time. \[NewLine]\[NewLine] Example: PainleveTest[{D[u[x, t], t] + 6*u[x, t]*D[u[x, t], x] + D[u[x, t], {x, 3}] == 0}, {u[x, t]}, {x, t}, Kruskal->x, WorkLog->\"worklog.m\"]" alphaSolver::freedom = "Freedom exists in this system, as `1` are both dominate powers in expression `2`." alphaSolver::fail = "This system failed the Painleve test while determinining the negative integer values of alpha[n]." alphaSolver::noSoln = "Alpha Solver did not find any negative integer values for alpha, and will test the values from -3 to -1 to see if the package missed a solution due to inequalities." initialSolver::fail = "This system failed the Painleve test while determining the expression(s) for u[i,0]." resonanceSolver::fail = "This system failed the Painleve test while determining the resonances of the system." resonanceSolver::fail2 = "This system failed the Painleve test while determining the resonances of the system. As, `1`, are not non-negative integers." PainleveTest::symbol = "One of the equations had a power that was of type Symbol." PainleveTest::exp = "All the terms in the equations must be polynomial." constantsSolver::fail = "This system failed the Painleve test while determining the constants of the system." constantsSolver::true = "This system is free at resonance `1`." Unprotect[alphaSolver, alphaSolve, resonanceSolver, resonanceSolve, constantsSolve, PainleveTest]; (* : Title: pdeForm *) (* : Summary : Displays functions and expressions in standard mathematical notation. *) (* : Authors : Initially done by Tracy Otto and Tony Miller (1995), later updated first by U. Goktas, then by Mike Colagrosso, and finally rewritten by Douglas Baldwin *) pdeForm[expr_] := expr /. {Derivative[n__][u_[a_,b_]][x__] :> SequenceForm[ Subscript[u, SequenceForm["(",a,",",b,")"] ], Subscript[ SequenceForm @@ Flatten[ Table[#[[1]], {#[[2]]} ]& /@ Transpose[{{x}, {n}}] ] ] ], Derivative[n__][F[a_]][x__] :> SequenceForm[ Subscript[u,a], Subscript[ SequenceForm @@ Flatten[ Table[#[[1]], {#[[2]]}]& /@ Transpose[{{x}, {n}}] ] ] ], Derivative[n__][F_[a_]][x__] :> SequenceForm[ Subscript[F,a], Subscript[ SequenceForm @@ Flatten[ Table[#[[1]], {#[[2]]}]& /@ Transpose[{{x}, {n}}] ] ] ], Derivative[n__][g_][x__] :> SequenceForm[g, Subscript[ SequenceForm @@ Flatten[ Table[#[[1]], {#[[2]]}]& /@ Transpose[{{x}, {n}}] ] ] ], F[a_][__] :> Subscript[u,a], g[__] :> g, u[a_,b_][__] :> Subscript[u, SequenceForm["(",a,",",b,")"] ], u[a_,b_] :> Subscript[u, SequenceForm["(",a,",",b,")"] ] }; (* : Title: toPure *) (* : Author : Douglas Baldwin *) (* : Summary : Converts u[i,j][x,y,...]->expr to u[i,j]:>Function[{x,y,z,t,...}, expr]. *) toPure = (# /. (u_[a_,b_][var__]->temp__):> (u[a,b]:>Function[{var}, temp]))&; (* : Title: alphaSolver: *) (* : Summary : Determines the dominant behavior of the system. *) (* : Input : system of equations, the dependent variables, and the independent variables. *) (* : Output : A list of sub-list of alpha solutions, e.g. {{alpha[1]->-2, alpha[2]->-2}, {alpha[1]->-1, alpha[2]->-1}} *) debug = debug1 = debug2 = debug3 = debug4 = False; alphaSolver[eqList_List, funcs_List, vars_List] := Block[{newFuncs = Array[F, Length[funcs]], funcsRules, eqnList0, numberOfEquations, replaceRules, eqnList, alphaList0, alphaList, alphaListStructure, perm, comboList, alphaRules, comboList2, alphaSoln2, freeVariable, freeVariableMin, alphaSolnFree, alphaSoln, alphaSoln0, checkList2}, Off[Solve::svars]; Off[Solve::ifun]; funcsRules = Table[ Head[funcs[[i]]] -> newFuncs[[i]], {i, Length[funcs]} ]; eqnList0 = (eqList /. (a__ == b__) :> a-b) /. funcsRules; numberOfEquations = Length[eqList]; replaceRules = F[n_] :> Function[vars, Evaluate[u0[n]*g[Sequence @@ vars]^alpha[n]] ]; eqnList = If[Head[#]===Plus, List @@ #, # ]& /@ ExpandAll[eqnList0 /. replaceRules]; alphaList0 = Exponent[#, g[Sequence @@ vars] ]& /@ #& /@ eqnList; If[debug, Print["The possible lead terms are:"]; Print[alphaList0]]; alphaList = Map[ If[Head[#]===Plus, List @@ #, If[Head[#]===alpha || Head[#]===Times, {#}, # ] ]&, alphaList0, 2]; If[debug, Print["The lead terms where List has been applied:"]; Print[alphaList]]; alphaListStructure = Table[ Union[# /. {a_Integer, b__} :> {b}& /@ alphaList[[i]] ], {i, numberOfEquations} ]; If[debug, Print["The underlying alphaListStructure is: "]; Print[alphaListStructure]]; alphaList = Table[ Flatten[ Table[ Cases[#, Flatten[ {_, alphaListStructure[[i,j]]} ] | alphaListStructure[[i,j]] ]& /@ {alphaList[[i]]}, {j, Length[alphaListStructure[[i]] ]} ], 1 ], {i, numberOfEquations} ]; If[debug, Print["The alphaList sorted by alphaListStructure: "]; Print[alphaList]]; alphaList = {Min[# /. {a_, ___}:>If[IntegerQ[a], a, 0]& /@ #]}& /@ #& /@ alphaList; alphaList = Plus @@ #& /@ #& /@ Table[Table[Join[alphaList[[i,j]], alphaListStructure[[i,j]]], {j, Length[alphaList[[i]] ]}], {i, numberOfEquations}]; If[debug, Print["The alphaList is: "];Print[alphaList]]; alphaList = Table[If[Length[ alphaList[[i]] ]===1, PadRight[alphaList[[i]], 2], alphaList[[i]] ], {i, numberOfEquations}]; perm := Permutations[Range[Length[alphaList[[j]]]]]; 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[alphaList[[j]] ]}], {j, numberOfEquations}], 2]; If[debug1, Print["comboList1: "];Print[comboList]]; alphaRules = Table[Union[Flatten[Table[ Solve[Equal @@ Flatten[ Table[ alphaList[[comboList[[i, 1]], comboList[[i, k]] ]], {k, 2, Length[comboList[[i]] ]}]], alpha[j]], {i, 1, Length[comboList]}]]], {j, 1, numberOfEquations}]; If[debug1, Print["alphaRules = "];Print[alphaRules]]; comboList2 = Partition[Flatten[ Fold[Table, Array[beta, numberOfEquations], Table[{beta[i], Length[alphaRules[[i]] ]}, {i,numberOfEquations}] ] ], numberOfEquations]; If[debug1, Print["comboList2: "];Print[comboList2]]; alphaSoln = Table[Solve[ Table[Equal @@ alphaRules[[i, comboList2[[j, i]] ]], {i, numberOfEquations}], Array[alpha, numberOfEquations]], {j, Length[comboList2]}]; alphaSoln = Union[Table[Flatten[alphaSoln[[i]]], {i, Length[alphaSoln]}]]; If[debug2, Print["The resulting solutions are:"];Print[alphaSoln]]; alphaSoln2 = Complement[alphaSoln, Flatten[Union[ Table[ If[NumericQ[ Evaluate[alpha[j] /. alphaSoln[[i]] ]], alphaSoln[[i]], Null ], {i, Length[alphaSoln]}, {j, numberOfEquations}]], 1]]; alphaSoln = Complement[alphaSoln, Flatten[Union[ Table[ If[IntegerQ[ Evaluate[alpha[j] /. alphaSoln[[i]] ]], Null, alphaSoln[[i]] ], {i, Length[alphaSoln]}, {j, numberOfEquations}]], 1]]; alphaSoln = Complement[alphaSoln, Flatten[Union[ Table[ If[Negative[ Evaluate[alpha[j] /. alphaSoln[[i]] ]], Null,alphaSoln[[i]] ], {i, Length[alphaSoln]}, {j, numberOfEquations}]], 1]]; alphaSoln0 = alphaSoln; If[debug2, Print["Of these, we are left with:"]; Print[alphaSoln]]; checkList2= Table[DeleteCases[Union[Flatten[Table[ Table[If[ Evaluate[alphaList0[[i,k]] /. alphaSoln[[j]] ]== Min[alphaList0[[i]]/.alphaSoln[[j]]], alphaList0[[i,k]]], {k, Length[alphaList0[[i]]]}], {j, Length[alphaSoln]}]]], Null], {i,numberOfEquations}]; If[debug4, Print["The alpha terms that are of highest order are: "]; Print[checkList2]]; freeVariable = DeleteCases[Flatten[Table[ If[Length[Cases[checkList2[[i]], ___Integer+__*alpha[j]|___Integer+alpha[j] ] ]>1 && numberOfEquations>1, alpha[j], Null], {j, numberOfEquations}, {i, numberOfEquations}]], Null]; If[Length[freeVariable]!=0, freeVariableMin = Min[Table[freeVariable[[j]] /. alphaSoln[[i]], {i, Length[alphaSoln]}, {j, Length[freeVariable]}] ]]; If[debug4, Print["The free variable(s) is/are: "];Print[freeVariable]; Print["The minimum value of either being: "];Print[freeVariableMin]]; If[Length[freeVariable]!=0, alphaSolnFree = Partition[Flatten[ Fold[Table, Table[alpha[i]->beta[i], {i, numberOfEquations}], Table[{beta[k],freeVariableMin,-1},{k,numberOfEquations}]]], numberOfEquations]]; If[Length[freeVariable]!=0, alphaSolnFree = Complement[DeleteCases[Table[If[Plus @@ Flatten[Table[ If[Length[Cases[alphaList0[[i]] /. alphaSolnFree[[j]], Min[alphaList0[[i]] /. alphaSolnFree[[j]]] ]]>=2, 1, 0], {i, numberOfEquations}]]>=numberOfEquations, alphaSolnFree[[j]], Null], {j, Length[alphaSolnFree]}], Null], alphaSoln]]; Do[If[Length[Cases[checkList2[[i]], ___Integer+__*alpha[j]|___Integer+alpha[j] ] ]>1 && numberOfEquations>1, Message[alphaSolver::freedom, Cases[checkList2[[i]], ___Integer+__*alpha[j]|___Integer+alpha[j] ], i] ], {j, numberOfEquations}, {i, numberOfEquations}]; alphaSoln = DeleteCases[ Table[ If[Plus @@ Flatten[ Table[If[Length[Cases[alphaList0[[i]] /. alphaSoln[[j]], Min[alphaList0[[i]] /. alphaSoln[[j]] ] ]]>=2, 1, 0], {i, numberOfEquations}]]>=numberOfEquations, alphaSoln[[j]], Null ], {j, Length[alphaSoln]}], Null]; If[debug1, Print["alphaSoln at point ZZZZ is: "];Print[alphaSoln]]; If[Length[alphaSolnFree]!=0, alphaSoln = Join[alphaSoln, alphaSolnFree]]; On[Solve::svars]; On[Solve::ifun]; If[debug4, Print["AlphaSoln: "];Print[alphaSoln]]; If[ Length[alphaSoln] === 0, Message[alphaSolver::noSoln]; alphaSoln = Partition[Flatten[ Fold[Table, Table[alpha[i]->beta[i], {i, numberOfEquations}], Table[{beta[i], -3, -1}, {i,numberOfEquations}] ] ], numberOfEquations]; alphaSoln = DeleteCases[ Table[ If[Plus @@ Flatten[ Table[If[Length[Cases[alphaList0[[i]] /. alphaSoln[[j]], Min[alphaList0[[i]] /. alphaSoln[[j]] ] ]]>=2, 1, 0], {i, numberOfEquations}]]>=numberOfEquations, alphaSoln[[j]], Null ], {j, Length[alphaSoln]}], Null], alphaSoln = alphaSoln]; If[Length[alphaSoln] == 0, Message[alphaSolver::fail]; Abort[ ], alphaSoln] ]; (* Initial Equation Solver: *) debug2 = debug3 = False; initialSolver[eqList_List, funcs_List, vars_List]:= Block[{alphaSoln, u0Soln}, alphaSoln = alphaSolver[eqList, funcs, vars]; u0Soln = initialSolve[eqList, funcs, vars, #]& /@ alphaSoln; u0Soln = DeleteCases[u0Soln /. {{__},{___,u[_,0][__]->0,___}}:>Null, Null, Infinity]; If[Length[u0Soln] == 0, Message[initialSolver::fail];Abort[ ], Return[u0Soln]] ]; initialSolve[eqList_List, funcs_List, vars_List, alphaSoln_List, prnt_:True]:= Block[{newFuncs = Array[F, Length[funcs]], funcsRules, eqnList0, numberOfEquations, replaceRules, eqnList, uSoln0eqnList, uSoln0}, Off[Solve::svars]; funcsRules = Table[Head[funcs[[i]]]->newFuncs[[i]], {i, Length[funcs]}]; eqnList0 = (eqList /. (a__ == b__):>a-b) /. funcsRules; numberOfEquations = Length[eqList]; If[prnt, Print["Leading order analysis determines: "]; Print[alphaSoln /. (alpha[a_]->b_):>(Subscript[\[Alpha], a]->b)]]; replaceRules = F[n_] :> Function[vars, Evaluate[u[n,0][Sequence @@ vars]*g[Sequence @@ vars]^alpha[n]]]; eqnList = (eqnList0 /. replaceRules) /. alphaSoln; eqnList = Expand[Divide[#, (PolynomialGCD @@ If[Head[#] === Plus, List @@ #, {#}])]] & /@ eqnList; eqnList0 = eqnList; uSoln0eqnList = Coefficient[#, g[Sequence @@ vars], 0]& /@ eqnList0; uSoln0eqnList = Expand[Divide[#, (PolynomialGCD @@ If[Head[#] === Plus, List @@ #, {#}])]] & /@ uSoln0eqnList; uSoln0eqnList = (# == 0)& /@ uSoln0eqnList; uSoln0 = Solve[uSoln0eqnList, Sort[Table[u[i,0][Sequence @@ vars], {i, numberOfEquations}], Greater]]; On[Solve::svars]; Return[{alphaSoln, #}& /@ uSoln0]]; (* Resonance Finder *) debug5 = debug6 = False; resonanceSolver[eqList_List, funcs_List, vars_List, prnt_:True]:= Block[{alphaSoln, resSoln}, alphaSoln = alphaSolver[eqList, funcs, vars]; resSoln = {#1, #2, resonanceSolve[eqList, funcs, vars, #1, #2, prnt] }& @@ #& /@ initialSolve[eqList, funcs, vars, #, prnt]& /@ alphaSoln; resSoln = Flatten[resSoln, 1]; resSoln = DeleteCases[resSoln //. {{__}, {__}, Null}->Null, Null]; resSoln = DeleteCases[resSoln //. {___, List[], ___}->Null, Null]; If[Length[resSoln] == 0, Message[resonanceSolver::fail];Abort[ ], resSoln] ]; resonanceSolve[eqList_List, funcs_List, vars_List, alphaSoln_, u0Soln_, prnt_:True]:= Block[{newFuncs = Array[F, Length[funcs]], funcsRules, eqnList, numberOfEquations, replaceRules, (* eqnReplaceRules,termList, *) rList0, rList, rListMin, matrixQ, detQ, resSoln}, funcsRules = Table[Head[funcs[[i]]]->newFuncs[[i]], {i, Length[funcs]}]; eqnList0 = (eqList /. (a__ == b__):>a-b) /. funcsRules; numberOfEquations = Length[eqList]; replaceRules = Flatten[Join[alphaSoln, toPure[{u0Soln}], {F[n_] :> Function[vars, Evaluate[u[n,0][Sequence @@ vars]*g[Sequence @@ vars]^alpha[n]+ u[n,r][Sequence @@ vars]*g[Sequence @@ vars]^(alpha[n]+r)]]}]]; If[debug5, Print["The replacement rules for this function are:"]; Print[replaceRules]]; eqnList = eqnList0 //. replaceRules; eqnList = Expand[Divide[#, (PolynomialGCD @@ If[Head[#] === Plus, List @@ #, {#}])]]& /@ eqnList; matrixQ = Table[Table[ Coefficient[ Coefficient[eqnList[[i]], u[j,r][Sequence @@ vars]*g[Sequence @@ vars]^r], g[Sequence @@ vars], 0], {i, numberOfEquations}], {j, numberOfEquations}]; If[prnt, StylePrint["The matrix for the resonances is:", "Text"]; Print[MatrixForm[pdeForm /@ matrixQ]]]; detQ = Factor[ExpandAll[Det[matrixQ]]]; If[prnt, StylePrint["The characteristic equation for the resonances is:", "Text"]; Print[pdeForm[detQ]==0]]; resSoln = Solve[detQ==0, r]; If[FreeQ[resSoln, r->-1], Return[Null]]; resSoln = DeleteCases[ If[FreeQ[resSoln, (u[_,0][Sequence @@ vars]|Derivative[__][g][__])], resSoln, Null], Null]; resSoln = Flatten[{resSoln}]; If[prnt, StylePrint["The resonances are: ", "Text"]; Print[resSoln]; If[Length[DeleteCases[Select[resSoln /. (r->a_):>a, Negative], -1]]>0, CellPrint[Cell["This case has negative resonances, other than r\[Rule]-1, which will be ignored.", "Message"]]]]; (* Removed on 8-8-00 to show which cases were ignored. *) (* If[Length[resSoln]>=1, resSoln, Return[Null]]; resSoln = If[Length[resSoln] === Length[Cases[#[[2]]& /@ (List @@ #& /@ resSoln), _Integer]], resSoln, If[prnt, Message[resonanceSolver::fail2, DeleteCases[resSoln, r->_Integer]]];resSoln];*) Return[resSoln] ]; (* : Title : Constants of Integration Solver *) (* : Summary : Completes the last step of the ARS algorithm. *) debug7 = debug8 = debug9 = debug10 = debug11 = debug12 = False; PainleveTest[eqList_List, funcs_List, vars_List, opt___] := Block[{eqnList, newFuncs = Array[F, Length[funcs]], funcsRules, resSoln, uSoln, time}, If[!FreeQ[eqList, Power[__, _Symbol]], Message[PainleveTest::symbol];Abort[]]; If[!FreeQ[eqList, Power[E, __]], Message[PainleveTest::exp];Abort[]]; funcsRules = Table[Head[funcs[[i]]]->newFuncs[[i]], {i, Length[funcs]}]; eqnList = (eqList /. (a__==b__):>(a-b)) /. funcsRules; time = TimeUsed[ ]; StylePrint["Painleve analysis of the equation(s):", "Text"]; Print[TableForm[pdeForm[(# == 0)&/@eqnList]]]; CellPrint[Cell[ TextData[{"Using the Laurent series of the solution ", Cell[BoxData[FormBox[RowBox[{SubscriptBox["u", "i"], ToBoxes[vars] /. {"{"->"(", "}"->")"}}], TraditionalForm]]], ":"}], "Text"]]; Print[SequenceForm[Subscript[u,i], StringReplace[ToString[vars], {"{" -> "(", "}" -> ")", ", " -> "," }], " = ",Power[g,Subscript[\[Alpha],i]],StringReplace[ToString[vars], {"{"->"(","}"->")",", "->","}], HoldForm[Sum[" ", {k,0,\[Infinity]}]], HoldForm[Subscript[u,"(i,k)"]],StringReplace[ToString[vars], {"{"->"(","}"->")",", "->","}],Power[g,k], StringReplace[ToString[vars], {"{"->"(","}"->")",", "->","}]]]; resSoln = resonanceSolver[eqnList, funcs, vars, False]; StylePrint["Leading order analysis determines "<>ToString[Length[resSoln]]<> " possible case(s). ", "Text"]; Print[TableForm[ {#[[1]] /. (alpha[a_]->b_):>(Subscript[\[Alpha],a]->b), pdeForm[#[[2]]], #[[3]]}& /@ resSoln]]; resSoln = Select[resSoln, (Length[#[[3]]]===Length[Select[#[[3]] /. (r->a_):>a, IntegerQ]] && NonNegative[Max[Select[#[[3]] /. (r->a_):>a, IntegerQ]]])&]; If[Length[resSoln]===0, Message[resonanceSolver::fail];Abort[]]; StylePrint["The "<>ToString[Length[resSoln]]<>" case(s) whose resonances are all integers and whose resonances are all non-negative. These cases will now be evaluated one at a time.", "Text"]; Print[TableForm[ {#[[1]] /. (alpha[a_]->b_):>(Subscript[\[Alpha],a]->b), pdeForm[#[[2]]], #[[3]]}& /@ resSoln]]; uSoln = {#[[1]], #[[3]], constantsSolve[eqnList, funcs, vars, #[[1]], #[[2]], #[[3]], opt]}& /@ resSoln; time = Round[(TimeUsed[ ] - time)*100]/100.0; If[!FreeQ[{opt}||{}, WorkLog->_], Save[WorkLog /. {opt}, time]]; If[Length[uSoln] == 0, Message[constantsSolver::fail]; Abort[]]; StylePrint["Painleve analysis of this system determined: ", "Text"]; (StylePrint["The dominate behaviour is: ","Text"]; Print[#[[1]] /. alpha[i_]:>Subscript[\[Alpha], i]]; StylePrint["The resonances are: ", "Text"];Print[#[[2]]]; StylePrint["The constants of integration are: ", "Text"]; Print[pdeForm[#[[3]]]])& /@ uSoln; Return[uSoln]]; PainleveTest[eqList_List, funcs_, vars_, opt___] := PainleveTest[eqList, If[Head[funcs]===List, funcs, {funcs}], If[Head[vars]===List, vars, {vars}], opt]; constantsSolve[eqList_List, funcs_List, vars_List, alphaSoln_List, u0Soln_List, resSoln_List, opt___] := Block[{numberOfEquations, maxRes, resList, fail, eqnList, replaceRules, minGPower, uSoln, uSolnMatrix1, uSolnMatrix2, uSolnEqnList, uSoln0eqnList, temp2, eqnReplaceRules, termList, minPhiPower, phiPower, eqnSolnList, eqnSolnList0, uSolni, solution, ultUSoln, eqnList0, tempUSoln, kruskal, kruskalG, conditions}, Off[\[Infinity]::indet, Power::infy]; Unprotect[$MessageList]; Do[$MessageList = Append[$MessageList, HoldForm[Solve::svars]],{3}]; numberOfEquations = Length[eqList]; StylePrint["_______________________________________________________________", "Text"]; StylePrint["In this case, leading order analysis determined: ", "Text"]; Print[alphaSoln /. (alpha[a_]->b_):>(Subscript[\[Alpha],a]->b)]; resList = (r /. #)& /@ resSoln; resList = Cases[resList, _Integer]; maxRes = Max[r /. #& /@ resSoln]; If[!NonNegative[maxRes], StylePrint["This case does not have a positive integer resonance.", "Text"]; Return[Null]]; If[debug7, Print["For which the maximum integer value is: ", maxRes]]; fail = {}; If[FreeQ[{opt}||{}, Kruskal], kruskalG = {}, kruskalG = {g :> Function[vars, Evaluate[(Kruskal /. {opt})- h[Sequence @@ (Complement[vars, {Kruskal /. {opt}}])]]], Derivative[var__][u[_,_]][Sequence @@ vars]:>0 /; {var}[[Position[vars, (Kruskal /. {opt})][[1,1]]]] != 0}]; uSoln = toPure[u0Soln /. kruskalG]; CellPrint[Cell[ TextData[{"The values(s) of ", Cell[BoxData[FormBox[SubscriptBox["u", "(i,0)"], TraditionalForm]]], " are: "}], "Text"]]; Print[pdeForm[Table[u[i,0] -> (u[i,0] /. uSoln) /. Function[vars, a_]:>a, {i,numberOfEquations}]]]; resonanceSolve[eqList, funcs, vars, alphaSoln, u0Soln]; replaceRules = F[n_] :> Function[vars, Sum[Evaluate[u[n,m][Sequence @@ vars]*g[Sequence @@ vars]^(m+alpha[n])], {m, 0, maxRes}]]; If[debug7, Print["The replacement rules for this function are:"]; Print[replaceRules]]; eqnList = eqList /. replaceRules; If[debug8, Print["The equation list, eqnList, after replacement is:"]; Print[eqnList]]; eqnList = eqnList /. alphaSoln; If[debug8, Print["The equation list after the alpha[i] are fixed is:"]; Print[eqnList]]; eqnList = Expand[Divide[#, (PolynomialGCD @@ If[Head[#] === Plus, List @@ #, {#}])]] & /@ eqnList; If[debug8, Print["The equation list after expanding and factoring is:"]; Print[eqnList]]; eqnList0 = eqnList; conditions = {}; Do[ uSolnEqnList = Coefficient[#, g[Sequence @@ vars], ii]& /@ eqnList0; uSolnEqnList = uSolnEqnList //. uSoln; If[debug12, Print["uSolnEqnList is: "]; (Print[pdeForm[#]];Print[pdeForm[(# /. kruskalG)]])&[uSolnEqnList]]; uSolnMatrix1 = Table[Table[ Expand[Coefficient[uSolnEqnList[[i]], u[j, ii][Sequence @@ vars]]], {j, numberOfEquations}], {i, numberOfEquations}]; uSolnMatrix2 = Table[-(Plus @@ DeleteCases[ If[Head[#]===Plus, List @@ #, {#}]&[Expand[uSolnEqnList[[i]] ]], (__*u[_,ii][Sequence @@ vars]|u[_,ii][Sequence @@ vars])]), {i, numberOfEquations}]; uSolnMatrix1 = uSolnMatrix1 /. kruskalG; uSolnMatrix2 = uSolnMatrix2 /. kruskalG; If[MemberQ[resList, ii], StylePrint["The system at resonance level k="<>ToString[ii]<>" is:", "Text"], StylePrint["The system at non-resonance level k="<>ToString[ii]<>" is:", "Text"]]; Print[MatrixForm[pdeForm[uSolnMatrix1]], MatrixForm[pdeForm[Table[{u[i,ii]}, {i, numberOfEquations}]]], " = "];Print[If[FreeQ[(Length[#]>7)&/@#, True], MatrixForm[#], #]& [pdeForm[uSolnMatrix2]]]; tempUSoln = If[FreeQ[uSolnMatrix1, Table[0, {numberOfEquations}]], Solve[uSolnMatrix1.Table[u[i,ii][Sequence @@ vars], {i, numberOfEquations}] == uSolnMatrix2, Table[u[numberOfEquations-i+1,ii][Sequence @@ vars], {i, numberOfEquations}]], If[FreeQ[Table[Append[#1[[i]], -#2[[i]]], {i, numberOfEquations}]&[uSolnMatrix1,uSolnMatrix2], Table[0, {numberOfEquations+1}]], If[numberOfEquations==1, fail = Append[fail, ii]; (uSoln = Flatten[Join[uSoln, Flatten[CellPrint[Cell["Compatibility conditions:", "Message"]]; StylePrint[pdeForm[#[[2]]], "Message"]; conditions = Flatten[Join[conditions, {#[[2]]}]];#[[1]]]]]; #[[1]])&[mySolve[ uSolnMatrix1.Table[u[i,ii][Sequence @@ vars], {i, numberOfEquations}] == uSolnMatrix2, ii]], {}], If[numberOfEquations==1, Print["The compatibility condition is satisfied."];{}, Solve[uSolnMatrix1.Table[u[i,ii][Sequence @@ vars], {i, numberOfEquations}] == uSolnMatrix2, Table[u[numberOfEquations-i+1,ii][Sequence @@ vars], {i, numberOfEquations}]]]]]; tempUSoln = Map[Factor[Expand[#]]&, tempUSoln, Infinity]; If[Length[Flatten[tempUSoln]]>0, If[ii===1, StylePrint["The solutions at level "<>ToString[ii]<>" is: ", "Text"], StylePrint["The solutions at level "<>ToString[ii]<>" are: ", "Text"]]; Print[pdeForm[tempUSoln]]]; If[debug11, Print["tempUSoln using Solve is:", toPure[tempUSoln]]]; If[FreeQ[resList, ii], If[numberOfEquations === Length[Flatten[tempUSoln]], uSoln = Flatten[Join[uSoln, toPure[tempUSoln]]]; If[debug11, Print["uSoln is: ",uSoln]], fail = Append[fail, ii]; uSoln = Flatten[Join[uSoln, Flatten[(CellPrint[Cell["Compatibility conditions:", "Message"]]; StylePrint[pdeForm[#[[2]]], "Message"]; conditions = Flatten[Join[conditions, {#[[2]]}]];#[[1]])&[mySolve[ uSolnMatrix1.Table[u[i,ii][Sequence @@ vars], {i, numberOfEquations}] == uSolnMatrix2, ii]]]]]], If[Det[uSolnMatrix1] === 0, If[numberOfEquations - Length[Cases[resList, ii]] === Length[Flatten[tempUSoln]], uSoln = Flatten[Join[uSoln, toPure[tempUSoln]]]; temp2 = Length[Cases[resList, ii]]; If[temp2>1, StylePrint["There are "<>ToString[temp2]<>" free functions.", "Text"], If[numberOfEquations != 1, StylePrint["There is one free function.", "Text"]]], fail = Append[fail, ii]; uSoln = Flatten[Join[uSoln, Flatten[(CellPrint[Cell["Compatibility conditions:", "Message"]]; StylePrint[pdeForm[#[[2]]], "Message"]; conditions = Flatten[Join[conditions, {#[[2]]}]];#[[1]])&[mySolve[ uSolnMatrix1.Table[u[i,ii][Sequence @@ vars], {i, numberOfEquations}] == uSolnMatrix2, ii]]]]]], fail = Append[fail, ii]; uSoln = Flatten[Join[uSoln, Flatten[(CellPrint[Cell["Compatibility conditions:", "Message"]]; StylePrint[pdeForm[#[[2]]], "Message"]; conditions = Flatten[Join[conditions, {#[[2]]}]];#[[1]])&[mySolve[ uSolnMatrix1.Table[u[i,ii][Sequence @@ vars], {i, numberOfEquations}] == uSolnMatrix2, ii]]]]]]], {ii, maxRes}]; ultUSoln = Join[ Table[u[j,i] -> Evaluate[u[j,i][Sequence @@ vars] /. uSoln], {i,0, maxRes}, {j, numberOfEquations}], {conditions}]; ultUSoln = DeleteCases[ultUSoln, List[]|{}, Infinity]; ii = 0; ultUSoln = If[Length[ultUSoln]==1, If[Length[ultUSoln[[1]] ]==1, Flatten[ultUSoln], Flatten[ultUSoln, 1]], ultUSoln] /. (u[a_,b_]->u[_,_][Sequence @@ vars]):>If[MemberQ[fail, b],(u[a,b]->$Undetermined), (ii++;u[a,b]->C[ii][Sequence @@ vars])]; Protect[$MessageList]; On[\[Infinity]::indet, Power::infy]; equations = eqList; alphas = alphaSoln; resonances = resSoln; constants = ultUSoln; If[!FreeQ[{opt}||{}, WorkLog->_], Save[WorkLog /. {opt}, equations, alphas, resonances, constants, ultUSoln]]; Return[ultUSoln]]; (* : Title : mySolve *) (* : Summary : This ``Solver'' attempts to determine the conditions when the Solve function fails to determine the constants of integration. *) mySolveDebug = False; mySolve[system_, level_]:= Block[{newSystem, numberOfEquations,allVar, param, iVar, iVarSoln, iVarSystem, iVarSystem2}, newSystem = DeleteCases[system, True]; newSystem = (#[[1]]-#[[2]])& /@ Transpose[(system /. Equal->List)]; newSystem = Factor[Expand[newSystem]]; If[mySolveDebug, Print["The newSystem is: "];Print[newSystem]]; iVarSoln = {}; iVarSystem = newSystem; iVarSystem = FixedPoint[( allVar = Union[DeleteCases[ DeleteCases[Flatten[# //. {Equal -> List, Times->List, Plus->List, Power->List, Rational->List}], _Integer],_Rational]]; iVar = Cases[allVar, u[_,level][__]]; If[Length[iVar]>0 && Length[#]>1, iVarSoln = Join[iVarSoln, toPure[{Flatten[{(Solve[# == 0, iVar[[1]] ]& /@ DeleteCases[If[FreeQ[#,iVar[[1]]], Null, #]& /@ #, Null])[[1]]}]}]]; iVarSystem = DeleteCases[(# == 0)& /@ Factor[Expand[Flatten[# /. iVarSoln]]], True]; (iVarSystem /. Equal[a__,0]:>a), iVarSystem])&, newSystem]; Return[{iVarSoln, DeleteCases[(# == 0)& /@ (iVarSystem /. Equal[a__,b__]:>a-b), True, Infinity]}] ]; PainleveTest[]:= Block[{choice}, StylePrint["Welcome to the Painleve Test Package for Mathematica 4.0:\[NewLine]"<> "By Douglas Baldwin and Willy Hereman\[NewLine]"<> "\[Copyright] 2000, Colorado School of Mines", "Subtitle"]; StylePrint["The built-in test cases are:", "Section"]; StylePrint["Single Equations:", "Subsection"]; StylePrint[" 1) Fisher Equation;\[NewLine]"<> " 2) Cylindrical KdV Equation;\[NewLine]"<> " 3) Burgers' Equation;\[NewLine]"<> " 4) Korteweg de Vries Equation;\[NewLine]"<> " 5) Modified Korteweg de Vries Equation;\[NewLine]"<> " 6) Ito Equation;\[NewLine]"<> " 7) Boussinesq Equation;\[NewLine]"<> " 8) Non-Linear Sine-Gorden Equation\[NewLine]"<> " 9) Liouville Non-Linear Klein-Gorden Equation", "Text"]; StylePrint["Systems of Equations:", "Subsection"]; StylePrint[" 10) Hirota Satsuma w/ parameter;\[NewLine]"<> " 11) Hirota Satsuma w/o parameter;\[NewLine]"<> " 12) One-Dimensional Non-Linear Schrodinger Equations;\[NewLine]"<> " 13) Reduced Maxwell Bloch Equations;\[NewLine]"<> " 14) Lorenz Model of a Dissipative System;\[NewLine]"<> " 15) Hu System;\[NewLine]\[NewLine]"<> " x) Exit", "Text"]; choice = Input["Enter your choice: "]; Switch[choice, x, Abort[], 1, StylePrint["The Fisher Equation:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t], {x, 2}]+c*D[u[x,t], x]-u[x,t]^2+u[x,t]==0}, {u[x,t]}, {x,t}, Kruskal->x]]]; PainleveTest[{D[u[x,t], {x, 2}]+c*D[u[x,t], x]-u[x,t]^2+u[x,t]==0}, {u[x,t]}, {x,t}, Kruskal->x], 2, StylePrint["The Cylindrical KdV Equation:", "Section"]; Print[HoldForm[PainleveTest[{a[t]*D[u[x,t], x]+D[u[x,t], {x, 4}]+ 6*u[x,t]*D[u[x,t], {x, 2}]+6*D[u[x,t], x]^2+ D[u[x,t], {x, 1}, {t, 1}]}, {u[x,t]}, {x,t}]]]; PainleveTest[{a[t]*D[u[x,t], x]+D[u[x,t], {x, 4}]+ 6*u[x,t]*D[u[x,t], {x, 2}]+6*D[u[x,t], x]^2+ D[u[x,t], {x, 1}, {t, 1}]}, {u[x,t]}, {x,t}], 3, StylePrint["The Burgers' Equation:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t], t]+u[x,t]*D[u[x,t], x]- sigma*D[u[x,t], {x,2}]}, {u[x,t]}, {x,t}]]]; PainleveTest[{D[u[x,t], t]+u[x,t]*D[u[x,t], x]- sigma*D[u[x,t], {x,2}]}, {u[x,t]}, {x,t}], 4, StylePrint["The Korteweg de Vries Equation:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t], t]+ 6*u[x,t]*D[u[x,t], x]+D[u[x,t], {x,3}]}, {u[x,t]}, {x,t}]]]; PainleveTest[{D[u[x,t], t]+ 6*u[x,t]*D[u[x,t], x]+D[u[x,t], {x,3}]}, {u[x,t]}, {x,t}], 5, StylePrint["The Modified Korteweg de Vries Equation:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t], t]-3* u[x,t]^2*D[u[x,t], x]+ 2*sigma^2*D[u[x,t], {x,3}]},{u[x,t]}, {x,t}]]]; PainleveTest[{D[u[x,t], t]-3* u[x,t]^2*D[u[x,t], x]+2*sigma^2* D[u[x,t], {x,3}]},{u[x,t]}, {x,t}], 6, StylePrint["The Ito Equation:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t], t]+a*u[x,t]^2*D[u[x,t], x]+ b*D[u[x,t], x]*D[u[x,t], {x,2}] + c*u[x,t]*D[u[x,t], {x, 3}]+ D[u[x,t], {x,5}]} /. {a->2, b->6, c->3}, {u[x,t]}, {x,t}, Kruskal->t]]]; PainleveTest[{D[u[x,t], t]+a*u[x,t]^2*D[u[x,t], x]+ b*D[u[x,t], x]*D[u[x,t], {x,2}] + c*u[x,t]*D[u[x,t], {x, 3}]+ D[u[x,t], {x,5}]} /. {a->2, b->6, c->3}, {u[x,t]}, {x,t}, Kruskal->t], 7, StylePrint["The Boussinesq Equation:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t], {t,2}]+2*u[x,t]*D[u[x,t], {x,2}]+ 2*D[u[x,t], x]^2+(1/3)*D[u[x,t], {x,4}]}, {u[x,t]}, {x,t}]]]; PainleveTest[{D[u[x,t], {t,2}]+2*u[x,t]*D[u[x,t], {x,2}]+ 2*D[u[x,t], x]^2+(1/3)*D[u[x,t], {x,4}]}, {u[x,t]}, {x,t}], 8, StylePrint["The Non-Linear Sine-Gorden Equation:", "Section"]; Print[HoldForm[PainleveTest[{u[x,t]*(D[u[x,t], {t, 2}]-D[u[x,t], {x,2}])- D[u[x,t], t]^2+D[u[x,t], x]^2 - (1/2)*u[x,t]^3+(1/2)*u[x,t]}, {u[x,t]}, {x,t}]]]; PainleveTest[{u[x,t]*(D[u[x,t], {t, 2}]-D[u[x,t], {x,2}])- D[u[x,t], t]^2+D[u[x,t], x]^2 - (1/2)*u[x,t]^3+(1/2)*u[x,t]}, {u[x,t]}, {x,t}], 9, StylePrint["Liouville Non-Linear Klein-Gorden Equation:", "Section"]; Print[HoldForm[PainleveTest[{u[x,t]*(D[u[x,t], {t, 2}]-D[u[x,t], {x,2}])- D[u[x,t], t]^2+D[u[x,t], x]^2 - u[x,t]^3}, {u[x,t]}, {x,t}]]]; PainleveTest[{u[x,t]*(D[u[x,t], {t, 2}]-D[u[x,t], {x,2}])- D[u[x,t], t]^2+D[u[x,t], x]^2 - u[x,t]^3}, {u[x,t]}, {x,t}], 10, StylePrint["The Hirota Satsuma w/ parameter:", "Section"]; Print[HoldForm[PainleveTest[{-aa*D[u[x,t], {x, 3}]- 6*aa*u[x,t]*D[u[x,t], x]+6*v[x,t]*D[v[x,t], x]+D[u[x,t], t], D[v[x,t], {x, 3}]+3*u[x,t]*D[v[x,t], x]+D[v[x,t], t]}, {u[x,t],v[x,t]}, {x,t}, Kruskal->x]]]; PainleveTest[{-aa*D[u[x,t], {x, 3}]-6*aa*u[x,t]*D[u[x,t], x]+6* v[x,t]*D[v[x,t], x]+D[u[x,t], t], D[v[x,t], {x, 3}]+ 3*u[x,t]*D[v[x,t], x]+D[v[x,t], t]}, {u[x,t],v[x,t]}, {x,t}, Kruskal->x], 11, StylePrint["The Hirota Satsuma w/o parameter:", "Section"]; Print[HoldForm[PainleveTest[{-aa*D[u[x,t], {x, 3}]- 6*aa*u[x,t]*D[u[x,t], x]+6*v[x,t]*D[v[x,t], x]+D[u[x,t], t], D[v[x,t], {x, 3}]+3*u[x,t]*D[v[x,t], x]+D[v[x,t], t]} /. aa->1/2, {u[x,t],v[x,t]}, {x,t}, Kruskal->x]]]; PainleveTest[{-aa*D[u[x,t], {x, 3}]-6*aa*u[x,t]*D[u[x,t], x]+ 6*v[x,t]*D[v[x,t], x]+D[u[x,t], t], D[v[x,t], {x, 3}]+ 3*u[x,t]*D[v[x,t], x]+D[v[x,t], t]} /. aa->1/2, {u[x,t],v[x,t]}, {x,t}, Kruskal->x], 12, StylePrint["The One-Dimensional Non-Linear Schrodinger Equations:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t], t]+D[v[x,t], {x,2}]+ a*v[x,t]*(u[x,t]^2+v[x,t]^2), D[v[x,t], t]-D[u[x,t], {x,2}]- a*u[x,t]*(u[x,t]^2+v[x,t]^2)}, {u[x,t], v[x,t]}, {x,t}, Kruskal->x]]]; PainleveTest[{D[u[x,t], t]+D[v[x,t], {x,2}]+ a*v[x,t]*(u[x,t]^2+v[x,t]^2), D[v[x,t], t]-D[u[x,t], {x,2}]- a*u[x,t]*(u[x,t]^2+v[x,t]^2)}, {u[x,t], v[x,t]}, {x,t}, Kruskal->x], 13, StylePrint["The Reduced Maxwell Bloch Equations:", "Section"]; Print[HoldForm[PainleveTest[{D[u[x,t],x]+mu*v[x,t], D[v[x,t], x]- e[x,t]*w[x,t]-mu*u[x,t], D[w[x,t], x]+e[x,t]*v[x,t], D[e[x,t], t]+ v[x,t]}, {u[x,t],v[x,t],w[x,t],e[x,t]}, {x,t}, Kruskal->x]]]; PainleveTest[{D[u[x,t],x]+mu*v[x,t], D[v[x,t], x]-e[x,t]*w[x,t]- mu*u[x,t], D[w[x,t], x]+e[x,t]*v[x,t], D[e[x,t], t]+v[x,t]}, {u[x,t],v[x,t],w[x,t],e[x,t]}, {x,t}, Kruskal->x], 14, StylePrint["The Lorenz Model of a Dissipative System:", "Section"]; Print[HoldForm[PainleveTest[{v[x,t]-D[u[x,t], t], u[x,t]-u[x,t]*w[x,t]- D[v[x,t], t], u[x,t]*v[x,t]-D[w[x,t], t]}, {u[x,t], v[x,t], w[x,t]}, {x,t}, Kruskal->x]]]; PainleveTest[{v[x,t]-D[u[x,t], t], u[x,t]-u[x,t]*w[x,t]-D[v[x,t], t], u[x,t]*v[x,t]-D[w[x,t], t]}, {u[x,t], v[x,t], w[x,t]}, {x,t}, Kruskal->x], 15, StylePrint["The Hu System:", "SubSection"]; Print[HoldForm[PainleveTest[{D[u[x,t], t]- D[u[x,t], {x, 2}]- 2*u[x,t]*D[u[x,t], x]-D[v[x,t], x], D[v[x,t], t]- D[v[x,t], {x,3}]- 6*D[u[x,t], x]*D[v[x,t], x]}, {u[x,t], v[x,t]}, {x,t}, Kruskal->x]]]; PainleveTest[{D[u[x,t], t]- D[u[x,t], {x, 2}]- 2*u[x,t]*D[u[x,t], x]- D[v[x,t], x], D[v[x,t], t]- D[v[x,t], {x,3}]- 6*D[u[x,t], x]*D[v[x,t], x]}, {u[x,t], v[x,t]}, {x,t}, Kruskal->x], ___, Abort[]]]; Protect[alphaSolver, alphaSolve, resonanceSolver, resonanceSolve, constantsSolve, PainleveTest]; EndPackage[]