Reduce ODE into 1st order

Question

I have an explicit set of differential equations:

$ \ddot{x}=f(x,\dot{x})$

I would like to reduce it in the following way:

$ \dot{y} = g(y)$

by substitutions as shown here: wikipedia. I have done this heuristically for a certain system, using rules and CoefficientArrays. Yet I would like to have a function that would make up the new family of unknown functions automatically, and get this reduction done in an efficient and clean way -whatever may be the number of equations, i.e. the number of unknown functions. I don't need it, yet it could be useful to also try and develop a code for a differential equations of higher order.

Example (2 equations)

$a_{11} \ddot{x}_1+a_{12} \ddot{x}_2 +b_{11} \dot{x}_1 +b_{12} \dot{x}_2 +c_{11} \dot{x}_1 x_2 +c_{12}x_1 \dot{x}_2+d_{11}x_1+d_{12}x_2+f_1=0 $ $a_{21} \ddot{x}_1+a_{22} \ddot{x}_2 +b_{21} \dot{x}_1 +b_{22} \dot{x}_2 +c_{21} \dot{x}_1 x_2 +c_{22}x_1 \dot{x}_2+d_{21}x_1+d_{22}x_2+f_2=0 $

(* Corresponding code *)
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};

We define:

$x_1 =y_1$

$\dot{x}_1=y_2$

$x_2 =y_3$

$\dot{x}_2=y_4$

And obtain the following first order system of equations:

$a_{11} \dot{y}_2+a_{12} \dot{y}_4 +b_{11} y_2 +b_{12} y_4 +c_{11} y_2 y_3 +c_{12}y_1 y_4+d_{11}y_1+d_{12}y_3+f_1=0 $

$a_{21} \dot{y}_2+a_{22} \dot{y}_4 +b_{21} y_2 +b_{22} y_4 +c_{21} y_2 y_3 +c_{22}y_1 y_4+d_{21}y_1+d_{22}y_3+f_2=0$

$\dot{y}_1=y_2$

$\dot{y}_3=y_4$

UPDATE: Mathematica already has a function for this

Indeed, I found out this is a duplicate question here, and that indeed Mathematica already has a function for this. Yet I would like to not delete this having xzczd worked on it and, especially, having received such an interesting answer. I consider his answer a real lesson of coding and thank him again.

Answer 1

Introduction

I thought I would present a way to leverage built-in functions to do this. I've known for a long time that NDSolve sets up an ODE problem as system of first-order equations, so the basic code must be in there. Apparently, though, it's not easily found. There is the tantalizingly named Internal`ProcessEquations`FirstOrderize, which sounds perfect. It takes 4 to 6 arguments, and eventually I guessed how to set up the OP's problem.

fop = Internal`ProcessEquations`FirstOrderize[Thread[f1 == 0], {t}, 1, {x1, x2}]; 
Column[fop, Dividers -> All]

The new system is given by joining the first two components of the output:

newsys = Join @@ fop[[1 ;; 2]]

The new variables are stored in

fop[[3]]
(*  {{x1, NDSolve`x1$56$1}, {x2, NDSolve`x2$84$1}}  *)

And the relationship to the original problem is given in the rules in the last component:

fop[[4]]
(*  {NDSolve`x1$56$1 -> Derivative[1][x1], NDSolve`x2$84$1 -> Derivative[1][x2]}  *)

If you don't like the NDSolve` module variables, there is another utility that can be discovered via state@"VariableTransformation" in the NDSolve state data components. There's no documentation, AFAIK, but you can generate examples by evaluating NDSolve`StateData objects. Its form is

Internal`ProcessEquations`FirstOrderReplace[expr, indepVars, n, depVars, newVars]

(I've only ever seen n == 1 in the third position in an ODE.) For example, for the OP's system,

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}
*)

Note the list of lists in the new variable names. Each instance of x1 and x1' is replaced by X and XP respectively; x1'' is replaced by XP'. Likewise for the other variable x2.

General-purpose function

Here is a function that allows renaming of the variables as the OP wishes. It's a bit tricky not to rename the left hand sides with FirstOrderReplace; I do it by inactivating Derivative temporarily.

(* 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 *)
   Flatten@ fop[[1 ;; 2]],
   (* rename *)
   toNewVar = With[{newvars = MapIndexed[newsym, fop[[3]], {2}]},
     Internal`ProcessEquations`FirstOrderReplace[#, {t}, 1, vars, newvars] &];
   Flatten@ {toNewVar[fop[[1]] /. Last[fop]],
     Activate[toNewVar[Inactivate[Evaluate@fop[[2]], Derivative]] /. 
       toNewVar[fop[[4]]]]}
   ]
  ]

Examples

OP's Example: The automatic renaming function uses Unique["y"], which will add a number to "y", whatever number is next.

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]}
*)

One can use the option "NewSymbolGenerator" to specify how you want the symbols generated. It should be a function, which will be applied to the NDSolve variables in fop[[3]] with MapIndexed at level {2}.

firstOrderize[{x1'[t]^2 == x2'[t] + x1[t], x2''[t] == -x1[t]}, {x1, x2}, t,
 "NewSymbolGenerator" -> (Symbol[{"a", "b"}[[First@#2]] <> ToString@Last@#2] &)]
(*  {a1'[t]^2 == a1[t] + b2[t], b2'[t] == -a1[t], b1'[t] == b2[t]}  *)

One way to get numbering from 1 to 4 every time:

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]}
*)

Answer 2

Here's my approach, I think it's tidier and more general:

ClearAll@to1storder

Options[to1storder] = {"form" -> (#[#2] &), "standardize" -> False};

to1storder[eq_List, func_List, argu_, OptionsPattern[]] := 
 Module[{maxorder, mapthread, lhsae, lhsde, detoae}, 
  maxorder = #[[First@Ordering[#, -1]]] &@
      Union@Cases[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];
  {If[OptionValue["standardize"], 
      Solve[#, lhsae] /. Rule -> Equal // Flatten, #] &[
      eq /. detoae] /. Thread[lhsae -> lhsde], 
    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[Flatten@{eq}, Flatten@{func}, argu, o];

If you want to put the term with highest order of derivative to the left hand side of the equation, set option "standardize" -> True.

You may find the definition of detoae a little confused. To understand it, just notice that b /. (n_: 0) b -> n returns 0. For more information, check the document of Default.

The following is an example that your OrderReduce can't handle correctly:

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

Answer 3

I have written the following code for this problem. The function reads as inputs the right-hand-side of the equation in implicit form and a list of it's unknown functions explicitly assigning their dependency to the time variable. Any suggestion to make this more tidy or efficient is highly welcomed.

    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 = 
    Flatten@{x, Table[D[var[[i]], {t, 1}], {i, 1, numbeq}], 
      Table[D[var[[i]], {t, 2}], {i, 1, numbeq}]};
   subvar = 
    Flatten@{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}]};
   Flatten@{f /. subvar, 
     Table[D[sysvar[[i]], {t, 1}] - sysvar[[i + numbeq]], {i, 1, 
       numbeq}]}
   , Print["Error: dimensional mismatch."]]]

If I apply this to the example above,

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

I correctly get:

As I said, any better alternative is highly welcomed!