# How to guide Eliminate or GroebnerBasis to reduce a set of simple odes to a single ode (which can be done by Laplace Transform)

Recently, I am trying to use Eliminate or GroebnerBasis to simplify a system of ODEs. I don't want the solution of ODEs. What I expected is the final single ODE, as shown in red box in the figure below after simplify odes. The figure was borrowed from here.

I should emphasize that the final single ode I desired is a classical conclusion in the theory of rheology, which can obtained using Laplace transform technique.

The number of odes is conroled by parameter n, i.e., n=2,3,4,5,...

I can get the final single ODE only when n=2, and the computation seems take forever when n=3. How to get the final single ODE for larger value of n = 3,4,5,...?

Below is my code. Anys suggestions are appreciated!

Remove["Global*"] // Quiet;

n = 2;

oriEq = {
Table[{s[i][t] == k[i] e[i][t] + \[Eta][i] e[i]'[t],
s[t] == s[i][t]}, {i, n}],
e[t] == Sum[e[i][t], {i, n}]
} // Flatten;

addEq = Table[D[oriEq, {t, i}], {i, n}] // Flatten;

keep0 = Table[s[t], {i, n}];
keep = {keep0, Table[D[keep0, {t, i}], {i, n}]} // Flatten;

dele0 = Table[{s[i][t], e[i][t]}, {i, n}];
dele = {dele0, Table[D[dele0, {t, i}], {i, n + 1}]} // Flatten

solAnd = Eliminate[allEq, dele] // Simplify

gb = GroebnerBasis[allEq, keep, dele] // Simplify


• What are you trying to obtain? Anyway, the allEq involves 2 oriEq, which looks like a mistake. Jan 6 at 15:22
• @xzczd Thanks for your attention. I corrected the code, but it still doesn't work. I am trying to derive odes of generalized kelvin voigt model. This was explained in this lecture m.youtube.com/watch?v=0cjKIPj5cNs. Jan 6 at 15:46
• @xzczd The procedure of deriving odes of generalized kelvin-voit model is the same as that of generalized maxwell model which is shown here en.wikipedia.org/wiki/Generalized_Maxwell_model. You can find that a single ode can be obtained. Jan 6 at 16:17
• @xzczd The general form of single ode was add in my post. I should emphasize that what I desired is a classical result concluded in the theory of rheology. Thanks :) Have nice day sir. Jan 7 at 2:31
• So, you're only interested in solving the problem with Eliminate? Jan 7 at 5:29

It is a matter of setting options to help GroebnerBasis. I show below for the case n=4.

