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

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1.