Solución:
Introducción
Pensé que presentaría una forma de aprovechar las funciones integradas para hacer esto. He sabido desde hace mucho tiempo que NDSolve establece un problema de EDO como un sistema de ecuaciones de primer orden, por lo que el código básico debe estar ahí. Sin embargo, aparentemente no se encuentra fácilmente. Ahí está el tentadoramente nombrado Internal`ProcessEquations`FirstOrderize, que suena perfecto. Se necesitan de 4 a 6 argumentos y, finalmente, adiviné cómo configurar el problema del OP.
fop = Internal`ProcessEquations`FirstOrderize[Thread[f1 == 0], {t}, 1, {x1, x2}];
Column[fop, Dividers -> All]
El nuevo sistema se da uniendo los dos primeros componentes de la salida:
newsys = Join @@ fop[[1 ;; 2]]
Las nuevas variables se almacenan en
fop[[3]]
(* {{x1, NDSolve`x1$56$1}, {x2, NDSolve`x2$84$1}} *)
Y la relación con el problema original se da en las reglas del último componente:
fop[[4]]
(* {NDSolve`x1$56$1 -> Derivative[1][x1], NDSolve`x2$84$1 -> Derivative[1][x2]} *)
Si no te gusta el NDSolve` variables del módulo, hay otra utilidad que se puede descubrir a través de [email protected]"VariableTransformation" en el NDSolve componentes de datos estatales. No hay documentación, AFAIK, pero puede generar ejemplos evaluando NDSolve`StateData objetos. Su forma es
Internal`ProcessEquations`FirstOrderReplace[expr, indepVars, n, depVars, newVars]
(Solo he visto n == 1 en la tercera posición en una EDO.) Por ejemplo, para el sistema OP,
Internal`ProcessEquations`FirstOrderReplace[
Thread[f1 == 0], {t}, 1, {x1, x2}, {{X, XP}, {Y,YP}}]
(*
{ff1 + d11 X[t] + b11 XP[t] + d12 Y[t] + c11 XP[t] Y[t] + b12 YP[t] +
c12 X[t] YP[t] + a11 XP'[t] + a12 YP'[t] == 0,
ff2 + d21 X[t] + b21 XP[t] + d22 Y[t] + c21 XP[t] Y[t] + b22 YP[t] +
c22 X[t] YP[t] + a21 XP'[t] + a22 YP'[t] == 0}
*)
Tenga en cuenta la lista de listas en los nuevos nombres de variables. Cada instancia de x1 y x1' es reemplazado por X y XP respectivamente; x1'' es reemplazado por XP'. Lo mismo ocurre con la otra variable x2.
Función de propósito general
Aquí hay una función que permite cambiar el nombre de las variables como desee el OP. Es un poco complicado no cambiar el nombre de los lados izquierdos con FirstOrderReplace; Lo hago inactivando Derivative temporalmente.
(* With arbitrary symbol renaming *)
ClearAll[firstOrderize];
Options[firstOrderize] = {"NewSymbolGenerator" -> (Unique["y"] &)};
firstOrderize[sys_, vars_, t_, OptionsPattern[]] :=
Module[{fop, newsym, toNewVar},
newsym = OptionValue["NewSymbolGenerator"];
fop = Internal`ProcessEquations`FirstOrderize[sys, {t}, 1, vars];
If[newsym === Automatic,
(* don't rename *)
Fla[email protected] fop[[1 ;; 2]],
(* rename *)
toNewVar = With[{newvars = MapIndexed[newsym, fop[[3]], {2}]},
Internal`ProcessEquations`FirstOrderReplace[#, {t}, 1, vars, newvars] &];
[email protected] {toNewVar[fop[[1]] /. Last[fop]],
Activate[toNewVar[Inactivate[[email protected][[2]], Derivative]] /.
toNewVar[fop[[4]]]]}
]
]
Ejemplos de
Ejemplo de OP: la función de cambio de nombre automático utiliza Unique["y"], que agregará un número a "y", el siguiente número.
firstOrderize[Thread[f1 == 0], {x1, x2}, t]
(*
{ff1 + d11 y3[t] + b11 y4[t] + d12 y5[t] + c11 y4[t] y5[t] +
b12 y6[t] + c12 y3[t] y6[t] + a11 y4'[t] + a12 y6'[t] == 0,
ff2 + d21 y3[t] + b21 y4[t] + d22 y5[t] + c21 y4[t] y5[t] +
b22 y6[t] + c22 y3[t] y6[t] + a21 y4'[t] + a22 y6'[t] == 0,
y3'[t] == y4[t],
y5'[t] == y6[t]}
*)
Uno puede usar la opción "NewSymbolGenerator" para especificar cómo desea que se generen los símbolos. Debería ser una función, que se aplicará a la NDSolve variables en fop[[3]] con MapIndexed a nivel {2}.
firstOrderize[{x1'[t]^2 == x2'[t] + x1[t], x2''[t] == -x1[t]}, {x1, x2}, t,
"NewSymbolGenerator" -> (Symbol[{"a", "b"}[[[email protected]#2]] <> [email protected]@#2] &)]
(* {a1'[t]^2 == a1[t] + b2[t], b2'[t] == -a1[t], b1'[t] == b2[t]} *)
Una forma de obtener la numeración de 1 para 4 cada vez:
Module[{n = 0},
firstOrderize[{x1''[t] == x2[t] + x1[t], x2''[t] == -x1[t]}, {x1, x2}, t,
"NewSymbolGenerator" -> (Symbol["y" <> ToString[++n]] &)]
] // Sort
(*
{y1'[t] == y2[t],
y2'[t] == y1[t] + y3[t],
y3'[t] == y4[t],
y4'[t] == -y1[t]}
*)
Este es mi enfoque, creo que es más ordenado y más general:
[email protected]
Options[to1storder] = {"form" -> (#[#2] &)};
to1storder[eq_List, func_List, argu_, OptionsPattern[]] :=
Module[{maxorder, mapthread, lhsae, lhsde, detoae},
maxorder = #[[[email protected][#, -1]]] &@
[email protected][eq, Derivative[i_][#][_] :> i, Infinity] & /@ func;
mapthread = MapThread[#, {func, maxorder}] &;
lhsae = mapthread[#[#2]@argu &];
lhsde = mapthread[#[#2 - 1]'@argu &];
detoae = ((f : Alternatives @@ func) (i_: 0) | Derivative[i_][f_])[a_] :> f[i][a];
{((*Solve[*)eq /. detoae(*,#]*)/. Thread[# -> lhsde](* /. Rule -> Equal*)) &@lhsae,
mapthread[Table[#[n - 1]'@argu == #[n]@argu, {n, 1, #2 - 1}] &]} /. (f :
Alternatives @@ func)[i_] :> OptionValue["form"][f, i]]
to1storder[eq_, func_, argu_, o : OptionsPattern[]] :=
to1storder[[email protected]{eq}, [email protected]{func}, argu, o];
Si desea colocar el término con el orden de derivada más alto al lado izquierdo de la ecuación, simplemente agregue el código en la anotación.
Puede encontrar la definición de detoae un poco confuso. Para entenderlo, solo nota que b /. (n_: 0) b -> n devoluciones 0. Para más información, consulte el documento de Default.
El siguiente es un ejemplo de que su OrderReduce no se puede manejar correctamente:
test = With[{x = x[t], y = y[t]}, {D[x, t, t, t]^2 == -((G m x)/(x^2 + y^2)^(3/2)),
D[y, t, t] == -((G m y)/(x^2 + y^2)^(3/2))}]
to1storder[test, {x, y}, t, "form" -> Subscript] // Flatten // TableForm
He escrito el siguiente código para este problema. La función lee como entradas el lado derecho de la ecuación en forma implícita y una lista de sus funciones desconocidas asignando explícitamente su dependencia a la variable de tiempo. Cualquier sugerencia para hacer esto más ordenado o eficiente es muy bienvenida.
QDim[a_, b_] := TrueQ[Length[a] == Length[b]]
OrderReduce[f__, var__] := With[{}, If[QDim[f, var] == True,
numbeq = Length[f];
dimsys = 2 numbeq;
sysvar = Array[Subscript[nvar, #][t] &, {dimsys}];
syscom =
[email protected]{x, Table[D[var[[i]], {t, 1}], {i, 1, numbeq}],
Table[D[var[[i]], {t, 2}], {i, 1, numbeq}]};
subvar =
[email protected]{Table[syscom[[i]] -> sysvar[[i]], {i, 1, dimsys}],
Table[D[var[[i]], {t, 2}] ->
D[sysvar[[numbeq + i]], {t, 1}], {i, 1, numbeq}],
Table[D[var[[i]], {t, 1}] -> sysvar[[numbeq + i]], {i, 1,
numbeq}]};
[email protected]{f /. subvar,
Table[D[sysvar[[i]], {t, 1}] - sysvar[[i + numbeq]], {i, 1,
numbeq}]}
, Print["Error: dimensional mismatch."]]]
Si aplico esto al ejemplo anterior,
f1 = {a11 D[x1[t], {t, 2}] + a12 D[x2[t], {t, 2}] +
b11 D[x1[t], {t, 1}] + b12 D[x2[t], {t, 1}] +
c11 D[x1[t], {t, 1}] x2[t] + c12 x1[t] D[x2[t], {t, 1}] +
d11 x1[t] + d12 x2[t] + ff1,
a21 D[x1[t], {t, 2}] + a22 D[x2[t], {t, 2}] +
b21 D[x1[t], {t, 1}] + b22 D[x2[t], {t, 1}] +
c21 D[x1[t], {t, 1}] x2[t] + c22 x1[t] D[x2[t], {t, 1}] +
d21 x1[t] + d22 x2[t] + ff2};
OrderReduce[f1, {x1[t], x2[t]}] // TableForm
Obtengo correctamente:
Como dije, ¡cualquier alternativa mejor es bienvenida!