(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2"; MacintoshStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L2, e8, 24, "Times"; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L2, e6, 18, "Times"; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, L2, e6, 14, "Times"; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, L2, a20, 18, "Times"; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, L2, a15, 14, "Times"; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, L2, a12, 12, "Times"; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 10, "Times"; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier"; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, L2, 12, "Courier"; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, L2, 10, "Geneva"; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = leftheader, inactive, L2, 12, "Times"; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, L2, 12, "Times"; fontset = leftfooter, inactive, L2, 12, "Times"; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 10, "Times"; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2, 12, "Times"; paletteColors = 128; currentKernel; ] :[font = info; inactive; preserveAspect; startGroup] ********************************************************************** Jet --- A collection of basic jet space routines. Written in Mathematica ********************************************************************** ;[s] 7:0,0;72,1;75,2;78,1;186,3;197,1;199,0;270,-1; 4:2,12,10,Courier,0,12,0,0,65535;3,17,13,Courier,0,17,65535,0,0;1,28,21,Courier,0,28,65535,0,0;1,17,13,Courier,2,17,0,65535,0; :[font = message; inactive; wordwrap; dontPreserveAspect] Peter J. Olver University of Minnesota December, 1996 ;[s] 3:0,0;15,1;110,0;111,-1; 2:2,17,13,Courier,0,18,65535,0,0;1,17,13,Courier,0,18,21844,21844,21844; :[font = input; inactive; wordwrap; dontPreserveAspect; endGroup] These routines implement basic computations in jet spaces, for both commutative and noncommutative variables. The package handles differential functions, differential equations, differential invariants, total and partial derivatives, symmetry groups, variational derivatives, prolonged vector fields, and more ... ;[s] 3:0,1;1,2;342,3;344,-1; 4:0,12,10,Courier,1,12,0,0,0;1,17,13,Courier,1,18,0,0,0;1,14,11,Courier,1,14,0,0,0;1,17,13,Courier,0,18,0,0,0; :[font = info; inactive; preserveAspect] ********************************************************************** Basic Functions and Abbreviations ********************************************************************** ;[s] 3:0,0;72,1;110,0;181,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; initialization; wordwrap; dontPreserveAspect] *) ex = ExpandAll; cf = Coefficient; ff = FullForm; ca = ClearAll; mf = MatrixForm; tg = Together; fa = Factor; pe = PowerExpand; num = Numerator; den = Denominator; tfp[f_] := f //Together //Factor //PowerExpand; epf[f_] := f //ExpandAll //PowerExpand //Factor ; pf[f_] := f //PowerExpand //Factor; exfr[f_] := ExpandAll[Numerator[f]]/ ExpandAll[Denominator[f]]; exfrt[f_] := exfr[Together[f]]; exff[f_] := Factor[PowerExpand[ExpandAll[Numerator[f]]]]/ Factor[PowerExpand[ExpandAll[Denominator[f]]]]; exfft[f_] := exff[Together[f]]; nu[f_] := Numerator[Together[f]]; Allclear := ClearAll["`*"]; solv[f_] := Solve[f][[1]]; solve[f_,a_,x_] := Module[{s,l},s = Solve[f==a,x]; l = Length[s]; Which[l==0,{},l==1,s[[1]],l==2,s[[2]],l>2,s[[1]]]]; solve[f_,x_] := solve[f,0,x]; solve0[f_,x_] := solve[f,0,x]; solve1[f_,x_] := solve[f,1,x]; solvem1[f_,x_] := solve[f,-1,x]; solveI[f_,x_] := solve[1/f,0,x]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] NoSpell --- Turn off Mathematica spelling error checker Spell --- Turn on Mathematica spelling error checker ;[s] 9:0,3;5,1;13,2;17,3;67,1;74,2;78,3;126,1;127,3;129,-1; 4:0,12,10,Courier,1,12,0,0,0;3,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;4,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Nospell := (Off[General::spell1]; Off[General::spell];); Spell := (On[General::spell1]; On[General::spell];); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Permute[l,i,j] --- Interchange list elements l[[i]] and l[[j]] ;[s] 8:0,3;5,1;19,2;23,3;51,1;58,3;64,1;71,3;73,-1; 4:0,12,10,Courier,1,12,0,0,0;3,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;4,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Permute[l_,i_,j_] := (m = l; k = m[[i]]; m[[i]] = m[[j]]; m[[j]] = k;m) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] sublist[j1,n1,j2,n2,...] = {j1 -> n1, j2 -> n2, ...} ;[s] 2:0,2;5,1;58,-1; 3:0,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) sublist[j_,n_] := {j->n}; sublist[j_,n_,k__] := Prepend[sublist[k],j->n]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] PosIntQ[y] --- True iff y is a positive integer PosIntListQ[y] --- True iff y is a list of positive integers NonposIntQ[y] --- True iff y is not a positive integer NonZeroQ[y] --- True iff y is not 0 NonZeroNumQ[y] --- True if y is a rational number not 0 False only if y is 0 NonZeroNumberQ[y] --- True iff y is a rational number not 0 NonNumberQ[y] --- True iff y is not a rational number (or integer) ;[s] 66:0,3;5,1;15,2;19,3;21,1;25,3;31,1;33,4;34,3;56,1;57,3;62,1;76,2;80,3;82,1;86,3;92,1;94,4;95,3;126,1;127,3;132,1;145,2;149,3;151,1;155,3;161,1;163,4;164,3;190,1;206,2;210,3;212,1;216,3;222,1;224,4;225,3;232,1;253,2;257,3;260,1;264,3;269,1;271,4;272,3;298,1;328,3;338,1;340,4;341,3;345,1;370,2;374,3;377,1;381,3;386,1;388,4;389,3;415,1;436,2;440,3;443,1;447,3;452,1;454,4;455,3;495,-1; 5:0,12,10,Courier,1,12,0,0,0;25,12,10,Courier,0,12,0,0,0;7,14,10,Palatino,1,12,0,0,0;26,14,10,Palatino,2,12,0,0,0;8,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) PosIntQ[y_] := If[y > 0,IntegerQ[y],False,False]; PosIntListQ[y_List] := And @@ PosIntQ /@ y; PosIntListQ[y_] := False; NonposIntQ[y_] := If[y >= 0,!IntegerQ[y],True,True]; NonZeroQ[y_] := Not[y === 0]; NonZeroNumQ[n_Integer] := n != 0; NonZeroNumQ[r_Rational] := True; NonZeroNumberQ[n_Integer] := n != 0; NonZeroNumberQ[r_Rational] := True; NonZeroNumberQ[f_] := False; NonNumberQ[n_Integer] := False; NonNumberQ[r_Rational] := False; NonNumberQ[f_] := True; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SubsetQ[l,m] --- True iff l is a subset of m StrictSubsetQ[l,m] --- True iff l is a strict subset of m ;[s] 18:0,2;4,1;17,0;20,2;23,1;27,2;34,1;36,2;52,1;55,2;59,1;78,0;81,2;84,1;88,2;95,1;97,2;120,1;123,-1; 3:2,14,10,Palatino,3,12,0,0,0;8,12,10,Courier,0,12,0,0,0;8,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SubsetQ[l_List,m_List] := And @@ (MemberQ[m,#]& /@ l); StrictSubsetQ[l_List,m_List] := And @@ Append[MemberQ[m,#]& /@ l,l=!=m]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] DeleteSupersets[l] --- Deletes elements of list l which are supersets of other elements ;[s] 6:0,2;4,1;23,0;26,2;55,1;57,2;154,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) DeleteSupersets[l_List] := Module[{s,x,m}, m = Union[Sort /@ l]; s[x_] := ! Or @@ (StrictSubsetQ[#,x]& /@ m); Select[m,s]]; SuperAppend[l_,x_] := DeleteSupersets[Append[l,x]]; SuperUnion[l__] := DeleteSupersets[Join[l]]; (* :[font = info; inactive; preserveAspect] ********************************************************************** Linear Equation Routines ********************************************************************** ;[s] 3:0,0;72,1;101,0;172,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] newvar[x,y] --- Make new variable xy by concentating x and y newvarc[x,y] --- Do same, and clear any old definitions thereof ;[s] 15:0,2;11,0;12,3;17,1;38,2;41,1;60,2;62,1;68,2;69,1;70,2;82,0;83,3;88,1;134,2;136,-1; 4:2,12,10,Courier,1,12,0,0,0;5,14,10,Palatino,2,12,0,0,0;6,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) newvar[x_,y_] := ToExpression[ToString[x] <> ToString[y]]; newvarc[x_,y_] := Module[{s,z}, Off[General::spell1]; s = ToString[x] <> ToString[y]; Clear[Evaluate[s]]; z = ToExpression[s]; On[General::spell1]; z ] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] sub0[{x,y,z}] --- {x->0,y->0,z->0} cf0[z,x,y,...] --- Set x,y,... all to 0 in z lcf[z,x,y,...] --- Give list of coefficients of linear function: z = a + b x + c y + ... gives lcf = {a, b, c,...} plcf[z,x,y,...] --- Print list of coefficients of linear function elim[f,g,z] --- Eliminate z from linear functions f = a z + b, g = c z + d, elim = a d - b c ;[s] 31:0,3;5,1;19,2;23,3;26,1;42,4;43,1;60,2;65,3;73,1;80,4;81,3;90,1;91,4;92,3;97,1;116,2;120,3;182,1;205,4;206,3;215,1;253,2;257,3;309,1;321,2;325,3;338,1;339,4;340,3;378,1;422,-1; 5:0,12,10,Courier,1,12,0,0,0;11,12,10,Courier,0,12,0,0,0;5,14,10,Palatino,1,12,0,0,0;10,14,10,Palatino,2,12,0,0,0;5,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) sub0[l__] := Rule[#,0]& /@ {l}; cf0[z_,x__] := z /. sub0[x]; lcf[z_,x__] := Join[{cf0[z,x]},cf[z,#]& /@ {x}]; plcf[z_,x__] := printlist[lcf[z,x]]; elim[f_,g_,z_] := (g /. z -> 0) cf[f,z] - (f /.z->0) cf[g,z] //ex; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] nterm[f] --- Eliminate constants from monomials ;[s] 4:0,3;5,1;13,2;17,3;54,-1; 4:0,12,10,Courier,1,12,0,0,0;1,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;2,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) nterm[f__] := f; nterm[n_Integer f_] := f; nterm[r_Rational f_] := f (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Reexpand[e] --- Re-expand list of equations discarding repetitions and trivial (zero) equations ;[s] 4:0,2;11,0;12,3;15,1;125,-1; 4:1,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,2,12,0,0,0;1,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Reexpand[e_] := Select[Union[ex[e]],NonZeroQ] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] CoefList[f,{l}] --- List of coefficients of terms in list l in f ;[s] 8:0,3;5,1;21,2;25,3;66,1;67,4;68,3;74,1;76,-1; 5:0,12,10,Courier,1,12,0,0,0;3,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) CoefList[f_,l_List] := Coefficient[f,#]& /@ l (* :[font = info; inactive; preserveAspect] ********************************************************************** Total derivatives and derivative variables ********************************************************************** ;[s] 3:0,0;72,1;136,0;207,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] dd[x][f] --- The basic format for total derivatives x is a list of independent variables, f is a function or a dependent variable dd[x,x,y,z,z,z][f] represents D2xDyD3zf dd[x,x,y,z,z,z][u] represents uxxyzzz ;[s] 24:0,3;4,1;12,2;16,3;65,1;66,3;119,1;120,3;179,1;198,3;208,1;211,4;212,5;213,1;214,5;215,1;216,4;217,5;218,1;227,3;228,1;247,3;257,1;260,5;267,-1; 6:0,20,14,Courier,34,12,0,0,0;10,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;7,14,10,Palatino,2,12,0,0,0;2,20,14,Courier,32,12,0,0,0;4,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) dd[x_][x_] := 1; dd[x_,y__][x_] := 0; dd[x__][n_Integer] := 0; dd[x__][r_Rational] := 0; dd[x_][f_+g_] := dd[x][f] + dd[x][g]; dd[x_][f_ g_] := dd[x][f] g + f dd[x][g]; dd[x_][f_^n_] := n f^(n-1) dd[x][f]; dd[x_][f_/g_] := dd[x][f]/g - f dd[x][g]/g^2; dd[x__][dd[y__][f_]] := Sort[dd[x,y]][f]; dd[x__][f_+g_] := dd[x][f] + dd[x][g]; dd[x_,y__][f_ g_] := dd[y][dd[x][f] g + f dd[x][g]]; dd[x_,y__][f_^n_] := n dd[y][f^(n-1) dd[x][f]]; dd[x_,y__][f_/g_] := dd[y][dd[x][f]/g - f dd[x][g]/g^2]; dd[x_][a_] := ( pd[x][a] + Through[(pd /@ Unt[order[a]])[a]] . (dd[x] /@ Unt[order[a]]))/; order[a] >=0; dd[x_][pd[y__][a_]] := ( pd[x,y][a] + Through[(pd /@ Unt[order[a]])[pd[y][a]]] . (dd[x] /@ Unt[order[a]]))/; order[a] >=0; dd[x_,y__][a_] := dd[y][ pd[x][a] + Through[(pd /@ Unt[order[a]])[a]] . (dd[x] /@ Unt[order[a]])]/; order[a] >=0; dd[x_,y__][pd[z__][a_]] := dd[y][ pd[x,z][a] + Through[(pd /@ Unt[order[a]])[pd[z][a]]] . (dd[x] /@ Unt[order[a]])]/; order[a] >=0; dd[x__][f_,0] := f //ex; dd[x__][f_,1] := dd[x][f] //ex; dd[x__][f_,n_Integer] := dd[x][dd[x][f,n-1]] //ex; dd[x__][l_List] := dd[x] /@ l; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] order[f] --- Gives the order of a total derivative: order[dd[x,x,y,z,z,z][u]] = 6 ;[s] 5:0,3;4,1;13,2;17,3;75,1;105,-1; 4:0,20,14,Courier,34,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;2,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) order[dd[x__][u_]] := Length[{x}]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] regformat --- The regular (default) printing format for derivatives as subscripts Ñ uxxyzzz altformat --- Alternative printing format Ñ subscripts indicate the order of derivative, so uxxyzzz becomes u2,0,1,3 when x,t,y,z are the independent variables ;[s] 18:0,1;9,2;13,3;132,1;134,4;140,3;141,1;150,2;154,3;265,1;267,4;273,3;286,1;288,4;308,3;315,1;322,3;353,4;374,-1; 5:0,20,14,Courier,34,12,0,0,0;6,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;6,14,10,Palatino,2,12,0,0,0;4,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) regformat := (Format[dd[x__][f_]] := SequenceForm[f,Subscript[SequenceForm[x]]]/; AtomQ[f]); altformat := (Format[dd[x__][f_]] := SequenceForm[f,Subscript[SequenceForm[ Infc[Count[{x},#]& /@ IndepVars]]]]/; AtomQ[f]); regformat (* :[font = info; inactive; preserveAspect] ********************************************************************** Jet coordinates ********************************************************************** ;[s] 3:0,0;72,1;92,0;163,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] flatlist --- Generic flat function flatsym --- Generic orderless flat function holdlist --- Generic hold all function ;[s] 10:0,3;2,1;11,2;15,3;41,1;49,2;53,3;89,1;98,2;102,3;130,-1; 4:0,20,14,Courier,34,12,0,0,0;3,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,1,12,0,0,0;4,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SetAttributes[flatsym, {Orderless,Flat}]; SetAttributes[flatlist, {Flat}]; SetAttributes[holdlist,HoldAll]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SetDx[x] --- Defines convenient abbreviations for the total derivatives, now abbreviated Dx, Dy, etc., except for total derivative w.r.t. t, which is Dtt so as not to conflict with Mathematica Dt ;[s] 11:0,3;4,1;13,2;17,3;106,1;114,3;166,1;169,3;180,1;184,3;236,1;239,-1; 4:0,20,14,Courier,34,12,0,0,0;5,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;5,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SetDx[x_] := (StringDx[x] = "D" <> ToString[x]; Dtotal[x] = ToExpression[StringDx[x]]; Release[Dtotal[x]] := dd[x];); SetDx[t] := (StringDx[t] = "Dtt"; Dtotal[t] = Dtt; Dtt := dd[t]); SetAttributes[SetDx,Listable]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] setderivs[u,n] --- Defines convenient abbreviations for the derivatives of u up to order n: For example, if there are two independent variable, x,y, then setderivs[u,2] will define the new variables ux, uy, uxx, uxy, uyx, uyy ;[s] 16:0,3;4,1;19,2;23,3;92,1;93,3;108,1;111,3;173,1;178,3;182,1;184,3;185,1;200,2;201,3;240,1;270,-1; 4:0,20,14,Courier,34,12,0,0,0;7,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;7,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) setderivs[u_,n_] := (Stringderivs[u] = (StringJoin @@ ToString /@ Prepend[#,u])& /@ listvart[n]; Release[ToExpression[Stringderivs[u]]] = (dd @@ # )[u]& /@ listvart[n]); SetAttributes[setderivs,Listable]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] jetset[iv,dv,k] --- Defines jet coordinates up to order k using independent variables iv and dependent variables dv jetset[iv,dv] --- Same as jetset[iv,dv,DefaultMaxDeriv] Includes the following functions: For illustration, assume iv = {x,y}, dv = {u} IndepVars --- independent variables: {x,y} DepVars --- dependent variables: {u} symvars[n] --- list of all symmetric combinations of n independent variables: symvars[2] = {{x,x},{x,y},{y,y}} symvart[n] --- list of all symmetric combinations of up to n independent variables: symvart[2] = {{x},{y},{x,x},{x,y},{y,y}} listvars[n] --- list of all combinations of n independent variables: listvars[2] = {{x,x},{x,y},{y,x},{y,y}} listvart[n] --- list of all combinations of up to n independent variables: listvart[2] = {{x},{y},{x,x},{x,y},{y,x},{y,y}} Un[n] --- list of all derivatives of order n Un[2] = {uxx, uxy, uyy} Unt[n] --- list of all derivatives up to order n Unt[2] = {u, ux, uy, uxx, uxy, uyy} Unt[n,m] --- list of all derivatives up from order n to order m Unt[2,3] = {uxx, uxy, uyy, uxxx, uxxy, uxyy, uyyy} Unx[n] --- list of all jet coordinates up to order n Unx[2] = {x, y, u, ux, uy, uxx, uxy, uyy} ;[s] 108:0,3;4,1;20,2;24,3;63,1;64,3;107,1;109,3;138,1;142,3;144,1;158,2;162,3;177,1;206,3;207,1;219,3;291,1;312,3;315,1;336,2;340,3;367,1;385,2;389,3;414,1;434,2;438,3;480,1;481,3;576,1;622,2;626,3;675,1;676,3;770,1;825,2;829,3;862,1;863,3;954,1;1008,2;1012,3;1051,1;1052,3;1143,1;1200,2;1204,3;1241,1;1242,3;1267,1;1277,0;1279,1;1282,0;1284,1;1287,0;1289,1;1300,2;1304,3;1344,1;1345,3;1370,1;1384,0;1385,1;1388,0;1389,1;1392,0;1394,1;1397,0;1399,1;1402,0;1404,1;1417,2;1421,3;1463,1;1464,3;1478,1;1479,3;1501,1;1514,0;1516,1;1519,0;1521,1;1524,0;1526,1;1529,0;1532,1;1535,0;1538,1;1541,0;1544,1;1547,0;1550,1;1561,2;1565,3;1609,1;1610,3;1635,1;1655,0;1656,1;1659,0;1660,1;1663,0;1665,1;1668,0;1670,1;1673,0;1675,1;1710,-1; 4:20,20,14,Courier,64,12,0,0,0;48,12,10,Courier,0,12,0,0,0;12,14,10,Palatino,1,12,0,0,0;28,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) jetset[iv_List,dv_List,k_Integer] := ( Off[General::spell1]; Off[General::spell]; Clear[symvars,symvart,symvarz, listvars,listvart,listvarz,Un,Unt,Unx]; If[Length[Stringvars] > 0, Clear @@ Stringvars]; IndepVars = iv; DepVars = dv; NIndepVars = Length[iv]; NDepVars = Length[dv]; StringIndepVars = ToString /@ iv; StringDepVars = ToString /@ dv; ivar[iv]; symvarz[1] = flatsym /@ iv; symvarz[n_] := (symvarz[n] = Union @@ Outer[flatsym,symvarz[n-1],iv]); symvars[n_] := (symvars[n] = symvarz[n] /. flatsym -> List); listvarz[1] = flatlist /@ iv; listvarz[n_] := (listvarz[n] = Flatten[Outer[flatlist,listvarz[n-1],iv]]); listvars[1] = listvarz[1] /. flatlist -> List; listvars[n_] := (listvars[n] = listvarz[n]/. flatlist -> List); listvart[1] = listvars[1]; listvart[n_] := (listvart[n] = Join[listvart[n-1],listvars[n]]); Un[0] = dv; Un[n_] := Un[n] = Flatten[Through[ Apply[dd,symvars[n],1][#]]& /@ dv]; Unt[0] = dv; Unt[n_] := Unt[n] = Join[Unt[n-1],Un[n]]; Unt[n_,n_] := Un[n]; Unt[n_,m_] := (Unt[n,m] = Join[Unt[n,m-1],Un[m]] ) /; m > n; Unx[0] = Join[iv,dv]; Unx[n_] := Unx[n] = Join[Unx[n-1],Un[n]]; setderivs[dv,k]; SetDx[iv]; TotalDerivs = Dtotal /@ iv; Stringvars = Join[StringIndepVars,StringDepVars, StringDx /@ iv, Flatten[Stringderivs /@ dv]]; On[General::spell1]; On[General::spell]) ; jetset[iv_,dv_List,k_] := jetset[{iv},dv,k]; jetset[iv_List,dv_,k_] := jetset[iv,{dv},k]; jetset[iv_,dv_,k_] := jetset[{iv},{dv},k]; jetset[iv_,dv_] := jetset[iv,dv,DefaultMaxDeriv]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] exponents[n,{w}] --- List of exponents in algebraic monomials of degree n using list of weights w ;[s] 7:0,3;4,1;21,2;25,3;112,1;113,3;140,1;142,-1; 4:0,20,14,Courier,34,12,0,0,0;3,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) exponents[k_,{w_,z__}] := Flatten[Table[ Prepend[#,i]& /@ exponents[k-i w,{z}], {i,0,k/w}],1]; exponents[k_,{w_}] := If[Mod[k,w]==0,{{k/w}},{}]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] monomials[m] --- List of differential monomials of weight m monomials[k,n] --- List of differential monomials of weight k in u and n derivatives in x ;[s] 16:0,3;3,1;16,2;20,3;65,1;67,3;70,1;85,2;89,3;164,1;165,3;172,1;173,3;180,1;182,3;202,1;204,-1; 4:0,20,14,Courier,34,12,0,0,0;7,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;7,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) monomials[k_,0] := Inner[Power,DepVars,#,Times]& /@ exponents[k,DepWeight] monomials[0,n_]:= {}; monomials[k_,n_]:= nterm /@ Flatten[Dx[#,n]& /@ monomials[k,0] /. Plus -> List]; monomials[m_] := (iw = IndepWeight; iq = Quotient[m,iw]; Flatten[Table[monomials[m - k iw,k],{k,0,iq}]]) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] monos[k,n] --- List of differential monomials of weight k in u and total weight m monos[n] --- List of differential monomials of total weight m arranged by weight in u ;[s] 13:0,1;12,2;16,3;91,1;92,3;99,1;100,3;121,1;133,2;137,3;188,1;206,3;230,1;232,-1; 4:0,20,14,Courier,34,12,0,0,0;6,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;5,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) monos[k_,m_] := If[Mod[m-k,IndepWeight] ==0, monomials[k,(m - k)/IndepWeight],{},{}]; monos[m_] := Table[monos[k,m+#], {k,minDepWeight,m+#}]& /@ DepWeight; (* :[font = info; inactive; preserveAspect] ********************************************************************** Functions and partial derivatives ********************************************************************** ;[s] 3:0,0;72,1;110,0;181,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] pd[x][f] --- The basic format for partial derivatives x list of independent and dependent variables and their derivatives, f a function pd[x,x,u,u,u][f] represents ¶6f/¶x2¶u3 pd[x,x,dd[y,z][u]][f] represents ¶6f/¶x2¶uyz ;[s] 26:0,3;4,1;12,2;16,3;68,1;69,3;153,1;154,3;184,1;201,3;211,1;214,4;215,1;219,4;220,1;222,4;224,3;239,1;261,3;271,1;274,4;275,1;279,4;280,1;282,5;284,4;286,-1; 6:0,12,10,Courier,1,12,0,0,0;11,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;7,14,10,Palatino,2,12,0,0,0;6,20,14,Courier,32,12,0,0,0;1,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SetAttributes[pd,Orderless]; Format[pd[x__][f_]] := SequenceForm[f,Subscript[SequenceForm[x]]]; pd[x_][x_] := 1; pd[x__][y_Symbol] := 0 /; NonposIntQ[order[y]] pd[x__][n_Integer] := 0; pd[x__][r_Rational] := 0; pd[x_][f_+g_] := pd[x][f] + pd[x][g]; pd[x_][f_ g_] := pd[x][f] g + f pd[x][g]; pd[x_][f_^n_] := n f^(n-1) pd[x][f]; pd[x_][f_/g_] := pd[x][f]/g - f pd[x][g]/g^2; pd[x__][pd[y__][f_]] := pd[x,y][f]; pd[x__][f_+g_] := pd[x][f] + pd[x][g]; pd[x_,y__][f_ g_] := pd[y][pd[x][f] g + f pd[x][g]]; pd[x_,y__][f_^n_] := n pd[y][f^(n-1) pd[x][f]]; pd[x_,y__][f_/g_] := pd[y][pd[x][f]/g - f pd[x][g]/g^2]; pd[y___,dd[x__][g_]][f_] := 0 /; Length[{x}] > order[f]; pd[x__][dd[y__][f_]] := 0; pd[x__][l_List] := pd[x] /@ l; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] nodep[f,x] --- Make f not depend on variable x dep[f,x] --- Make f depend on variable x again ;[s] 14:0,1;16,3;19,2;28,1;29,3;30,2;54,1;73,3;76,2;85,1;86,3;87,2;108,1;111,2;117,-1; 4:0,12,10,Courier,1,12,0,0,0;5,12,10,Courier,0,12,0,0,0;5,14,10,Palatino,2,12,0,0,0;4,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) nodep[x_,x_] := Null; nodep[f_,x_] := (f /: pd[y___,x,z___][f] = 0;) SetAttributes[nodep,Listable]; dep[f_,x_] := (f /: pd[y___,x,z___][f] = .;) SetAttributes[dep,Listable]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] ivar[x] --- Make x an independent variable ;[s] 7:0,3;5,1;12,2;16,3;25,1;26,4;27,3;52,-1; 5:0,12,10,Courier,1,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) ivar[x_] := (x /: dd[x][x] = 1; x /: dd[y__][x] = 0;); SetAttributes[ivar,Listable]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] depn[f,n] --- Make f depend on derivative variables up to order n ;[s] 8:0,3;4,1;13,2;17,3;27,1;28,4;29,3;75,1;77,-1; 5:0,12,10,Courier,1,12,0,0,0;3,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) depn[f_,n_] := (Clear[f]; order[f] = n; f/: pd[y___,dd[x__][g_]][f] := 0 /; Length[{x}] > n;); SetAttributes[depn,Listable]; dep0[f_] := depn[f,0]; dep1[f_] := depn[f,1]; dep2[f_] := depn[f,2]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] depu[f] --- Make f depend only on indep. and dependent variables ;[s] 7:0,2;5,1;13,3;16,2;26,1;27,3;28,2;75,-1; 4:0,12,10,Courier,1,12,0,0,0;2,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;2,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) depu[f_] := depn[f,0]; SetAttributes[depu,Listable]; depu[f_,g__] := (depn[f,0]; depu[g]); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] depx[f] --- Make f depend only on independent variables ;[s] 7:0,3;5,1;12,2;16,3;26,1;27,4;28,3;66,-1; 5:0,12,10,Courier,1,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) depx[f_] := (depu[f]; f/: pd[y___,u_][f] := 0 /; MemberQ[DepVars,u];); SetAttributes[depx,Listable]; depx[f_,g__] := (depx[f]; depx[g]); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] const[a,b,...] --- Make a,b,... constant ;[s] 7:0,3;5,1;19,2;23,3;33,1;40,4;41,3;51,-1; 5:0,12,10,Courier,1,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) const[f_] := (pd[y___][f] := 0;order[f] = 0;); const[f_,g__] := (const[f]; const[g];); SetAttributes[const,Listable]; constant = const; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] constants[a,n] --- Construct the list {a1,a2,a3,...,an} of constants constants[a,n,m] --- Construct the list {am1,am2,...,amn} of constants ;[s] 13:0,3;5,1;19,2;23,3;47,1;64,4;65,3;138,1;154,2;158,3;182,1;199,4;200,3;314,-1; 5:0,12,10,Courier,1,12,0,0,0;4,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;5,14,10,Palatino,2,12,0,0,0;2,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) constants[a_,n_] := Module[{c,s}, Off[General::spell1];Off[General::spell]; c = Table[s = ToString[a] <> ToString[i]; Clear[s]; ToExpression[s],{i,n}]; const[c];On[General::spell1]; On[General::spell]; c]; constants[a_,n_,m__] := Module[{c,s}, Off[General::spell1];Off[General::spell]; c = Flatten[Table[s = ToString[a] <> StringJoin @@ ToString /@ {m} <> ToString[i]; Clear[Evaluate[s]]; ToExpression[s],{i,n}]]; const[c]; On[General::spell1]; On[General::spell];c]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] setdegree[f,n] - Make f a polynomial of degree at most n ;[s] 8:0,3;2,1;16,2;18,3;28,1;29,4;30,3;64,1;66,-1; 5:0,12,10,Courier,1,12,0,0,0;3,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) setdegree[f_,n_] := (deg[f] = n; pd[x__][f] := 0 /; Length[{x}] > n;) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] dependencies[f] --- List of current variables f depends on printdep[f] --- Print current variables f depends on ;[s] 13:0,3;4,1;19,2;23,3;54,1;55,4;56,3;71,1;83,2;87,3;116,1;117,4;118,3;131,-1; 5:0,12,10,Courier,1,12,0,0,0;4,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;5,14,10,Palatino,2,12,0,0,0;2,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect] *) dependencies[a_] := If[order[a] >=0, pdaQ[x_] := NonZeroQ[pd[x][a]]; Select[Unx[order[a]],pdaQ],IndepVars,IndepVars] (* :[font = input; initialization; preserveAspect; endGroup] *) printdep[a_] := Module[{d = dependencies[a]}, If[Length[d] == 0, Print[a," is constant"], Print[a," depends on ",d]]] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] constantQ[f] --- True if and only if f is constant ;[s] 10:0,3;4,1;16,2;21,3;23,1;27,2;29,3;48,1;49,4;50,3;64,-1; 5:0,12,10,Courier,1,12,0,0,0;3,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;4,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) constantQ[n_Integer] := True; constantQ[r_Rational] := True; constantQ[z_Plus] := constantQ /@ And @@ z; constantQ[z_Times] := constantQ /@ And @@ z; constantQ[x_^n_] := constantQ[x]; constantQ[x_] := (Length[dependencies[x]] == 0); (* :[font = info; inactive; preserveAspect] ********************************************************************** Vector fields ********************************************************************** ;[s] 3:0,0;72,1;90,0;161,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] vq[q,u][z] - Apply prolongation of evolutionary vector field q¶u to z Both q and u can be lists, of the same length ;[s] 12:0,3;10,0;13,2;64,3;66,4;67,2;72,1;73,2;116,3;118,2;122,3;126,2;159,-1; 5:1,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,0,12,0,0,0;5,14,10,Palatino,2,12,0,0,0;4,12,10,Courier,0,12,0,0,0;1,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) ClearAll[vq]; vq[q_List,u_List][z_] := Plus@@Through[Inner[vq,q,u,List][z]]; vq[q_List] := vq[q,DepVars]; vq[q_,u_][u_] := q; vq[q_,u_][y_Symbol] := 0 /; NonposIntQ[order[y]]; vq[q_,u_][n_Integer] := 0; vq[q_,u_][r_Rational] := 0; vq[q_,u_][f_+g_] := vq[q,u][f] + vq[q,u][g]; vq[q_,u_][f_ g_] := vq[q,u][f] g + f vq[q,u][g]; vq[q_,u_][f_^n_] := n f^(n-1) vq[q,u][f]; vq[q_,u_][f_/g_] := vq[q,u][f]/g - f vq[q,u][g]/g^2; vq[q_,u_][dd[x__][u_]] := dd[x][q]; vq[q_,u_][dd[y__][v_]] := 0; vq[q_,u_][a_] := q pd[u][a] /; order[a] ==0; vq[q_,u_][a_] := ( q pd[u][a] + Through[(pd /@ Through[Apply[dd,symvars[order[a]] ,1][u]])[a]] . Through[Apply[dd, symvars[order[a]],1][q]]) /; order[a] > 0; vq[q_,u_][pd[y__][a_]] := q pd[u,y][a] /; order[a] ==0; vq[q_,u_][pd[y__][a_]] := ( q pd[u,y][a] + Through[(pd /@ Through[Apply[dd,symvars[order[a]],1][u]])[pd[y][a]]] . Through[Apply[dd,symvars[order[a]],1][q]])/; order[a] > 0; vq[q__][l_List] := vq[q] /@ l; vq[q_] := vq[q,DepVars[[1]]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] char[a,f] or char[a,f,x,u] --- Characteristic of vector field a¶x + f¶u a and x can be lists, of the same length f and u can be lists, of the same length ;[s] 16:0,2;27,0;31,1;106,2;108,3;109,2;114,3;115,1;156,2;158,1;162,2;166,1;238,2;240,1;244,2;248,1;316,-1; 4:1,12,10,Courier,1,12,0,0,0;6,14,10,Palatino,2,12,0,0,0;7,12,10,Courier,0,12,0,0,0;2,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) char[a_List,b_List,x_List,u_List] := b - ((a . #)& /@ (Through[(dd /@ x)[#]]& /@ u)); char[a_List,b_List] := char[a,b,IndepVars,DepVars]; char[a_,b_,x_,u_] := char[{a},{b},{x},{u}] char[a_,b_] := char[{a},{b},IndepVars,DepVars];; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] vp[a,f][z] or vp[a,f,x,u][z] --- Apply prolongation of (generalized) vector field v = a¶x + f¶u to z a and x can be lists, of the same length f and u can be lists, of the same length ;[s] 22:0,3;11,2;16,3;31,0;36,2;89,3;93,2;94,3;96,4;97,3;102,4;103,2;109,1;110,2;153,3;155,2;159,3;163,2;200,3;202,2;206,3;210,2;243,-1; 5:1,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,0,12,0,0,0;9,14,10,Palatino,2,12,0,0,0;9,12,10,Courier,0,12,0,0,0;2,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) vp[a_List,b_List,x_List,u_List][z_] := vq[char[a,b,x,u],u][z] + a . (dd[#][z]& /@ x) //ex; vp[a_List,b_List][z_] := vp[a,b,IndepVars,DepVars][z]; vp[a_,b_,x_,u_][z_] := vp[{a},{b},{x},{u}][z]; vp[a_,b_][z_] := vp[{a},{b},IndepVars,DepVars][z]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] vcom[q,s] --- Commutator [q¶u, s¶u] of evolutionary vector fields vcom[a,f,b,g] --- Commutator [a¶x + f¶u, b¶x + g¶u] of (generalized) vector fields ;[s] 22:0,2;10,0;14,1;30,2;33,3;34,2;38,3;39,2;40,1;72,2;86,0;90,1;106,2;109,3;110,2;115,3;116,2;120,3;121,2;126,3;127,2;128,1;200,-1; 4:2,12,10,Courier,1,12,0,0,0;4,14,10,Palatino,2,12,0,0,0;10,12,10,Courier,0,12,0,0,0;6,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) vcom[q_,s_] := vq[q][s] - vq[s][q] //ex; vcom[a_,f_,b_,g_] := vp[a,f][{b,g}] - vp[b,g][{a,f}] //ex; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] vfstr --- For printing vector field coefficients ;[s] 4:0,2;4,1;10,0;13,2;54,-1; 3:1,14,10,Palatino,3,12,0,0,0;1,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) vfstr[a_,x_] := SequenceForm["¶",Subscript[ToString[x]], " : ", ToString[a]]; vfstr[l_List,m_List] := Module[{v}, Thread[v[l,m]] /. v -> vfstr]; (* :[font = info; inactive; preserveAspect] ********************************************************************** :[font = message; inactive; wordwrap; dontPreserveAspect] Variational Derivatives ;[s] 3:0,1;25,2;26,0;31,-1; 3:1,14,11,Courier,0,14,65535,0,0;1,17,13,Courier,0,17,65535,0,0;1,13,10,Courier,0,13,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] derivs[f,u] --- List of derivatives of u occurring in f. ;[s] 7:0,3;11,0;15,2;42,3;43,2;58,1;59,2;61,-1; 4:1,14,10,Palatino,1,12,0,0,0;1,14,10,Palatino,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;2,12,10,Courier,0,12,0,0,0; :[font = input; initialization; wordwrap; dontPreserveAspect; endGroup] *) derivs[f_,u_] := Union[ Cases[Level[f,{-3}],dd[__][u]], Cases[Level[f,{-2}],dd[__][u]]]; derivs[f_List,l_List] := derivs[#,l]& /@ f; derivs[f_,l_List] := Flatten[derivs[f,#]& /@ l]; derivs[f_List,u_] := derivs[#,u]& /@ f; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] maxorder[f,u] --- Maximal order of derivatives of u occurring in f. ;[s] 7:0,3;13,0;17,2;53,3;54,2;69,1;70,2;72,-1; 4:1,14,10,Palatino,1,12,0,0,0;1,14,10,Palatino,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;2,12,10,Courier,0,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) maxorder[f_,u_] := Max[0,order /@ derivs[f,u]]; maxorder[f_] := Max[0,order /@ derivs[f,DepVars]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] dorder[n] --- Order of derivative indexed by n dsign[n] --- Sign of derivative indexed by n,i.e. (-1)^dorder[n] ;[s] 10:0,0;1,2;10,3;14,1;50,2;62,3;66,1;100,2;102,1;109,2;125,-1; 4:1,12,10,Courier,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;4,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; wordwrap; dontPreserveAspect; endGroup] *) dindex[x_]:=1; dindex[{x_,n_}]:=n; dorder[n__] := Plus @@ dindex /@ {n}; dsign[n__] := (-1)^dorder[n]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Euler[f] --- Variational derivative of fwith respect to current dependent variables (Euler operator) Euler[f,u] --- Variational derivative of fwith respect to u ;[s] 13:0,1;10,3;13,0;14,2;43,1;44,2;143,1;153,0;154,3;159,2;187,1;188,2;204,1;207,-1; 4:2,12,10,Courier,1,12,0,0,0;5,12,10,Courier,0,12,0,0,0;4,14,10,Palatino,2,12,0,0,0;2,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Eulerterm[dd[x__][u_],f_] := (-1)^Length[{x}] dd[x][pd[dd[x][u]][f]]; Euler[f_,u_] := pd[u][f] + Plus @@ (Eulerterm[#,f]& /@ derivs[f,u] )//ex; Euler[f_] := Module[{ez}, ez = Euler[f,#]& /@ DepVars; If[Length[DepVars] ==1,ez[[1]],ez]] (* :[font = info; inactive; preserveAspect] ********************************************************************** :[font = message; inactive; wordwrap; dontPreserveAspect] Symmetry Group Computations ;[s] 3:0,1;29,2;30,0;35,-1; 3:1,14,11,Courier,0,14,65535,0,0;1,17,13,Courier,0,17,65535,0,0;1,13,10,Courier,0,13,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] DeterminingEquations[p == q] Find determining equations for symmetry of p.d.e. p = q where p is a single jet variable, or solved system Symmetry vector field is defined by Vsym and depends on variables VfVars DetEq[p == q] Same as DeterminingEquations but vector field Vsym is set to generic point transformation at beginning ;[s] 14:0,2;28,1;88,2;94,1;115,2;117,1;208,2;224,1;251,2;271,1;290,2;311,1;343,2;348,1;400,-1; 3:0,12,10,Courier,1,12,0,0,0;7,14,10,Palatino,2,12,0,0,0;7,12,10,Courier,0,12,0,0,0; :[font = input; initialization; preserveAspect] *) SubtractEq[p_ == q_] := p - q; SubstituteEq[p_ == q_] := p -> q; (* :[font = input; initialization; preserveAspect; endGroup] *) DeterminingEquations[e_] := Module[{ez,eqorder,sz,vz,wz}, ez = SubtractEq[e]; eqorder = maxorder[ez]; vz = Vsym[ez]; sz = SubstituteEq[e]; wz = vz /. sz //ex; coefs[wz, Complement[Unx[eqorder],VfVars]]]; DeterminingEquations[l_List] := Module[{ez,eqorder,sz,vz,wz}, ez = SubtractEq /@ l; eqorder = Max[maxorder /@ ez]; vz = Vsym /@ ez; sz = SubstituteEq /@ l; wz = ExpandAll[vz /. sz]; Flatten[coefs[#, Complement[Unx[eqorder],VfVars]]& /@ wz]]; DetEq[e_] := (Vset; DeterminingEquations[e]); (* :[font = info; inactive; preserveAspect] ********************************************************************** :[font = message; inactive; wordwrap; dontPreserveAspect] Solution of Simple Overdetermined Systems of Partial Differential Equations ;[s] 3:0,1;89,2;90,0;95,-1; 3:1,14,11,Courier,0,14,65535,0,0;1,17,13,Courier,0,17,65535,0,0;1,13,10,Courier,0,13,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] coefs[f,l] --- List of nonzero coefficients of variables in list l occuring in f. ;[s] 9:0,3;10,0;11,4;14,0;15,2;67,3;68,2;125,1;126,2;128,-1; 5:2,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) coefs[f_,l_] := Select[Flatten[CoefficientList[ Numerator[Together[f]],l]],NonZeroQ]; coefs[0,l_] := {}; coefs[f_,{}] := {f}; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] coefn[f,m,n] --- List of coefficients of derivatives of u of orders between m and n occuring in f. ;[s] 13:0,3;12,0;13,4;16,0;17,2;58,3;59,2;121,3;124,2;130,3;133,2;147,1;148,2;150,-1; 5:2,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,0,12,0,0,0;5,14,10,Palatino,2,12,0,0,0;4,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) coefn[f_,m_,n_] := coefs[f,Unt[m,n]] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] updef[f,g] --- Make f have the same definitions as g. ;[s] 8:0,3;10,0;11,4;16,2;23,3;24,2;58,1;59,2;61,-1; 5:1,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) updef[f_,g_] := (UpValues[f] = UpValues[g] /. g -> f;) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] solvealg[e,z] --- Solve simple algebraic equation e = c z + g = 0 for z = -g/c (provided no derivatives of z occur in g) Also appends any differential equations which result from dependencies of z to PDEqs e.g. if z doesn't depend on x then ¶x(g/c) = 0 ;[s] 23:0,2;13,0;14,3;19,1;52,2;86,1;99,2;115,1;142,2;145,1;159,2;167,1;260,2;262,1;266,2;274,1;300,2;302,1;321,2;324,1;330,2;332,4;333,2;344,-1; 5:1,12,10,Courier,1,12,0,0,0;9,14,10,Palatino,2,12,0,0,0;11,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;1,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) solvealg[e_, z_Symbol] := Module[{s, c = cf[ex[e], z], v, w}, If[c === 0,,If[constantQ[c], s = z - ex[e/c]; If[FreeQ[s, z], v = Complement[Uvars, dependencies[z]]; Unknowns = DeleteCases[Unknowns,z]; z = s; PDEqs = Union[PDEqs, Select[pd[#][s]& /@ v, NonZeroQ]]]]]]; solvealg[e_, z_List] := (solvealg[e, #] & ) /@ z (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] solvepde[e] --- Solve simple linear p.d.e. e = 0 Possible equations: u = 0, ux = 0, uxx...x = 0, or multiples thereof. (If the multiple is not obviously constant, the program gives a warning message.) ;[s] 12:0,2;11,0;12,4;17,1;45,2;51,1;111,2;120,3;121,2;129,3;135,2;145,1;264,-1; 5:1,12,10,Courier,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;5,12,10,Courier,0,12,0,0,0;2,20,14,Courier,64,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) solvepde[0] := {}; solvepde[f_Symbol] := If[MemberQ[Unknowns,f], Unknowns = DeleteCases[Unknowns,f]; f = 0, Print["** Possible Contradiction:", f,"== 0 **"]; Print[" "]]; solvepde[pd[x_][f_]] := nodep[f,x]; solvepde[pd[x__][f_]] := Module[ {x0 = First[{x}],l = Length[{x}]-1,fz}, fz = newvarc[f,#]& /@ Range[0,l]; (order[#] = order[f])& /@ fz; updef[#,f]& /@ fz; nodep[#,x0]& /@ fz; Unknowns = Flatten[Unknowns /. f -> fz]; f = fz . (x0 ^ Range[0,l])] /; Equal @@ {x}; solvepde[n_Integer z_] := solvepde[z]; solvepde[r_Rational z_] := solvepde[z]; solvepde[z_Times]:= Module[{c,lz = List @@ z,zx = z//ex}, If[zx===z, c = Join[Flatten[Cases[lz,#]& /@ Unknowns], Flatten[Cases[lz,pd[__][_]]& /@ lz]]; If[Length[c]==0, Print["** A Possible Contradiction:", z," == 0 **"]; Print[" "], Print["Assuming equation ",z, " == 0 has non-zero coefficient."];Print[" "]; solvepde[c[[1]]]], solvepde[zx]]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SolvPDE --- Try to solve system of partial differential equations PDEqs involving original unknown functions OrigUnknowns with their string names in StringUnknowns, depending on variables Uvars current unknown functions Unknowns by recursively applying solvealg and solvepde until nothing more happens. Print intermediate results if PrintPDEQ = True ;[s] 21:0,2;7,0;8,3;11,1;74,2;80,1;132,2;152,1;179,2;195,1;234,2;240,1;243,2;244,1;280,2;290,1;320,2;332,1;335,2;347,1;420,2;438,-1; 4:1,12,10,Courier,1,12,0,0,0;9,14,10,Palatino,2,12,0,0,0;10,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SolvPDE := Module[{pdeqs = Null}, While[pdeqs =!= PDEqs, OrigUnknowns = OrigUnknowns //ex; If[PrintPDEQ, PrintPDE]; pdeqs = PDEqs; solvealg[#,Unknowns]& /@ PDEqs; PDEqs = Reexpand[PDEqs]; solvepde /@ PDEqs; PDEqs = Reexpand[PDEqs]; Vvars = Union @@ (dependencies /@ Unknowns); Wvars = Complement[Uvars,Vvars]; Uvars = Vvars; PDEqs = Flatten[coefs[#,Wvars]& /@ PDEqs]]; If[!PrintPDEQ,PrintPDE]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SolvePDE[e,u,s,v] --- Try to solve system of partial differential equations e involving unknown functions u, with their string names in s, depending on variables v by recursively applying solvealg and solvepde until nothing more happens. ;[s] 16:0,2;17,0;18,3;21,1;85,2;87,1;129,2;139,1;166,2;169,1;208,2;210,1;243,2;255,1;258,2;270,1;305,-1; 4:1,12,10,Courier,1,12,0,0,0;7,14,10,Palatino,2,12,0,0,0;7,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SolvePDE[e_,u_,s_,v_] := (Unknowns = OrigUnknowns = u; StringUnknowns = s; Uvars = v; PDEqs = e; SolvPDE); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] DiffPDE --- Prolong system of differential equations PDEqs by appending all nontrivial equations obtained by differentiating them once ;[s] 6:0,2;7,0;8,3;11,1;59,2;65,1;152,-1; 4:1,12,10,Courier,1,12,0,0,0;2,14,10,Palatino,2,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) DiffPDE := (PDEqs = Union[PDEqs, Select[Flatten[ Through /@ ((pd /@ Uvars) /@ PDEqs)],NonZeroQ]]) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SolveDiff --- Try to solve system of differential equations PDEqs by differentiating them ;[s] 6:0,2;9,0;10,3;13,1;66,2;79,1;106,-1; 4:1,12,10,Courier,1,12,0,0,0;2,14,10,Palatino,2,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SolveDiff := (DiffPDE; SolvPDE) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] PrintPDE --- Print system of differential equations PDEqs ;[s] 5:0,2;8,0;9,3;12,1;58,2;65,-1; 4:1,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,2,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) PrintPDE := (If[Length[PDEqs] > 0, Print[" Current Differential Equations"]; prstars; printlist[PDEqs]; prstars;, Print["--> Differential Equations Solved! <--"]; prstars]; Print["Current Values for Unknowns"]; prstars; printeq[StringUnknowns,OrigUnknowns]; prstars; printdep /@ Unknowns; prstars;) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SolveDeterminingEquations[e] --- Compute and try to solve determining equations for system e Symmetry vector field is defined by Vsym ;[s] 8:0,2;28,0;29,3;32,1;98,2;104,1;152,2;158,1;161,-1; 4:1,12,10,Courier,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;3,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SolveDeterminingEquations[e_] := (deqs = DeterminingEquations[e]; SolvePDE[deqs,VfCoefs,VfString,VfVars]) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SolveDetEq[e] --- Same as SolveDeterminingEquations but vector field Vsym is set to generic point transformation at beginning ;[s] 8:0,2;13,0;14,3;17,1;35,2;61,1;93,2;98,1;151,-1; 4:1,12,10,Courier,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0;3,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SolveDetEq[e_] := (Vset;SolveDeterminingEquations[e]) (* :[font = info; inactive; preserveAspect] ********************************************************************** Noncommutative Multiplication ********************************************************************** ;[s] 3:0,0;72,1;106,0;177,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Replace Mathematica non-commutative multiplication u ** v or NonCommutativeMultiply[u,v] by package non-commutative multiplication NonCommutative[u,v] formatted as u*v which allows for differentiations, etc. Occasionally, for technical reasons, a non Flat version is required, called NonCommutation[u,v] ;[s] 21:0,2;11,3;22,0;23,2;98,1;106,2;108,1;110,2;113,1;144,2;230,1;251,2;268,1;272,2;275,1;279,2;366,1;370,2;444,1;464,2;467,1;470,-1; 4:1,14,10,Palatino,3,12,0,0,0;9,12,10,Courier,0,12,0,0,0;10,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect] *) Unprotect[NonCommutativeMultiply]; NonCommutativeMultiply = NonCommutative; SetAttributes[NonCommutative,Flat]; Format[NonCommutative[x__]] := x; Format[NonCommutative[x_,y__]] := Infix[NonCommutative[x,y],"*"]; NonCommutative[a___,0,c___] := 0; NonCommutative[a___,b_ + c_,d___] := NonCommutative[a,b,d] + NonCommutative[a,c,d]; NonCommutative[a___,n_Integer b_,c___] := n NonCommutative[a,b,c]; NonCommutative[a___,r_Rational b_,c___] := r NonCommutative[a,b,c]; NonCommutative[a___,c_ b_,d___] := c NonCommutative[a,b,d] /; order[c] == 0; NonCommutative[a___,c_^n_ b_,d___] := c^n NonCommutative[a,b,d] /; order[c] == 0; (* :[font = input; initialization; preserveAspect] *) NonCommutation[f_ + g_] := NonCommutation[f] + NonCommutation[g]; NonCommutation[n_Integer f_] := n NonCommutation[f]; NonCommutation[r_Rational f_] := r NonCommutation[f]; NonCommutation[f_ g_] := NonCommutative[ NonCommutation[f], NonCommutation[g]]; NonCommutation[f_^n_] := NonCommutative @@ Table[NonCommutation[f],{n}]; NonCommutation[f_] := f; (* :[font = input; initialization; preserveAspect] *) dd[x_][NonCommutatively[f_,g_]] := NonCommutatively[dd[x][f],g] + NonCommutatively[f,dd[x][g]]; dd[x_][NonCommutatively[f_,g__]] := NonCommutatively[dd[x][f],g] + NonCommutatively[f,dd[x][NonCommutatively[g]]] ; dd[x_,y__][NonCommutatively[f__]] := dd[y][dd[x][NonCommutatively[f]]]; dd[x__][NonCommutative[f__]] := dd[x][NonCommutatively[f]] /. NonCommutatively -> NonCommutative; (* :[font = input; initialization; preserveAspect] *) pd[x_][NonCommutatively[f_,g_]] := NonCommutatively[pd[x][f],g] + NonCommutatively[f,dd[x][g]]; pd[x_][NonCommutatively[f_,g__]] := NonCommutatively[pd[x][f],g] + NonCommutatively[f,pd[x][NonCommutatively[g]]] ; pd[x_,y__][NonCommutatively[f__]] := pd[y][pd[x][NonCommutatively[f]]]; pd[x__][NonCommutative[f__]] := pd[x][NonCommutatively[f]] /. NonCommutatively -> NonCommutative; (* :[font = input; initialization; preserveAspect; endGroup] *) vq[q_,u_][NonCommutatively[f_,g_]] := NonCommutatively[vq[q,u][f],g] + NonCommutatively[f,vq[q,u][g]]; vq[q_,u_][NonCommutatively[f_,g__]] := NonCommutatively[vq[q,u][f],g] + NonCommutatively[f,vq[q,u][NonCommutatively[g]]]; vq[q_,u_][NonCommutative[f__]] := vq[q,u][NonCommutatively[f]] /. NonCommutatively -> NonCommutative; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] NCpermutations[m] --- List of different non-commutative versions of a term. For example NCpermutations[a b^2] = {a*b*b,b*a*b,b*b*a} ;[s] 5:0,3;3,1;21,2;25,3;171,1;215,-1; 4:0,20,14,Courier,34,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;2,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) NCpermutations[f_[x__]] := Permutations[NonCommutation[f[x]]]; NCpermutations[f_] := f (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] NCmonomials[m] --- List of non-commutative differential monomials of weight m NCmonomials[k,n] --- List of non-commutative differential monomials of weight k in u and n derivatives in x NCmonos[k,m] --- List of non-commutative differential monomials of weight k in u and total weight m ;[s] 20:0,1;16,2;20,3;133,1;153,2;157,3;269,1;270,3;277,1;278,3;285,1;287,3;307,1;323,2;327,3;437,1;438,3;445,1;446,3;467,1;469,-1; 4:0,20,14,Courier,34,12,0,0,0;9,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,1,12,0,0,0;8,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) NCmonomials[n__] := Flatten[NCpermutations /@ monomials[n]]; (* :[font = input; initialization; preserveAspect; endGroup] *) NCmonos[k_,m_] := If[Mod[m-k,IndepWeight] ==0, NCmonomials[k,(m - k)/IndepWeight],{},{}]; NCmonos[m_] := Table[NCmonos[k,m+#], {k,minDepWeight,m+#}]& /@ DepWeight; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] NCtranspose[f] --- Transpose noncommutative expression f ;[s] 4:0,1;16,2;20,3;60,1;63,-1; 4:0,20,14,Courier,34,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;1,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) NCtranspose[f_ + g_] := NCtranspose[f] + NCtranspose[g]; NCtranspose[n_Integer f_] := n NCtranspose[f]; NCtranspose[r_Rational f_] := r NCtranspose[f]; NCtranspose[f_] := Reverse[f]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Commutative[f] --- Replace noncommutative multiplication by commutative multiplication ;[s] 3:0,1;16,2;20,3;150,-1; 4:0,20,14,Courier,34,12,0,0,0;1,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;1,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Commmutative[f_] := f /. NonCommutative -> Times; (* :[font = info; inactive; preserveAspect] ********************************************************************** Classification of integrable evolution equations via higher order symmetries ********************************************************************** ;[s] 3:0,0;72,1;160,0;231,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect] Q --- evolution equation Qw --- weight of equation Qm --- monomials in equation collected by weight Qmf --- monomials in equation Qcf --- list of coefficients in equation Qc --- list of coefficients in equation arranged by degree Ql --- list of components in equation arranged by degree Qa --- remaining unnormalized coefficients in equation S --- Symmetry Sw --- weight of symmetry Sm --- monomials in symmetry collected by weight Smf --- monomials in symmetry Scf --- list of coefficients in symmetry Sc --- list of coefficients in symmetry arranged by degree Sl --- list of components in symmetry arranged by degree Sb --- remaining unnormalized coefficients in symmetry QSw --- weight of commutator QSm --- monomials in commutator ;[s] 55:0,1;2,0;5,2;26,1;32,0;35,2;56,1;62,0;65,2;109,1;115,0;118,2;142,1;149,0;152,2;187,1;193,0;196,2;250,1;256,0;259,2;311,1;317,0;320,2;371,1;376,0;379,2;390,1;396,0;399,2;420,1;426,0;429,2;473,1;480,0;483,2;507,1;514,0;517,2;552,1;558,0;561,2;615,1;621,0;624,2;676,1;682,0;685,2;736,1;743,0;746,2;769,1;776,0;779,2;805,1;809,-1; 3:18,14,10,Palatino,3,12,0,0,0;19,12,10,Courier,0,12,0,0,0;18,14,10,Palatino,2,12,0,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; rightWrapOffset = 558; startGroup] CQS[m,n] --- Initialize to general commutative evolution equation of order m and general symmetry of order n CQS[m,q,n] --- Initialize to commutative evolution equation q of order m and general symmetry of order n CQS[m,n,r] --- Initialize to general commutative evolution equation of order m and symmetry r of order n CQS[m,q,n,r] --- Initialize to commutative evolution equation q of order m and symmetry r of order n NCQS[m,n] --- Initialize to general non-commutative evolution equation of order m and general symmetry of order n NCQS[m,q,n] --- Initialize to non-commutative evolution equation q of order m and general symmetry of order n NCQS[m,n,r] --- Initialize to general non-commutative evolution equation of order m and symmetry r of order n NCQS[m,q,n,r] --- Initialize to non-commutative evolution equation q of order m and symmetry r of order n CQ[m] --- Initialize to general commutative evolution equation of order m CQ[m,q] --- Initialize to commutative evolution equation q of order m NCQ[m] --- Initialize to general non-commutative evolution equation of order m NCQ[m,q] --- Initialize to non-commutative evolution equation q of order m CS[n] --- Initialize to general commutative symmetry of order n CS[n,r] --- Initialize to commutative symmetry r of order n NCS[n] --- Initialize to general non-commutative symmetry of order n NCS[n,r] --- Initialize to non-commutative symmetry r of order n In all cases, if q or r is a list of positive integers, then take terms of those degrees of homogeneity in the dependent vars from the general evolution equation or symmetry. ;[s] 97:0,1;9,0;12,2;80,1;103,2;136,1;149,0;152,2;202,1;225,2;237,1;239,2;272,1;285,0;288,2;356,1;379,2;394,1;396,2;408,1;423,0;426,2;476,1;499,2;511,1;513,2;528,1;530,2;542,1;554,0;557,2;629,1;652,2;685,1;699,0;702,2;756,1;779,2;791,1;793,2;826,1;840,0;843,2;902,1;923,2;936,1;938,2;953,1;955,2;967,1;983,0;986,2;1040,1;1063,2;1075,1;1077,2;1092,1;1094,2;1106,1;1114,0;1117,2;1185,1;1196,0;1199,2;1249,1;1251,2;1263,1;1272,0;1275,2;1347,1;1358,0;1361,2;1415,1;1417,2;1429,1;1437,0;1440,2;1498,1;1508,0;1511,2;1551,1;1553,2;1565,1;1574,0;1577,2;1639,1;1650,0;1653,2;1697,1;1699,2;1711,1;1714,2;1733,1;1735,2;1740,1;1741,2;1905,1;1909,-1; 3:16,14,10,Palatino,3,12,0,0,0;41,12,10,Courier,0,12,0,0,0;40,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) CQS[m_Integer,n_Integer] := (CQ[m]; CS[n]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := monos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); CQS[m_Integer,q_,n_Integer] := (CQ[m,q]; CS[n]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := monos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); CQS[m_Integer,n_Integer,r_] := (CQ[m]; CS[n,r]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := monos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); CQS[m_Integer,q_,n_Integer,r_] := (CQ[m,q]; CS[n,r]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := monos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) NCQS[m_Integer,n_Integer] := (NCQ[m]; NCS[n]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := NCmonos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); NCQS[m_Integer,q_,n_Integer] := (NCQ[m,q]; NCS[n]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := NCmonos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); NCQS[m_Integer,n_Integer,r_] := (NCQ[m]; NCS[n,r]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := NCmonos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); NCQS[m_Integer,q_,n_Integer,r_] := (NCQ[m,q]; NCS[n,r]; QSw = Qw + Sw; QSnorm = {{},Qa,Sb,{},{}}; QSm[i_] := NCmonos[i+#,QSw + #]& /@ DepWeight; ncomp=0; SaveData;); (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) CQ[m_] := (Qw = m; Qm = monos[m]; Qcomp; ); CQ[m_,q_List] := (Qw = m; Qm = monos[m]; Qcompl[q];); CQ[m_,q_] := CQ[m,{q}]; (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) NCQ[m_] := (Qw = m; Qm = NCmonos[m]; Qcomp; ); NCQ[m_,q_List] := (Qw = m; Qm = NCmonos[m]; Qcompl[q]; ); NCQ[m_,q_] := NCQ[m,{q}]; (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) CS[m_] := (Sw = m; Sm = monos[m]; Scomp; ); CS[m_,q_List] := (Sw = m; Sm = monos[m]; Scompl[q]; ); CS[m_,q_] := CS[m,{q}]; (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522; endGroup] *) NCS[m_] := (Sw = m; Sm = NCmonos[m]; Scomp; ); NCS[m_,q_List] := (Sw = m; Sm = NCmonos[m]; Scompl[q]; ); NCS[m_,q_] := NCS[m,{q}]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; rightWrapOffset = 558; startGroup] Qcomp --- Determine form of evolution equation Qcompl[q] --- Determine form of evolution equation from list q Qcompq[q] --- Determine form of specific evolution equation q Scomp --- Determine form of symmetry Scompl[q] --- Determine form of symmetry from list q Scompq[q] --- Determine form of specific symmetry q ;[s] 22:0,1;6,0;9,2;48,1;59,0;62,2;113,1;125,0;128,2;178,1;179,2;180,1;186,0;189,2;218,1;229,0;232,2;273,1;285,0;288,2;328,1;329,2;331,-1; 3:6,14,10,Palatino,3,12,0,0,0;8,12,10,Courier,0,12,0,0,0;8,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) Qcomp := (ClearAll[a]; Qarep = {}; Qmf = Qm//Flatten; Qc = Table[Table[constants[a,Length[Qm[[i,j]]],i,j], {j,Length[Qm[[i]]]}],{i,NDepVars}]; Qcf = Flatten[Qc]; Qa = Select[Qcf,NonNumberQ]; Ql = Table[(Thread[dotit[Qc[[i]],Qm[[i]]]]) /.dotit -> Dot,{i,NDepVars}]; Q = Apply[Plus, Ql,{1}]); (* :[font = input; initialization; preserveAspect] *) Qcompl[q_] := If[PosIntListQ[q], Qm = EmptyList[#,q]& /@ Qm; Qcomp, Qcompq[q]; If[q == Q,, printline["** Warning - some terms have the wrong weight! **"], printline["** Warning - some terms have the wrong weight! **"]]]; (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) Qcompq[q_] := Module[{cc},Qmf = Qm//Flatten; Qc = Thread[cc[q,Qm]] /. cc -> Qgetcfs; Qcf = Union[Level[Qc,{-1}]]; Qa = Select[Qcf,NonNumberQ]; const[Qa]; Ql = Table[(Thread[dotit[Qc[[i]],Qm[[i]]]]) /.dotit -> Dot,{i,NDepVars}]; Q = Apply[Plus, Ql,{1}]]; (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) Scomp := (ClearAll[b];Sbrep = {}; Smf = Sm//Flatten; Sc = Table[Table[constants[b,Length[Sm[[i,j]]],i,j], {j,Length[Sm[[i]]]}],{i,NDepVars}]; Scf = Flatten[Sc]; Sb = Select[Scf,NonNumberQ]; Sl = Table[(Thread[dotit[Sc[[i]],Sm[[i]]]]) /.dotit -> Dot,{i,NDepVars}]; S = Apply[Plus, Sl,{1}]); (* :[font = input; initialization; preserveAspect] *) Scompl[q_] := If[PosIntListQ[q], Sm = EmptyList[#,q]& /@ Sm; Scomp, Scompq[q]; If[q == S,, printline["** Warning - some terms have the wrong weight! **"], printline["** Warning - some terms have the wrong weight! **"]]]; (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522] *) Scompq[q_] := Module[{cc},Qmf = Qm//Flatten; Sc = Thread[cc[q,Sm]] /. cc -> Qgetcfs; Scf = Union[Level[Sc,{-1}]]; Sb = Select[Scf,NonNumberQ]; const[Sb]; Sl = Table[(Thread[dotit[Sc[[i]],Sm[[i]]]]) /.dotit -> Dot,{i,NDepVars}]; S = Apply[Plus, Sl,{1}]]; (* :[font = input; initialization; preserveAspect; rightWrapOffset = 522; endGroup] *) Qgetcfs[e_,l_] := Map[Coefficient[e,#]&,l,{2}]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] EmptyList[l,q] --- Replace entries of list l which are not listed in q by an empty list ;[s] 10:0,2;4,1;19,0;23,2;49,1;50,2;77,1;78,2;82,1;114,2;131,-1; 3:1,14,10,Palatino,3,12,0,0,0;4,12,10,Courier,0,12,0,0,0;5,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) EmptyList[l_,q_] := ReplacePart[l,{}, List /@ Complement[Range[Length[l]],q]] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Qpr --- Print current equation Spr --- Print current symmetry QSpr --- Print current equation and current symmetry Qprtx --- Print TeX form of current equation Sprtx --- Print TeX form of current symmetry QSprtx --- Print TeX form of current equation and current symmetry ;[s] 19:0,2;4,1;8,0;11,2;41,1;45,0;48,2;78,1;83,0;86,2;137,1;143,0;146,2;188,1;194,0;197,2;239,1;246,0;249,2;308,-1; 3:6,14,10,Palatino,3,12,0,0,0;6,12,10,Courier,0,12,0,0,0;7,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Qp[{e_,a_,b_,r_,s_}] := Do[printeq[DepVart[[i]],Collect[Q[[i]]/.r//ex,Qmf]], {i,NDepVars}]; Qpr := (printline["** Equation **"];Qp[QSnorm]); Sp[{e_,a_,b_,r_,s_}] := Do[printeq[DepVart[[i]],Collect[S[[i]]/.s//ex,Smf]], {i,NDepVars}]; Spr := (printline["** Symmetry **"];Sp[QSnorm]); QSp[l_] := (printline["** Equation **"];Qp[l]; printline["** Symmetry **"];Sp[l];); QSpr := (Qpr;Spr); Qptx[{e_,a_,b_,r_,s_}] := Do[ptxeq[DepVart[[i]], Collect[Q[[i]]/.r//ex,Qmf]],{i,NDepVars}]; Qprtx := (printline["** Equation **"];Qptx[QSnorm]); Sptx[{e_,a_,b_,r_,s_}] := Do[ptxeq[DepVart[[i]], Collect[S[[i]]/.s//ex,Smf]],{i,NDepVars}]; Sprtx := (printline["** Symmetry **"];Sptx[QSnorm]); QSptx[l_] := (printline["** Equation **"];Qptx[l]; printline["** Symmetry **"];Sptx[l];); QSprtx := (Qprtx;Sprtx); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] DepVart --- t derivatives of dependent variables ;[s] 6:0,2;4,1;12,0;15,2;20,1;22,2;57,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) DepVart := dd[t][#]& /@ DepVars; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] QSave[file] --- Save current data in file "file". ;[s] 6:0,2;4,1;16,0;19,2;50,1;56,2;58,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) QSave[f_] := Save[ToString[f],ncomp,niter, Qc,Qcf,Ql,Q,Qw,Qm,Qa, Sc,Scf,Sl,S,Sw,Sm,Sb, QSm,QS,QSw,QSnorm,QSbranch]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] SaveData --- Save data from current computation RestoreData[n] or Qres[n]--- Restore data from (end of) computation number n ( n = 0 for original data ) RestoreData[n,m] or Qres[n,m]--- Restore data from (end of) computation number n ( n = 0 for original data ) and go to branch number n ;[s] 29:0,2;2,1;11,0;14,2;54,1;69,2;71,1;80,0;83,2;153,1;154,2;160,1;161,2;163,1;171,2;190,1;210,2;212,1;223,0;226,2;296,1;297,2;303,1;304,2;306,1;314,2;333,1;343,2;369,1;399,-1; 3:3,14,10,Palatino,3,12,0,0,0;13,12,10,Courier,0,12,0,0,0;13,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) SaveData := ( printline["** Saving data from Computation #",ncomp]; data[ncomp] = {ncomp,niter, Qc,Qcf,Ql,Q,Qw,Qm,Qa, Sc,Scf,Sl,S,Sw,Sm,Sb, QSm,QS,QSw,QSnorm,QSbranch};) RestoreData[n_] := ({ncomp,niter, Qc,Qcf,Ql,Q,Qw,Qm,Qa, Sc,Scf,Sl,S,Sw,Sm,Sb, QSm,QS,QSw,QSnorm,QSbranch} = data[n];); RestoreData[n_,m_] := (RestoreData[n]; Setbr[m];); Qres = RestoreData (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] QSpreq --- Print remaining nonlinear equations ;[s] 4:0,2;3,1;10,0;13,2;52,-1; 3:1,14,10,Palatino,3,12,0,0,0;1,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) QSpreq := If[Length[QSnorm[[1]]] == 0, printline["There are no unsolved nonlinear equations."], printline["Unsolved Nonlinear Equations"]; printnum[QSnorm[[1]]]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] QSadd[c1,n1,c2,n2,...] --- Add additional normalizations c1->n1, c2->n2, ... for equation & symmetry coefficients. QSadda[{c1,c2,...},a] --- Add normalizations c1->a, c2->a, ... QSadd0[c1,c2,...] --- Add normalizations c1->0, c2->0, ... QSadd1[c1,c2,...] --- Add normalizations c1->1, c2->1, ... QSadd2[c1,c2,...] --- Add normalizations c1->2, c2->2, ... QSaddm1[c1,c2,...] --- Add normalizations c1->-1, c2->-1, ... QSaddm2[c1,c2,...] --- Add normalizations c1->-2, c2->-2, ... Normalize[x -> c,QSnorm] --- Add normalization x -> c to QSnorm ;[s] 32:0,2;2,1;25,0;28,2;89,1;109,2;150,1;172,0;175,2;200,1;218,2;220,1;238,0;241,2;266,1;303,0;306,2;331,1;368,0;371,2;396,1;434,0;437,2;462,1;502,0;505,2;530,1;576,0;579,2;602,1;611,2;613,1;623,-1; 3:8,14,10,Palatino,3,12,0,0,0;12,12,10,Courier,0,12,0,0,0;12,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) QSadd[x_,a_] := (QSnorm = NormalizeT[x,a,QSnorm];); QSadd[x_,a_,b__] := (QSadd[x,a]; QSadd[b];); QSadd[x_ -> a_] := (QSnorm = NormalizeT[x,a,QSnorm];); QSadd[x_ -> a_,b__] := (QSadd[x,a]; QSadd[b];); QSadd[l_List] := QSadd @@ l; QSadda[l_,a_] := Do[QSnorm = NormalizeT[l[[i]],a,QSnorm], {i,Length[l]}]; QSadd0[x__] := QSadda[{x},0]; QSadd1[x__] := QSadda[{x},1]; QSadd2[x__] := QSadda[{x},2]; QSaddm1[x__] := QSadda[{x},-1]; QSaddm2[x__] := QSadda[{x},-2]; Normalize[x_ -> c_,{e_,a_,b_,r_,s_}] := Module[{m,n,t}, t = (x -> c); m = Position[a,x]; n = Position[b,x]; If[Length[m] == 0, If[Length[n] == 0, If[NWarn, Print["Warning: variable ",x," cannot be normalized"]]; {e,a,b,r,s}, {Select[e/.t//ex,NonZeroQ],a,DeleteCases[b,x], r/.t//ex,Append[s/.t//ex,t]}], {Select[e/.t//ex,NonZeroQ],DeleteCases[a,x],b, Append[r/.t//ex,t],s/.t//ex}]]; Normalize[x_,c_,l_] := Normalize[x -> c,l]; NormalizeT[z__] := (NWarn = True; Normalize[z]); NormalizeF[z__] := (NWarn = False; Normalize[z]); NWarn = True; (* :[font = info; inactive; preserveAspect] Basic Solution routine ;[s] 1:0,1;31,-1; 2:0,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] QSolve[m] --- Solve commutator equations of order m using LSolve, first for coefficients of symmetry, then for coefficients of equation. QSolve[m,n] --- Solve commutator equations of order m for branch number n as determined by Branch ;[s] 15:0,1;10,0;13,2;54,1;73,2;81,1;87,2;198,1;210,0;213,2;254,1;273,2;293,1;295,2;314,1;321,-1; 3:2,14,10,Palatino,3,12,0,0,0;7,12,10,Courier,0,12,0,0,0;6,14,10,Palatino,2,12,0,0,0; :[font = input; inactive; preserveAspect] QSnorm = {e,a,b,r,s} --- current normalizations and equations to solve e --- Remaining nonlinear equations to be solved a --- Remaining unknown coefficients in equation b --- Remaining unknown coefficients in symmetry r --- Current normalizations for coefficients in equation s --- Current normalizations for coefficients in symmetry QSmon --- Current monomials in commutator QSeqs --- New equations arising from commutator QSjnorm = {f,a,b,r,s} --- new equations and normalizations f --- Current equations to be solved, including remaining unsolved nonlinear equations from previous step a, b, r, s --- Same as in QSnorm ncomp --- Current computation number (order) ;[s] 38:0,1;21,2;24,3;82,1;84,2;87,3;142,1;144,2;147,3;201,1;203,2;206,3;260,1;262,2;265,3;328,1;330,2;333,3;389,1;395,2;398,3;434,1;440,2;443,3;485,1;507,2;510,3;556,1;558,2;561,3;702,1;713,2;716,3;735,1;741,3;743,1;749,2;752,3;791,-1; 4:0,12,10,Courier,1,12,0,0,0;13,12,10,Courier,0,12,0,0,0;12,14,10,Palatino,3,12,0,0,0;13,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) QSolve[m_] := Module[{c,i},ncomp = m; QSmon = QSm[m]; If[Plus @@ Length /@ QSmon ==0, Print["No terms in this degree"], printline["** Starting Commutator Computation for Order ", ncomp]; QSeqc = Flatten[Thread[c[QS[m],QSmon]]/.c -> CoefList]; QSeqs = Select[QSeqc /. QSnorm[[5]] /. QSnorm[[4]] //ex, NonZeroQ]; i = Length[QSeqs]; If[i == 0, printline["No additional equations in this degree"], printline[i,"additional equations to be solved"]; QSjnorm = QSjoin[QSeqs,QSnorm]; QSnorm = LSolve[QSjnorm]; QSnorm[[1]] = Pruneq[QSnorm[[1]]]]; SaveData; i = Length[QSnorm[[1]]]; If[i>0,printline["** There are",i, "unsolved nonlinear equations **"]]];]; QSolve[m_,n_] := (Choosebr[n]; QSolve[m];); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] QSjoin[e,QSnorm] --- Join additional equations to QSnorm and eliminate redundancies. ;[s] 6:0,2;1,1;18,0;21,2;55,1;79,2;119,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) QSjoin[e_,{f_,a_,b_,r_,s_}] := {Union[e/.s/.r//ex,f],a,b,r,s} (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] LSolve[{e,a,r}] --- Solve system of equations e==0 for variables in a as much as possible. Only equations such that one of the variables in a has a constant coefficient are solved, r is current replacement list. The result is the list {f,b,s}, where f is the remaining unsolved nonlinear equations b is the remaining unnormalized variables r is the updated replacment list LSolve[{e,a,b,r,s}] --- Combination of LSolve[{e,b,s}] and LSolve[{e,a,r}] Lprintline[a] --- Prints line a when LSolvePrintFlag is True ;[s] 38:0,2;3,1;19,0;22,2;71,1;76,2;95,1;105,2;163,1;172,2;188,1;189,2;249,1;250,2;325,1;334,2;369,1;371,2;445,1;446,2;517,1;518,2;557,1;577,0;580,2;601,1;651,2;654,1;670,2;675,1;689,0;692,2;710,1;712,2;719,1;735,2;740,1;746,2;772,-1; 3:3,14,10,Palatino,3,12,0,0,0;17,12,10,Courier,0,12,0,0,0;18,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) LSolve[{e_,a_,r_}] := Module[{ee,aa,rr,ii=1,f,i,j,tf,l,r1,a1,c,con=False}, ee = e; aa = a; rr = r; Do[f = ee[[ii]]/.rr //ex; If[NumberQ[f],If[f==0,, Lprintline["** Contradiction: ",f," == 0"]; con=True,Break[]], tf = True; j = 1; l = Length[aa]; While[tf && (j <= l), a1 = aa[[j]]; If[FreeQ[f,a1], j = j+1, c = CoefficientList[f,a1]; If[Length[c] == 2 && NumberQ[c[[2]]], ee = Delete[ee,ii]; ii = ii-1; aa = DeleteCases[aa,a1]; r1 = a1 -> (- c[[1]]/c[[2]])//ex; Lprintline[r1]; rr = rr /. r1 //ex; rr = Append[rr,r1]; tf = False, j=j+1]]]]; ii=ii+1,{i,Length[ee]}]; If[con,{{1},aa,rr}, ee = Select[ee/.rr //ex,NonZeroQ]; {ee,aa,rr}]]; (* :[font = input; initialization; preserveAspect] *) LSolve[{e_,a_,b_,r_,s_}] := Module[{ee,ef,aa,bb,rr,ss}, {ee,bb,ss} = LSolve[{e,b,s}]; rr = r/.ss//ex; {ef,aa,rr} = LSolve[{ee,a,rr}]; ss = ss/.rr//ex; {ef,aa,bb,rr,ss}] (* :[font = input; initialization; preserveAspect] *) Lprintline[a__] := If[LSolvePrintFlag,printline[a]]; (* :[font = input; initialization; preserveAspect; endGroup] *) LSolvePrintFlag = True; (* :[font = message; inactive; preserveAspect; startGroup] Qh[k] Ñ terms in equation of degree k Qv[k] Ñ evolutionary vector field from terms in equation of degree k Sh[k] Ñ terms in symmetry of degree k Sv[k] Ñ evolutionary vector field from terms in symmetry of degree k QS[m] --- terms in commutator of degree m QSc[m] --- coefficients of terms in commutator of degree m Zsel[l,k] --- returns either l[[k]] or 0 if k > Length[l] ;[s] 27:0,1;6,2;36,1;47,2;108,1;119,2;149,1;160,3;161,2;221,1;232,3;235,2;257,1;258,2;267,1;279,3;282,2;320,1;321,2;330,1;345,3;348,2;365,1;374,2;376,1;381,2;383,1;401,-1; 4:0,12,10,Courier,0,12,65535,0,0;12,12,10,Courier,0,12,0,0,0;11,14,10,Palatino,2,12,0,0,0;4,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect] *) Qh[k_] := Zsel[#,k]& /@ Ql; Sh[k_] := Zsel[#,k]& /@ Sl; Qv[k_] := vq[Qh[k]]; Sv[k_] := vq[Sh[k]]; (* :[font = input; initialization; preserveAspect] *) QS[m_] := (Qc = Qc // ex; Sc = Sc //ex; Sum[Qv[i] /@ Sh[m-i] - Sv[i] /@ Qh[m-i], {i,0,m}] //ex) (* :[font = input; initialization; preserveAspect] *) QSc[m_] := Module[{c}, Flatten[Thread[c[QS[m],QSm[m]]]/. c -> CoefList]]; (* :[font = input; initialization; preserveAspect; endGroup] *) Zsel[l_List,k_] := If[k>=Length[l],0,l[[k+1]]] (* :[font = message; inactive; preserveAspect; fontColorRed = 0; startGroup] leadcoef[x] Ñ leading coefficient of expression x lead1[x] = x /leadcoef[x] Ñ so leading coefficient is now 1 Pruneq[e] Ñ eliminate redundate equations in list e == 0 ;[s] 7:0,0;12,1;49,0;80,1;112,0;127,1;168,0;178,-1; 2:4,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) leadcoef[x_Times] := If[NumberQ[x[[1]]],x[[1]],1]; leadcoef[x_Plus] := leadcoef[x[[1]]]; leadcoef[x_] := If[NumberQ[x],x,1]; lead1[x_] := x/leadcoef[x]//ex; (* :[font = input; initialization; preserveAspect; endGroup] *) Pruneq[e_] := Union[lead1 /@ Select[e,NonZeroQ]]; (* :[font = info; inactive; preserveAspect] Branching routines ;[s] 2:0,0;1,1;24,-1; 2:1,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] QFactor[l] --- Factor terms in list l QLFactor[{e,a}] --- Factor equations e, retaining replacements in a Factorz[e] --- Factor z over complex numbers ;[s] 18:0,2;1,1;12,0;15,2;41,1;43,2;44,1;60,0;63,2;85,1;88,2;119,1;121,2;122,1;133,0;136,2;149,1;151,2;172,-1; 3:3,14,10,Palatino,3,12,0,0,0;7,12,10,Courier,0,12,0,0,0;8,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) QFactor[l_] := Union[ Union /@ Distribute[Factorz /@ l,List]]; (* :[font = input; initialization; preserveAspect] *) QLFactor[{e_,a__}] := {#,a}& /@ QFactor[e] (* :[font = input; initialization; preserveAspect; endGroup] *) Factorz[f_] := Select[First /@ FactorList[f,GaussianIntegers -> True],NonNumberQ] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] uneqrepl[a != 0] --- normalizes nonzero variable to a -> 1 neqsolve[r] --- normalizes nonzero variables in list r to have value 1 and then solves remaining equations ;[s] 13:0,2;1,1;18,0;21,2;57,1;64,2;65,1;77,0;80,2;122,1;124,2;138,1;158,2;194,-1; 3:2,14,10,Palatino,3,12,0,0,0;5,12,10,Courier,0,12,0,0,0;6,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) uneqrepl[a_ != 0] := a -> 1; neqsolv[r_] := Module[{cs = Cases[r,a_ != b_],rs,ss,zs}, If[cs == {},rs = r; zs = {}, zs = uneqrepl /@ cs;rs = r/.zs]; ss = Solve[rs]; If[Length[ss] > 0,{ss[[1]],zs},{ss,zs}]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] RSolve[f] --- Adapts Mathematica's Reduce routine to solve system of nonlinear equations f == 0 Note: nonzero variables are normalized to have value 1 ;[s] 9:0,2;1,1;11,0;14,2;41,1;48,2;133,1;140,2;228,1;231,-1; 3:1,14,10,Palatino,3,12,0,0,0;4,12,10,Courier,0,12,0,0,0;4,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) RSolve[f_] := Module[{s},s = Reduce[f==0]; If[Head[s]===Or, If[RSflag,printlist[s]]; av = Union[neqsolv /@ List @@ s], If[RSflag,printline[s]]; av = Union[{neqsolv[s]}]]; If[RSflag,printnum[av]]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Setbrs --- Determine branches by trying to solve current system of nonlinear equations Selectbr[a] --- Determine if a is a nontrivial branch Addbr[{e,r}] --- Complete normalizations for new branch ;[s] 12:0,2;1,1;8,0;11,2;126,1;138,0;141,2;161,1;164,2;188,1;201,0;204,2;246,-1; 3:3,14,10,Palatino,3,12,0,0,0;4,12,10,Courier,0,12,0,0,0;5,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) Selectbr[a_] := Module[{fa = Flatten[a]}, {Select[ae /. fa//ex,NonZeroQ], Sort[Join[ar /. fa//ex,fa]], Join[an,a[[2]]]}]; (* :[font = input; initialization; preserveAspect; endGroup] *) Setbrs:= (RSolve[ae[[1]]]; If[av[[1,1]] === {}, If[Cflag,printline["** Contradiction **"]; printline[ae]; printline[ar]; printline[an]], brs = Join[brs,Selectbr /@ av]]; If[BuPrintFlag, printline["** There are now ",Length[brs], "unsolved branches **"]]); (* :[font = input; initialization; preserveAspect] *) Addbr[{e_,r_}] := Module[{i,q=QSnorm}, Do[q = NormalizeF[r[[i]],q],{i,Length[r]}]; QSjoin[e,q]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] StrictsubsetBQ DeleteSupersetBQ --- Used to eliminate brancing redundancies ;[s] 6:0,2;1,1;16,2;17,1;34,0;37,2;80,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) StrictSubsetBQ[{l_,x_},{m_,y_}] := And @@ Append[MemberQ[m,#]& /@ l, Sort[l]=!=Sort[m]]; DeleteSupersetsB[l_List] := Module[{s,x,m}, m = Union[Sort /@ l]; s[x_] := ! Or @@ (StrictSubsetBQ[#,x]& /@ m); Select[m,s]]; SuperAppendB[l_,x_] := DeleteSupersetsB[Append[l,x]]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Branch --- Determine branches by trying to solve current system of nonlinear equations Note: nonzero coefficients are normalized to have value 1 ;[s] 5:0,2;1,1;8,0;11,2;215,1;218,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Branch := (Off[Solve::svars]; ae = QSnorm[[1]]; ar = {}; an = {}; brs = {}; bran = {}; Setbrs; l=0; While[Length[brs] > 0, l = l + 1; If[BuPrintFlag, printline["Branch #",l]]; {ae,ar,an} = Last[brs]; brs = Drop[brs,-1]; If[Length[ae]==0, bran = SuperAppendB[bran,{ar,an}]; lbran = Length[bran]; If[BPrintFlag, printline["** New branch: now",lbran, "solved branches **"];Qpar], Setbrs]]; printline["** There are",Length[bran],"branches **"]; QSbranch = Addbr /@ bran; SaveData;); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Choosebr[n] --- Choose branch n for further computations ;[s] 6:0,2;1,1;13,0;16,2;33,1;37,2;63,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Choosebr[n_] := (If[n > 0 && n <= Length[QSbranch], QSnorm = QSbranch[[n]], printline["*** Illegal Branch ****"], printline["*** Illegal Branch ****"]];); (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Default branching print flags ;[s] 1:0,1;34,-1; 2:0,14,10,Palatino,3,12,0,0,0;1,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) BPrintFlag = False; BuPrintFlag = False; RSflag =False; LSolvePrintFlag = False; Cflag = True; (* :[font = input; inactive; wordwrap; dontPreserveAspect] Qpbrs --- Print equation branches Qpbr[k] --- Print equation branch number k Qpar --- Print the equation branch under consideration Spbrs --- Print symmetry branches Spbr[k] --- Print symmetry branch number k QSpbrs --- Print equation and symmetry branches QSpbr[k] --- Print equation and symmetry branch number k Qpbrstx --- Print TeX form of equation branches Qpbrtx[k] --- Print TeX form of equation branch number k Spbrstx --- Print TeX form of symmetry branches Spbrtx[k] --- Print TeX form of symmetry branch number k QSpbrstx --- Print TeX form of equation and symmetry branches QSpbrtx[k] --- Print TeX form of equation and symmetry branch number k ;[s] 51:0,2;4,1;10,0;13,2;44,1;52,0;55,2;89,1;90,2;95,1;100,0;103,2;156,1;162,0;165,2;196,1;204,0;207,2;241,1;242,2;247,1;254,0;257,2;301,1;310,0;313,2;360,1;361,2;366,1;374,0;377,2;420,1;430,0;433,2;479,1;480,2;485,1;493,0;496,2;539,1;549,0;552,2;598,1;599,2;604,1;613,0;616,2;672,1;683,0;686,2;745,1;747,-1; 3:13,14,10,Palatino,3,12,0,0,0;19,12,10,Courier,0,12,0,0,0;19,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) Qpbr[k_] := (printline["** Equation #",k,"**"]; Qp[QSbranch[[k]]];); Qpbrs := (Do[Qpbr[k],{k,Length[QSbranch]}];); Qpar := (printline["** Equation **"]; printeq[DepVart,Q /. ar];); Spbr[k_] := (printline["** Symmetry #",k,"**"]; Sp[QSbranch[[k]]];); Spbrs := (Do[Spbr[k],{k,Length[QSbranch]}];); QSpbr[k_] := (printline["** Equation #",k,"**"]; Qp[QSbranch[[k]]]; printline["** Symmetry #",k,"**"]; Sp[QSbranch[[k]]];); QSpbrs := (Do[QSpbr[k],{k,Length[QSbranch]}];); Qpbrtx[k_] := (printline["** Equation #",k,"**"]; Qptx[QSbranch[[k]]];); Qpbrtxs := (Do[Qpbrtx[k],{k,Length[QSbranch]}];); Spbrtx[k_] := (printline["** Symmetry #",k,"**"]; Sptx[QSbranch[[k]]];); Spbrstx := (Do[Spbrtx[k],{k,Length[QSbranch]}];); QSpbrtx[k_] := (printline["** Equation #",k,"**"]; Qptx[QSbranch[[k]]]; printline["** Symmetry #",k,"**"]; Sptx[QSbranch[[k]]];); QSpbrstx := (Do[QSpbrtx[k],{k,Length[QSbranch]}];); (* :[font = info; inactive; preserveAspect] ********************************************************************** Printing Routines ********************************************************************** ;[s] 3:0,0;72,1;94,0;165,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] TeXFormat[f] --- TeX form of f TeXMatrix[f] --- TeX matrix form of f ;[s] 11:0,3;5,1;17,2;21,3;35,1;38,3;43,1;55,2;59,3;80,1;88,3;90,-1; 4:0,12,10,Courier,1,12,0,0,0;4,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;5,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) TeXFormat[dd[x__][f_]] := SequenceForm[f,"_{",x,"}"]; TeXFormat[pd[x__][f_]] := SequenceForm["{\partial ",TexFormat[f], "\over \partial ",TexFormat[x],"}"]; TeXFormat[f_+g_] := TeXFormat[f] + TeXFormat[g]; TeXFormat[f_/g_] := SequenceForm["{", TeXFormat[Numerator[f/g]],"\over ", TeXFormat[Denominator[f/g]],"}"]; TeXFormat[f_ g_] := TeXFormat[f] TeXFormat[g]; TeXFormat[f_Times^n_] := SequenceForm["(",TeXFormat[f],")^{", TeXFormat[n],"}"]; TeXFormat[f_Plus^n_] := SequenceForm["(",TeXFormat[f],")^{", TeXFormat[n],"}"]; TeXFormat[f_^n_] := SequenceForm[TeXFormat[f],"^{", TeXFormat[n],"}"]; TeXFormat[r_Rational] := SequenceForm["{",Numerator[r], "\over ",Denominator[r]"}"]; TeXFormat[l_List] := SequenceForm["\{", Infx[TeXFormat /@ l," , "],"\}"]; TeXFormat[f_[g__]] := TeXFormat /@ f[g]; TeXFormat[f_ -> g_] := SequenceForm[TeXFormat[f], " \longrightarrow ",TeXFormat[g]]; TeXFormat[Equal[f_,g_]] := SequenceForm[TeXFormat[f], " = ",TeXFormat[g]]; TeXFormat[f_] := f; TeXMatrix[l_] := SequenceForm["\pmatrix {", Infx[TeXMatrixRow /@ l," \cr "],"\cr }"]; TeXMatrixRow[l_] := Infx[TeXFormat /@ l," & "]; tx = TeXFormat; txm = TeXMatrix; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] blank --- Print a blank line prstars --- Print a line of stars printline[l] --- Print l followed by a blank line printlist[l] --- Print elements in list l line by line printnum[l] --- Print elements in list l with numbering line by line printeq[x,y] --- Print x = y (line by line for lists) printmatrix[l] --- Print matrix m followed by a line of stars ptxline[l] --- Print Tex form of l followed by a blank line ptxlist[l] --- Print Tex form of elements in list l line by line ptxnum[l] --- Print Tex form of numbered elements in list l ptxeq[x,y] --- Print Tex form of x = y (line by line for lists) ptxmatrix[l] --- PrintTex form of matrix m followed by stars ;[s] 66:0,3;5,1;10,2;14,3;40,1;48,2;52,3;81,1;93,2;97,3;107,1;108,4;109,3;140,1;152,2;156,3;183,1;184,4;185,3;204,1;215,2;219,3;246,1;247,4;248,3;282,1;294,2;298,3;308,1;315,4;316,3;348,1;362,2;366,3;383,1;384,4;385,3;419,1;429,2;433,3;455,1;456,4;457,3;488,1;498,2;502,3;541,1;542,4;543,3;562,1;571,2;575,3;623,1;624,3;630,1;640,2;644,3;666,1;673,4;674,3;706,1;718,2;722,3;750,1;751,4;752,3;782,-1; 5:0,12,10,Courier,1,12,0,0,0;22,12,10,Courier,0,12,0,0,0;12,14,10,Palatino,1,12,0,0,0;23,14,10,Palatino,2,12,0,0,0;9,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) blank := Print[" "]; stars = "*********************************************************"; prstars := (blank; Print[stars]; blank;) printline[x_] := (Print[x];blank); printline[x__] := (Print[Infix[{x}," "]];blank;) printlist[l_] := Scan[printline,l]; printnum[l_] := (Do[Print["#",i," ",l[[i]]];blank, {i,Length[l]}];); printeq[x_,y_] := (Print[x," = ",y];blank;); printeq[l_List,m_List] := Module[{p}, Thread[p[l,m]] /. p -> printeq;]; printmatrix[m_] := (Print[MatrixForm[m]]; prstars;); ptxline[x_] := (Print[x//TeXFormat];blank); ptxline[x__] := (Print[Infix[{x//TeXFormat}," "]];blank;) ptxnum[l_] := (Do[Print["#",i," ",l[[i]]//TeXFormat];blank, {i,Length[l]}];); ptxlist[l_] := Scan[ptxline,l]; ptxeq[x_,y_] := (Print[x//TeXFormat," = ",y//TeXFormat]; blank;); ptxeq[l_List,m_List] := Module[{p}, Thread[p[l,m]] /. p -> ptxeq;]; ptxmatrix[m_] := (Print[TeXMatrix[m]]; prstars;); pl = printline; pm = printmatrix; pn = printnum; po = printlist; ptl = ptxline; ptm = ptxmatrix; ptn = ptxnum; pto = ptxlist; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Infx[l,h] --- Modifies function Infix to allow lists l of length 1 Infc[l] --- Use comma for infix operator ;[s] 14:0,3;5,1;14,2;18,3;39,1;44,4;45,3;63,1;65,3;77,1;78,3;84,1;91,2;95,3;127,-1; 5:0,12,10,Courier,1,12,0,0,0;5,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;6,14,10,Palatino,2,12,0,0,0;1,14,10,Palatino,3,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Infx[l_,h_] := If[Length[l]==1,Prefix[l,""], Infix[l,h],Infix[l,h]]; Infc[l_] := Infx[l,","] (* :[font = info; inactive; preserveAspect] ********************************************************************** Output session to a file ********************************************************************** ;[s] 3:0,0;72,1;101,0;172,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] startout[f] --- Begin to output of session to file f stopout --- Stop output to current output file OutputFile --- Current output file ;[s] 14:0,3;5,1;16,2;20,3;57,1;69,2;73,3;109,1;111,3;113,1;126,2;130,3;151,1;153,3;155,-1; 4:0,12,10,Courier,1,12,0,0,0;5,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,1,12,0,0,0;6,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) startout[f_] := (If[OutputFile =!= Null, stopoutput]; OutputFile = ToString[f]; o = OpenWrite[OutputFile,FormatType -> OutputForm]; $Echo = Append[$Echo,o]; $Output = Append[$Output,o];); stopout := ( Close[OutputFile]; $Echo = DeleteCases[$Echo, OutputStream[OutputFile,_]]; $Output = DeleteCases[$Output, OutputStream[OutputFile,_]];); OutputFile = Null; (* :[font = info; inactive; preserveAspect] ********************************************************************** Initialization Routines ********************************************************************** ;[s] 3:0,0;72,1;100,0;171,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] DefaultIndepVars --- List of default independent variables DefaultDepVars --- List of default dependent variables ;[s] 7:0,3;4,1;21,0;24,3;69,1;84,2;88,3;127,-1; 4:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;1,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) DefaultIndepVars = {x,t,y,z,X,Y,Z}; DefaultDepVars = {u,v,w,U,V,W}; (* :[font = input; inactive; wordwrap; dontPreserveAspect] DefaultVfiCoefs, DefaultVfdCoefs --- List of default vector field coefficients DefaultVfChars --- List of default vector field characteristics ;[s] 7:0,3;4,1;37,2;40,3;140,1;155,2;158,3;205,-1; 4:0,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,1,12,0,0,0;3,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Clear[a,b,c,e,ao,bo,co,f,g,h,fo,go,ho,q,r,s,qo,ro,so]; DefaultVfiCoefs = {a,b,c,e,ao,bo,co}; DefaultVfdCoefs = {f,g,h,fo,go,ho}; DefaultVfChars = {q,r,s,qo,ro,so}; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Vset --- Define default vector field of order 0 Vsetup[i,f,c] --- Define vector field of order 0 with i as independent variable coeffcients f as dependent variable coeffcients c as characteristics ;[s] 17:0,2;4,1;9,0;12,2;53,1;55,2;59,1;73,0;76,2;109,1;110,2;119,1;142,2;178,1;200,2;234,1;256,2;275,-1; 3:2,14,10,Palatino,3,12,0,0,0;7,12,10,Courier,0,12,0,0,0;8,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect] *) Vset := Module[{i,j}, Clear[a,b,c,e,ao,bo,co,f,g,h,fo,go,ho,q,r,s,qo,ro,so]; i = Length[IndepVars]; j = Length[DepVars]; VfiCoefs = Take[DefaultVfiCoefs,i]; VfdCoefs = Take[DefaultVfdCoefs,j]; VfCoefs = Join[VfiCoefs,VfdCoefs]; VfChars = Take[DefaultVfChars,j]; VfVars = VxVars = Unx[0]; VfOrder = 0; VfString = vfstr[VfCoefs,VxVars]; depn[VfCoefs,0]; Release[VfChars] = char[VfiCoefs,VfdCoefs,IndepVars,DepVars]; Vsym := vp[VfiCoefs,VfdCoefs,IndepVars,DepVars] ]; (* :[font = input; initialization; preserveAspect; endGroup] *) Vsetup[i_,f_,c_] := ( VfiCoefs = i; VfdCoefs = f; VfCoefs = Join[VfiCoefs,VfdCoefs]; VfChars = c; VfVars = VxVars = Unx[0]; VfOrder = 0; VfString = vfstr[VfCoefs,VxVars]; depn[VfCoefs,0]; Release[VfChars] = char[VfiCoefs,VfdCoefs,IndepVars,DepVars]; V := vp[VfiCoefs,VfdCoefs,IndepVars,DepVars]) (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Evset[n] --- Define default evolutionary vector field of order n ;[s] 5:0,2;4,1;13,0;16,2;70,1;72,-1; 3:1,14,10,Palatino,3,12,0,0,0;2,12,10,Courier,0,12,0,0,0;2,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) Evset[n_] := Module[{j}, Clear[q,r,s]; j = Length[DepVars]; VfChars = VfCoefs = Take[DefaultVfChars,j]; VfVars = Unx[n]; VxVars = DepVars; VfOrder = n; VfString = vfstr[VfCoefs,VxVars]; depn[VfCoefs,n]; Vsym := vp[0& /@ IndepVars,VfChars,IndepVars,DepVars]] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] quickjet[i,j,k] --- Define default jet variables and vector fields for i independent variables j dependent variables and kth order of derivatives defined quickjet[i,j] --- same as quickjet[i,j,DefaultMaxDeriv] quickjet[i] --- same as quickjet[i,1,DefaultMaxDeriv] ;[s] 22:0,2;15,0;16,3;19,1;87,2;91,1;116,2;118,1;157,2;158,1;191,2;204,0;205,3;208,1;221,2;250,1;251,2;262,0;263,3;266,1;279,2;308,1;310,-1; 4:3,12,10,Courier,1,12,0,0,0;8,14,10,Palatino,2,12,0,0,0;8,12,10,Courier,0,12,0,0,0;3,14,10,Palatino,1,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) quickjet[i_Integer,j_Integer,k_Integer] := ( Clear[x,t,y,z,u,v,w]; iv = Take[DefaultIndepVars,i]; dv = Take[DefaultDepVars,j]; jetset[iv,dv,k]; Vset; ) /; (0 <= i <= Length[DefaultIndepVars] && 0 <= j <= Length[DefaultDepVars] && k >=0 ); quickjet[i_Integer,j_Integer] := quickjet[i,j,DefaultMaxDeriv] quickjet[i_Integer] := quickjet[i,1,DefaultMaxDeriv] (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] weightjet[{w},m,k] --- Define default jet variables and weights for 1 independent variable j dependent variables where {w} = {w1,...,wj} are the weights of dependent variables m is the weight of independent variable k is the order of derivatives that can be easily entered weightjet[{w},m] --- same as weightjet[{w},m,DefaultMaxDeriv] weightjet[{w}] --- same as weightjet[{w},1,DefaultMaxDeriv] weightjet[n,m,k] --- same as weightjet[{n},m,k] weightjet[n,m] --- same as weightjet[{n},m,DefaultMaxDeriv] weightjet[n] --- same as weightjet[{n},1,DefaultMaxDeriv] ;[s] 48:0,2;18,0;19,3;22,1;84,2;88,1;125,2;127,1;168,2;176,4;177,2;183,4;184,2;185,1;240,2;241,1;295,2;296,1;355,2;371,0;372,3;375,1;388,2;420,1;421,2;435,0;436,3;439,1;452,2;484,1;485,2;501,0;502,3;505,1;518,2;536,1;537,2;551,0;552,3;555,1;568,2;600,1;601,2;613,0;614,3;617,1;630,2;662,1;664,-1; 5:6,12,10,Courier,1,12,0,0,0;16,14,10,Palatino,2,12,0,0,0;18,12,10,Courier,0,12,0,0,0;6,14,10,Palatino,1,12,0,0,0;2,20,14,Courier,64,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) weightjet[w_List,m_Integer,k_Integer] := ( j = Length[w]; If[k >=0 && j <= Length[DefaultDepVars], quickjet[1,j,k]; DepWeight = w; minDepWeight = Min[w]; IndepWeight = m, printline["** Error in Parameters **"], printline["** Error in Parameters **"]];); weightjet[w_List,m_Integer] := weightjet[w,m,DefaultMaxDeriv]; weightjet[w_List] := weightjet[w,1,DefaultMaxDeriv]; weightjet[n_Integer,m_Integer,k_Integer] := weightjet[{n},m,k]; weightjet[n_Integer,m_Integer] := weightjet[{n},m,DefaultMaxDeriv]; weightjet[n_Integer] := weightjet[{n},1,DefaultMaxDeriv]; (* :[font = message; inactive; preserveAspect] ********************************************************************** Set Default Initialization ********************************************************************** ;[s] 3:0,0;72,1;103,0;174,-1; 2:2,12,10,Courier,0,12,0,0,65535;1,17,13,Courier,0,17,65535,0,0; :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Number of default jet variables: NDefaultIndepVars = 2 independent variables x,t NDefaultDepVars = 1 dependent variable u DefaultMaxDeriv = 3 order of derivatives which can be easily entered. (Does not limit order of derivatives in computations.) ;[s] 11:0,1;41,2;65,1;90,2;94,1;102,2;124,1;152,2;154,1;162,2;184,1;348,-1; 3:0,12,10,Courier,1,12,0,0,0;6,14,10,Palatino,2,12,0,0,0;5,12,10,Courier,0,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) NDefaultIndepVars = 2; NDefaultDepVars = 1; DefaultMaxDeriv = 3; (* :[font = input; initialization; preserveAspect] *) quickjet[NDefaultIndepVars,NDefaultDepVars,DefaultMaxDeriv]; (* :[font = input; inactive; wordwrap; dontPreserveAspect; startGroup] Default value for printing intermediate p.d.e.s when solving them. ;[s] 1:0,1;68,-1; 2:0,12,10,Courier,1,12,0,0,0;1,14,10,Palatino,2,12,0,0,0; :[font = input; initialization; preserveAspect; endGroup] *) PrintPDEQ = False; (* ^*)