n = 4;
origDPolys = {Table[{s[i][t] - (k[i] e[i][t] + nu[i] e[i]'[t]),
s[t] - s[i][t]}, {i, n}], e[t] - Sum[e[i][t], {i, n}]} //
Flatten;
Table[D[origDPolys, {t, i}], {i, n}] // Simplify // Flatten;
keep0 = s[t];
keep = {keep0, Table[D[keep0, {t, i}], {i, n}]} // Flatten;
dele0 = Table[{s[i][t], e[i][t]}, {i, n}];
dele = {dele0, Table[D[dele0, {t, i}], {i, n + 1}]} // Flatten;

In[49]:= Timing[
gb = GroebnerBasis[allDPolys, keep, dele,
MonomialOrder -> EliminationOrder,
CoefficientDomain -> RationalFunctions]]
(* Out[49]= {0.737172, {(-e[t])*k[1]*k[2]*k[3]*
k[4] + (k[1]*k[2]*k[3] + k[1]*k[2]*k[4] + k[1]*k[3]*k[4] +
k[2]*k[3]*k[4])*s[t] -
k[2]*k[3]*k[4]*nu[1]*Derivative[1][e][t] -
k[1]*k[3]*k[4]*nu[2]*Derivative[1][e][t] -
k[1]*k[2]*k[4]*nu[3]*Derivative[1][e][t] -
k[1]*k[2]*k[3]*nu[4]*Derivative[1][e][t] +
(k[2]*k[3]*nu[1] + k[2]*k[4]*nu[1] + k[3]*k[4]*nu[1] +
k[1]*k[3]*nu[2] + k[1]*k[4]*nu[2] + k[3]*k[4]*nu[2] +
k[1]*k[2]*nu[3] + k[1]*k[4]*nu[3] + k[2]*k[4]*nu[3] +
k[1]*k[2]*nu[4] + k[1]*k[3]*nu[4] + k[2]*k[3]*nu[4])*
Derivative[1][s][t] -
k[3]*k[4]*nu[1]*nu[2]*Derivative[2][e][t] -
k[2]*k[4]*nu[1]*nu[3]*Derivative[2][e][t] -
k[1]*k[4]*nu[2]*nu[3]*Derivative[2][e][t] -
k[2]*k[3]*nu[1]*nu[4]*Derivative[2][e][t] -
k[1]*k[3]*nu[2]*nu[4]*Derivative[2][e][t] -
k[1]*k[2]*nu[3]*nu[4]*Derivative[2][e][t] +
(k[3]*nu[1]*nu[2] + k[4]*nu[1]*nu[2] + k[2]*nu[1]*nu[3] +
k[4]*nu[1]*nu[3] + k[1]*nu[2]*nu[3] + k[4]*nu[2]*nu[3] +
k[2]*nu[1]*nu[4] + k[3]*nu[1]*nu[4] + k[1]*nu[2]*nu[4] +
k[3]*nu[2]*nu[4] + k[1]*nu[3]*nu[4] + k[2]*nu[3]*nu[4])*
Derivative[2][s][t] -
k[4]*nu[1]*nu[2]*nu[3]*Derivative[3][e][t] -
k[3]*nu[1]*nu[2]*nu[4]*Derivative[3][e][t] -
k[2]*nu[1]*nu[3]*nu[4]*Derivative[3][e][t] -
k[1]*nu[2]*nu[3]*nu[4]*Derivative[3][e][t] +
(nu[1]*nu[2]*nu[3] + nu[1]*nu[2]*nu[4] + nu[1]*nu[3]*nu[4] +
nu[2]*nu[3]*nu[4])*Derivative[3][s][t] -
nu[1]*nu[2]*nu[3]*nu[4]*Derivative[4][e][t]}} *)


At n=5 the timing is over 5 seconds, so I would not expect this to work for much larger values.

• This is amazing. Thansk for your help. By the way,it seems that GroebnerBasis  is generally more powerful than Eliminate, because it has more options to tune. Jan 8 at 1:20
• Er… why does this work? Jan 8 at 2:46
• @xzczd It's operating under the PRinciples of Advanced MAgic (PRAMA). I had thought that was obvious. Jan 8 at 16:33
• Re Eliminate vs GroebnerBasis. The former is built on some really old code infrastructure, the latter is more modern albeit also perhaps in need of refinement at this point. There are some subtleties that make it nontrivial to rewrite Eliminate using GroebnerBasis, so I just try to promote the one and downplay the other. Jan 8 at 16:37

Well, I've no idea about how to speed up Eliminate, but the following is my work-around, still based on the idea of differentiation and equation solving, just with a bit of manual analysis:

Clear[Derivative, n, s];
s[i_][t] = s[t];
Derivative[n_][e[i_]][t] = D[(s[i][t] - k[i] e[i][t])/η[i], {t, n - 1}];
eqsum[n_] = e[t] == Sum[e[i][t], {i, n}];
n = 5;
tst = D[eqsum[n], {t, n}] /.
Solve[Table[D[eqsum@n, {t, index}], {index, 0, n - 1}], e[#][t] & /@
Range[n]][[1]](*// Simplify*); // AbsoluteTiming
(* {4.53276, Null} *)

• Yeah, I think your code just caputrues the essence of problem I am sovling. It is the eqsum` is the only key. Jan 7 at 9:58