5
$\begingroup$

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;

allEq = Join[oriEq, addEq];

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

The final single ODE I expected

$\endgroup$
10
  • $\begingroup$ What are you trying to obtain? Anyway, the allEq involves 2 oriEq, which looks like a mistake. $\endgroup$
    – xzczd
    Jan 6 at 15:22
  • $\begingroup$ @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. $\endgroup$
    – xinxin guo
    Jan 6 at 15:46
  • $\begingroup$ @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. $\endgroup$
    – xinxin guo
    Jan 6 at 16:17
  • $\begingroup$ @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. $\endgroup$
    – xinxin guo
    Jan 7 at 2:31
  • $\begingroup$ So, you're only interested in solving the problem with Eliminate? $\endgroup$
    – xzczd
    Jan 7 at 5:29
4
$\begingroup$

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;
addDPolys = 
  Table[D[origDPolys, {t, i}], {i, n}] // Simplify // Flatten;
allDPolys = Join[origDPolys, addDPolys];
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.

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – xinxin guo
    Jan 8 at 1:20
  • $\begingroup$ Er… why does this work? $\endgroup$
    – xzczd
    Jan 8 at 2:46
  • 1
    $\begingroup$ @xzczd It's operating under the PRinciples of Advanced MAgic (PRAMA). I had thought that was obvious. $\endgroup$ Jan 8 at 16:33
  • 2
    $\begingroup$ 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. $\endgroup$ Jan 8 at 16:37
2
$\begingroup$

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} *)
$\endgroup$
1
  • $\begingroup$ Yeah, I think your code just caputrues the essence of problem I am sovling. It is the eqsum is the only key. $\endgroup$
    – xinxin guo
    Jan 7 at 9:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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