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$


$x_2 =y_3$


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$



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.

  • 1
    $\begingroup$ Could you give an example set of equations? $\endgroup$
    – Feyre
    Dec 3, 2016 at 13:42
  • 1
    $\begingroup$ @Feyre I added an example, thanks for helping. $\endgroup$ Dec 3, 2016 at 14:16
  • $\begingroup$ I hope you don't mind I put an undocumented tag on the question, since my answer explains an undocumented function. (The tag wiki says it's for questions about undocumented functions, but I thought it might help people find Internal`ProcessEquations`FirstOrderize, who were searching for something like it.) $\endgroup$
    – Michael E2
    Dec 6, 2016 at 1:43
  • 1
    $\begingroup$ Oh, no @MichaelE2 feel welcome to do this. I think your answer is very useful. I must add that Kuba's answer to the linked question above (use of the function NonlinearStateSpaceModel (or AffineSpaceStateModel), is valid only under the condition of linearizing the highest derivative (or the term regarding, eventually the input). I suppose your approach too, right? I think your answer is particularly valid as it helps getting through this beyond "control" theory. $\endgroup$ Dec 6, 2016 at 8:52
  • 2
    $\begingroup$ "I suppose your approach too, right?" -- I'm not sure if you're asking whether my approach linearizes, but just to be clear, it does not: firstOrderize[{x1'''[t]^2 == x2'[t] + x1[t], x2''[t]^3 == -x1[t]}, {x1, x2}, t]. $\endgroup$
    – Michael E2
    Dec 6, 2016 at 12:05

3 Answers 3



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]

Mathematica graphics

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

(*  {{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:

(*  {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,

 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 *)
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]] /. 


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

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


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_] :> 
      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

Mathematica graphics

  • 1
    $\begingroup$ May I ask what does the notatin x[2], y[1] mean? $\endgroup$ Dec 4, 2016 at 10:04
  • 1
    $\begingroup$ Well then in your output why is there a x[2]'? That would be a second order derivation. The code should output a set of first order differential equations. Is it that I'm misunderstanding your output? $\endgroup$ Dec 4, 2016 at 10:46
  • 1
    $\begingroup$ I'm still not convinced. The output of your to1storder[] are the last three equations. Yet why don't the parameters G and m appear in your solution? $\endgroup$ Dec 4, 2016 at 11:54
  • 2
    $\begingroup$ @MirkoAveta Please check my edit. $\endgroup$
    – xzczd
    Dec 4, 2016 at 12:28
  • 1
    $\begingroup$ You're my hero! I didn't know one could create options!!! This is amazing. $\endgroup$ Dec 4, 2016 at 12:30

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, 
   Flatten@{f /. subvar, 
     Table[D[sysvar[[i]], {t, 1}] - sysvar[[i + numbeq]], {i, 1, 
   , 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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.