Skip to main content

differential equations - Reduce ODE into 1st order


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



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]


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


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

Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...