(* ******************************************************************** :Nombre: YaLie.m :Version: 2.0 :Título: Simetrías de ecuaciones diferenciales :Autor: JM Díaz Moreno Departamento de Matemáticas Universidad de Cádiz España josemanuel.diaz@uca.es :Palabras claves: Lie; Grupos de Lie; Álgebras de Lie; Simetrías; Ecuaciones diferenciales, Ecuaciones diferenciales ordinarias; Ecuaciones en derivadas parciales; :Referencias: P. Olver, Applications of Lie Groups to Differential Equations Springer, New York, 1986 G.W. Bluman and S. Kumei, Symmetries and Differential Equations Springer, New York, 1989 ******************************************************************** *) (* --- *) (* *) (* Variables globales *) (* *) (* Xi = Infinitesimales Xi formales *) (* Phi = Infinitesimales Phi formales *) (* Pivotes = Factores cancelados *) (* Verborrea = Opción de información *) (* *) (* Variables globales privadas *) (* *) (* EQs = Ecuaciones *) (* Vd = Variables dependientes *) (* Vdd = Variables dependientes con dependencia *) (* NVd = Número de variables dependientes *) (* *) (* Vi = Variables independientes *) (* NVi = Número de variables independientes *) (* *) (* RSust = Sustituciones Mod H(u) = 0 *) (* *) (* XiF = Lista de infinitesimales Xi formales *) (* PhiF = Lista de infinitesimales Phi formales *) (* *) (* ------------------------------------------------------------ *) BeginPackage["YaLie`"] Print[" "]; Print["====================================================="]; Print[" "]; Print[" YaLie "]; Print[" Version 2.0 (2000) "]; Print[" @ JM Díaz "]; Print[" Departamento de Matemáticas "]; Print[" Universidad de Cádiz "]; Print[" "]; Print["====================================================="]; Print[" "]; (* Funciones accesibles *) YaLie::usage = "YaLie[Eqs, VDs, VIs, VPs]" YaLieTeX::usage = "YaLieTeX[archivo, EDs]" Prolongación::usage = "Prolongación[V, VDs, VIs, k]" EDeterminantes::usage = "EDeterminantes[Pr, VDs]" Simplifica::usage = "Simplifica[EDeterminantes]" SInfinitesimales::usage = "SInfinitesimales[EDs, LXi, LPhi]" Pivotes::usage = "Pivotes es una lista con los factores cancelados" DReduce::usage = "ReduceD[F, Derivative[m][u][x]==expr]" DFormat::usage = "DFormat[True] o DFormat[False]" (* Símbolos formales *) Xi::usage="Xi" Phi::usage="Phi" Verborrea::usage= "" (* Mensajes *) LieError::Símbolos = "Xi y Phi no pueden estar asignadas" LieError::Argumentos = "Error en los argumentos" LieError::Solve = "Solve no ofrece los resultados esperados" LieError::NoPolinómica = "Probablemente existen ecuaciones determinantes no polinómica en las derivadas" Begin["`Private`"] (* ============================================================ *) (* Lie *) (* ============================================================ *) YaLie[Eqs_, VDs_, VIs_, VPs_, Opciones___Rule]:= YaLie[Eqs, VDs, VIs, VPs, {}, {}, Opciones] YaLie[Eqs_, VDs_, VIs_, VPs_, LXi_, LPhi_, Opciones___Rule]:= Module[{EDs, MaxDer, u, P, i, TInicial, Verbo}, TInicial = TimeUsed[]; Verbo = Verborrea /. {Opciones} /. Options[Lie]; (* Comprobación de los argumentos y variables globales privadas *) If[ValueQ[Xi], Message[LieError::Símbolos]; Return[{}]]; If[ValueQ[Phi], Message[LieError::Símbolos]; Return[{}]]; If[!MatchQ[VDs, {(__Symbol)}], Message[LieError::Argumentos]; Return[{}]]; If[!MatchQ[VIs, {(__Symbol)}], Message[LieError::Argumentos]; Return[{}]]; Vd = VDs; Vdd = Map[Apply[#, VIs]&, VDs]; NVd = Length[Vd]; Vi = VIs; NVi = Length[Vi]; If[!VectorQ[Eqs, (Head[#]===Equal)&], Message[LieError::Argumentos]; Return[{}]]; EQs = Map[(Part[#,1]-Part[#,2])&, Eqs]; RSust = Solve[Map[(#==0)&, EQs], VPs]; If[Length[RSust]=!=1, Message[LieError::Solve]; Return[{}]]; RSust = Part[RSust, 1]; (* Proceso principal *) MaxDer = Cases[EQs, Derivative[m___][u_ /; MemberQ[Vd,u]][___]:>Plus[m], -1]; If[MaxDer==={}, MaxDer=0, MaxDer=Max[MaxDer]]; XiF = Map[Apply[#, Join[Vi, Vdd]]&, Table[Xi[j], {j,NVi}]]; PhiF = Map[Apply[#, Join[Vi, Vdd]]&, Table[Phi[j], {j,NVd}]]; (* Prolongación, sustitución y ecuaciones determinantes *) EDs = {}; For[i=1, i<=Length[EQs], i++, If[Verbo === True, Print[":: Ecuación (",i,") (",TimeUsed[]-TInicial,")"]]; P = ProlongacionE[Part[EQs, i], MaxDer]; If[Verbo === True, Print[" Prolongación (",TimeUsed[]-TInicial,")"]]; P = P /. RSust; P = EDeterminantesE[P]; If[Verbo === True, Print[" Ecuaciones determinantes (",TimeUsed[]-TInicial,")"]]; If[P === False, Message[LieError::NoPolinómica]; Return[{}]]; EDs = Join[EDs, P]]; (* Simplificación *) EDs = Simplifica[EDs]; If[Verbo === True, Print[":: Simplificación (",TimeUsed[]-TInicial,")"]]; (* Sustituir los infinitesimales si procede *) If[Length[LXi] > 0 && Length[LPhi] > 0, XiF = LXi; PhiF = LPhi; EDs = SustituyeInfinitesimales[EDs, Join[Vi, Vd], LXi, LPhi]; If[Verbo === True, Print[":: Sustitución de infintesimales (",TimeUsed[]-TInicial,")"]]]; Map[(#==0)&, EDs]]; Options[Lie] = {Verborrea->False} (* ============================================================ *) (* Página 1. Prolongaciones *) (* ============================================================ *) Unprotect[Prolongación, ProlongacionE, ProlongacionVP, Eta, DD]; Prolongación[Eq_, VDs_, VIs_, k_Integer?Positive]:= Module[{}, (* Comprobación de los argumentos *) If[ValueQ[Xi], Message[LieError::Símbolos]; Return[{}]]; If[ValueQ[Phi], Message[LieError::Símbolos]; Return[{}]]; If[!MatchQ[VDs, {(__Symbol)}], Message[LieError::Argumentos]; Return[{}]]; Vd = VDs; Vdd = Map[Apply[#, VIs]&, VDs]; NVd = Length[Vd]; If[!MatchQ[VIs, {(__Symbol)}], Message[LieError::Argumentos]; Return[{}]]; Vi = VIs; NVi = Length[Vi]; XiF = Map[Apply[#, Join[Vi, Vdd]]&, Table[Xi[j], {j,NVi}]]; PhiF = Map[Apply[#, Join[Vi, Vdd]]&, Table[Phi[j], {j,NVd}]]; ProlongacionE[Eq, k]] ProlongacionE[eqt_, k_Integer?Positive]:= Module[{Vk, pr, ff, i, xx, mm, eta, x}, (* k-ésima prolongación de v parcialmente desarrollado *) Vk = ProlongacionVP[k, eta]; (* La aplica a la ecuación *) pr = Through[Vk[eqt]]; pr = pr /. {Phi[i_] :> PhiF[[i]], Xi[i_] :> XiF[[i]]}; pr = pr /. eta[i_, x_] -> Eta[i,x]; (* Sólo en Phi y Xi porque después hay que sustituir ModH *) pr /. {Derivative[mm___][ff:(Xi|Phi)[_]][xx___] :> Derivative[mm][ff][Apply[Sequence, {xx} /. ToRules[Vdd==Vd]]], (ff:(Xi|Phi)[_])[xx___]:>ff[Apply[Sequence, {xx} /. ToRules[Vdd==Vd]]]} ] (* ------------------------------------------------------------ *) (* 1.1 pr^k(v) parcialmente desarrollada *) (* ------------------------------------------------------------ *) ProlongacionVP[k_Integer?Positive, eta_]:= Module[{DList, S, DT, Dd, m, i, F}, DList = {{}}; S = Sum[Xi[i] Dd[#, Vi[[i]]], {i,1,NVi}] + Sum[Phi[i] DT[#, Vdd[[i]]], {i,1,NVd}]; For[m=1, m<=k, m++, (* Construye todas las posibles parciales *) DList = Union[Map[Sort[Flatten[#]]&, Distribute[List[DList, Vi], List]]]; (* Suma sobre las variables dependientes *) For[i=1, i<=NVd, i++, (* Suma sobre J *) S += Sum[eta[i, DList[[j]]]* DT[#, D[Vdd[[i]], Apply[Sequence, DList[[j]]]]], {j, 1, Length[DList]}]] ]; S /. {F_ DT[#, x_]->(F D[#,x]&), F_ Dd[#,x_]->(F DD[#,x]&)}] (* ------------------------------------------------------------ *) (* 1.2 Eta_i^J *) (* ------------------------------------------------------------ *) Eta[i_, x_List]:= D[PhiF[[i]] - Sum[XiF[[j]] D[Vdd[[i]], Vi[[j]]], {j,1,NVi}], Apply[Sequence, x]] + Sum[XiF[[j]] D[Vdd[[i]], Apply[Sequence,x], Vi[[j]]], {j,1,NVi}] (* ------------------------------------------------------------ *) (* 1.3 Derivada para las variables independientes *) (* ------------------------------------------------------------ *) DD[ecuacion_, x_]:= Module[{eqt, DList, DLiss, u, m, xx}, (* Variables dependientes como constantes *) eqt = ecuacion /. ToRules[Vdd==Vd]; (* Derivadas de las variables dependiendes como constantes *) DList = Cases[eqt, Derivative[m___][u_ /; MemberQ[Vd,u]][___], -1]; DLiss = Cases[eqt, Derivative[m___][u_ /; MemberQ[Vd,u]][___]->u[m], -1]; eqt = eqt /. ToRules[DList==DLiss]; (* Derivada con respecto a x *) eqt = D[eqt, x]; (* Restaurar las variables dependientes y las derivadas *) eqt = eqt /. ToRules[Vd==Vdd]; DLiss = DLiss /. ToRules[Vd==Vdd]; eqt = eqt /. ToRules[DLiss==DList]]; Protect[Prolongación, ProlongacionE, ProlongacionVP, Eta, DD]; (* ============================================================ *) (* Página 2. Ecuaciones determinantes *) (* ============================================================ *) Unprotect[EDeterminantes, EDeterminantesE, Incógnitas]; EDeterminantes[Pr_, VDs_, VIs_]:= Module[{}, (* Comprobación de los argumentos *) If[ValueQ[Xi], Message[LieError::Símbolos]; Return[{}]]; If[ValueQ[Phi], Message[LieError::Símbolos]; Return[{}]]; If[!MatchQ[VDs, {(__Symbol)}], Message[LieError::Argumentos]; Return[{}]]; If[!MatchQ[VIs, {(__Symbol)}], Message[LieError::Argumentos]; Return[{}]]; Vd = VDs; Vdd = Map[Apply[#, VIs]&, VDs]; EDeterminantesE[Pr]] EDeterminantesE[Prolg_]:= Module[{u, DList, aux, i}, DList = Union[Cases[Prolg, Derivative[___][u_ /; MemberQ[Vd, u]][___], -1]]; Print[DList]; (* If[!PolynomialQ[Prolg, DList], Return[False]]; *) aux = Prolg; For[i=1, i<=Length[DList], i++, aux = Flatten[Map[CoefficientList[#, DList[[i]]]&, {aux}]]]; aux = Map[Collect[#, Incógnitas[#]]&, aux]; Union[DeleteCases[aux, 0]] /. ToRules[Vdd==Vd]] (* ------------------------------------------------------------ *) (* 2.1 Incógnitas *) (* ------------------------------------------------------------ *) Incógnitas[Eq_]:= Module[{}, Union[Cases[{Eq}, ((Xi|Phi)[_])[___], -1], Cases[{Eq}, Derivative[___][(Xi|Phi)[_]][___], -1]]]; Protect[EDeterminantes, EDeterminantesE, Incógnitas]; (* ============================================================ *) (* Página 3. Simplificación *) (* ============================================================ *) Unprotect[Simplifica, SimplificaE, FormaEstándarSimples, ReduceDSimples, Incógnitas, EsDerivadaQ]; Simplifica[EcuacionesD_]:= Module[{}, If[ValueQ[Xi], Message[LieError::Símbolos]; Return[{}]]; If[ValueQ[Phi], Message[LieError::Símbolos]; Return[{}]]; Pivotes = {}; SimplificaE[EcuacionesD]] SimplificaE[EcuacionesD_]:= Module[{M, P, F, S}, Pivotes = {}; M = Select[EcuacionesD, (Length[Incógnitas[#]]==1)&]; P = Select[EcuacionesD, (Length[Incógnitas[#]]>1)&]; S = {}; F = {}; While[M =!= {}, S = FormaEstándarSimples[S, M]; P = ReduceDSimples[S, P]; M = Select[P, (Length[Incógnitas[#]]==1)&]; P = Select[P, (Length[Incógnitas[#]]>1)&]]; Join[S, P]] (* ------------------------------------------------------------ *) (* 3.1 FormaEstándarSimples *) (* ------------------------------------------------------------ *) FormaEstándarSimples[S_, F_]:= Module[{SF, incgs, incgP, factores, k}, incgs = Incógnitas[F]; factores = F /. Map[(#->1)&, incgs]; factores = Select[factores, Not[NumericQ[#]]&]; Pivotes = Join[Pivotes, factores]; incgs = Union[S, incgs]; SF = {}; For[k=1, k<=Length[incgs], k++, If[!Apply[Or, Map[EsDerivadaQ[Part[incgs,k], #]&, Delete[incgs,k]]], SF = Append[SF, Part[incgs,k]]]]; SF] (* ------------------------------------------------------------ *) (* 3.2 ReduceDSimples *) (* ------------------------------------------------------------ *) ReduceDSimples[S_, F_]:= Module[{P, incgsF, incgP}, P = F /. Map[(#->0)&, S]; incgsF = Incógnitas[P]; While[incgsF =!= {}, incgP = First[incgsF]; incgsF = Rest[incgsF]; If[Apply[Or, Map[EsDerivadaQ[incgP, #]&, S]], P = P /. incgP->0]]; P] (* ------------------------------------------------------------ *) (* 3.3 Incógnitas *) (* ------------------------------------------------------------ *) Incógnitas[Eq_]:= Module[{}, Union[Cases[{Eq}, ((Xi|Phi)[_])[___], -1], Cases[{Eq}, Derivative[___][(Xi|Phi)[_]][___], -1]]]; (* ------------------------------------------------------------ *) (* 3.3 EsDerivadaQ *) (* ------------------------------------------------------------ *) EsDerivadaQ[Derivative[m___][V_[i_]][z___], V_[i_][z___]]:=True EsDerivadaQ[Derivative[m___][V_[i_]][z___],Derivative[m___][V_[i_]][z___]]:=False EsDerivadaQ[Derivative[m___][V_[i_]][z___], Derivative[n___][V_[i_]][z___]]:= Apply[And, NonNegative[{m}-{n}]] EsDerivadaQ[f_, g_]:=False Protect[Simplifica, SimplificaE, FormaEstándarSimples, ReduceDSimples, Incógnitas, EsDerivadaQ]; (* ============================================================ *) (* Página 4. Sustitución de infinitesimales *) (* ============================================================ *) Unprotect[SInfinitesimales, SustituyeInfinitesimales]; SInfinitesimales[EDets_, LXi_List, LPhi_List]:= Module[{Vars}, If[ValueQ[Xi], Message[LieError::Símbolos]; Return[{}]]; If[ValueQ[Phi], Message[LieError::Símbolos]; Return[{}]]; Vars = Union[Cases[EDets ((Xi|Phi)[_])[v___]->{v}, -1], Cases[EDets, Derivative[___][(Xi|Phi)[_]][v___]->{v}, -1]]; If[Vars==={}, Message[LieError::Argumentos]; Return[{}]]; SustituyeInfinitesimales[EDets, First[Vars], LXi, LPhi]]; SustituyeInfinitesimales[Eqsdets_, Variables_, LXi_, LPhi_]:= Module[{i, sustinf}, (* Muy sutil *) sustinf = Join[Table[Xi[i] -> Function[Evaluate[Variables], Evaluate[LXi[[i]]]], {i,NVi}], Table[Phi[i] -> Function[Evaluate[Variables], Evaluate[LPhi[[i]]]], {i,NVd}]]; Expand[Eqsdets /. sustinf]]; Protect[SInfinitesimales, SustituyeInfinitesimales]; (* ============================================================ *) (* Página 5. Función auxiliar: sustituir derivadas *) (* ============================================================ *) Unprotect[DReduce, DerivadaQ]; DReduce[F_, Derivative[m___][u_][x__]==expr_]:= Module[{dF, Dm, incgs, dif, subst}, dF = F /. Derivative[m][u][x] -> expr; incgs = Union[Cases[dF, Derivative[___][u][___], -1]]; While[incgs =!= {}, Dm = First[incgs]; incgs = Rest[incgs]; dif = DerivadaQ[Dm, Derivative[m][u][x]]; If[dif =!= False, subst = D[expr, dif]; dF = dF /. D[Derivative[m][u][x], dif] -> subst; incgs = Union[incgs, Cases[subst, Derivative[___][u][___], -1]]]]; dF]; DerivadaQ[Derivative[m___][V_][z___], Derivative[n___][V_][z___]]:= If[Apply[And, NonNegative[{m}-{n}]], Apply[Sequence, Transpose[{{z}, {m}-{n}}]], False] Protect[DReduce, DerivadaQ]; (* ============================================================ *) (* Página 6. Formato breve para derivadas *) (* ============================================================ *) Unprotect[DFormat, VarListD]; DFormat[True]:= Module[{}, Unprotect[Derivative]; Format[Literal[Derivative[nn__][ff_][vars:(_Symbol ..)]]]:= Apply[SequenceForm, {ff, Subscript[VarListD[{nn},{vars}]]}]; Format[Literal[Derivative[nn__][ff_][vars:(_Symbol ..)]], TeXForm]:= SequenceForm[ff,"_{", VarListD[{nn},{vars}],"}"]; Protect[Derivative]]; DFormat[False]:= Module[{}, Unprotect[Derivative]; Format[Derivative[nn__][ff_][vars:(_Symbol ..)]]=.; Format[Derivative[nn__][ff_][vars:(_Symbol ..)], TeXForm]=.; Protect[Derivative]]; VarListD[nn_List, vars_List]:= Module[{i, result = {}}, Do[If[nn[[i]]>0, result = Join[result, Table[vars[[i]], {nn[[i]]}]]], {i,1,Length[nn]}]; Apply[SequenceForm,result]]; Protect[DFormat, VarListD]; (* ============================================================ *) (* Página 7: Escritura de archivo TeX *) (* ============================================================ *) Unprotect[YaLieTeX, eqnarray, eqnarrayleft, eqnarrayright, Xi, Phi]; YaLieTeX[Archivo_, Eds_]:= YaLieTeX[Archivo, Eds, XiF, PhiF]; YaLieTeX[Archivo_, Eds_, LXi_List, LPhi_List]:= Module[{tmp}, EqsDs = DeleteCases[Eds, True]; If[VectorQ[EqsDs, (Head[#]===Equal)&], EqsDs = Map[First, Eds]]; LXiPhi = Join[LXi, LPhi]; Splice["yalie.mx", Archivo, FormatType->TeXForm, PageWidth->Infinity]; EqsDs=.; LXiPhi=.] Format[eqnarrayleft[l_], TeXForm]:= Module[{r}, r = Flatten[Map[({#, " & = & 0 ", " \\\\ ", "\n"})&, l]]; r = Apply[Sequence, Drop[r, -1]]; SequenceForm[r]] /; Length[l]>0 Format[eqnarrayleft[l_], TeXForm]:={} /; Length[l]==0 Format[eqnarrayright[l_], TeXForm]:= Module[{r}, r = Flatten[Map[({"& & ", #, " \\\\\n "})&, l]]; r = Apply[Sequence, Drop[r, -1]]; TeXFormInf[Off]; SequenceForm[r]] /; Length[l]>0 Format[eqnarrayrightCF[l_], TeXForm]:={} /; Length[l]==0 Format[eqnarray[l_], TeXForm]:= Module[{r}, r = Flatten[Map[({#[[1]], " & = & ", #[[2]], " \\\\\n "})&, l]]; r = Apply[Sequence, Drop[r, -1]]; SequenceForm[r]] /; Length[l]>0 Format[eqnarray[l_], TeXForm]:={} /; Length[l]==0 Format[Literal[Xi[i_]], TeXForm]:=StringReplace["\\xi^{(k)}", "k"->ToString[i]]; Format[Literal[Xi[i_][___]], TeXForm]:=StringReplace["\\xi^{(k)}", "k"->ToString[i]]; Format[Literal[Phi[i_]], TeXForm]:=StringReplace["\\phi^{(k)}", "k"->ToString[i]]; Format[Literal[Phi[i_][___]], TeXForm]:=StringReplace["\\phi^{(k)}", "k"->ToString[i]]; Protect[YaLieTeX, eqnarray, eqnarrayleft, eqnarrayright, Xi, Phi]; End[]; EndPackage[];