how to solve a non linear system of ODEs

Question

I define this Lagrangian describing a physical system under investigation (two pendula connected by non-linear springs):

    Lag = 1/2 m (x1'[t]^2 + x2'[t]^2) - (m g)/(
    2 l) (x1[t]^2 + x2[t]^2) - (k/
      2 (x1[t]^2 + (x2[t] - x1[t])^2 + (x2[t] - x1[t])^2 + (L - 
      x2[t])^2) + α (x1[
     t]^3 + (x2[t] - x1[t])^3 + (x2[t] - x1[t])^3 + (L - 
      x2[t])^3) + β (x1[
     t]^4 + (x2[t] - x1[t])^4 + (x2[t] - x1[t])^4 + (L - 
      x2[t])^4));

where all coefficients are constant. I ask Mathematica to calculate the Euler equations as follows

    ee = EulerEquations[Lag, {x1[t], x2[t]}, t]

Hence I ask Mathematica to solve the corresponding non linear system of two second order ODEs, using NDSolve

    sol = NDSolve[
    Join[ee, {x1[0] == 0, x2[0] == 0, x1'[0] == 0.1, 
    x2'[0] == 0}], {x1[t], x2[t]}, {t, 0, 10000}][[1, All, 2]]

NDSolve works pretty fast. But after that I need to calculate the average energy of each body as a function of time and I use those two list of lists:

    u1=Table[NIntegrate[(m g)/(2 l) sol[[1]]^2 + 
    k/2 (sol[[1]]^2 + (sol[[2]] - sol[[1]])^2) + α (sol[[
    1]]^3 + (sol[[2]] - sol[[1]])^3) + β (sol[[
    1]]^4 + (sol[[2]] - sol[[1]])^4), {t, 0, i}], {i, n}]];
    uMedia1=Accumulate[u1];

It takes too much to plot sol[[1]] (or 2) and even longer to NIntegrate over each time interval (by the way Accumulate is very fast). As a consequence I guess I'm doing something wrong with that.

Is this the correct procedure to solve this problem? I appreciate your help. Marta

k = 1; l = 1; L = 2; m = 1; g = 9.8; α = 0.5; β = 0; λ = 0.1;

Answer 1

You can speed it up doing the integration with NDSolve. This has to be done only once to get an interpolation-function of the average energy.

ndsol2 = NDSolve[{en'[t] == (m g)/(2 l) sol[[1]]^2 + 
            k/2 (sol[[1]]^2 + (sol[[2]] - 
       sol[[1]])^2) + \[Alpha] (sol[[1]]^3 + (sol[[2]] - 
       sol[[1]])^3) + \[Beta] (sol[[1]]^4 + (sol[[2]] - 
       sol[[1]])^4), en[0] == 0}, en, {t, 0, 10000}, 
       MaxSteps -> \[Infinity]]

(*     {{en -> InterpolatingFunctionx[{0.,10000.}},<>]}}   *)

In order to get the average energy, I think you have to divide by t

Plot[(en[t]/t) /. First@ndsol2, {t, 0, 100}]

But since the x1 and x2 are highly oscillatory, I think you have to use higher workingprecision to get reliable results up to t=10000.

Answer 2

A tweak on @Akku14's nice answer:

It seems more efficient to simultaneously solve for x1, x2 and en like this:

sol2 = NDSolve[
  Join[ee, {x1[0] == 0, x2[0] == 0, x1'[0] == 0.1, x2'[0] == 0},
  {en'[t] == (m g)/(2 l) x1[t]^2 + k/2 (x1[t]^2 + (x2[t] - x1[t])^2) +
   α (x1[t]^3 + (x2[t] - x1[t])^3) + β (x1[t]^4 + (x2[t] - x1[t])^4), en[0] == 0}
  ], {x1, x2, en}, {t, 0, 10000}][[1]]

Takes 1.2 sec vs 0.98 sec (sol) plus 6.02 sec (ndsol2) doing them separately. It might also be more accurate, since it avoids the InterpolatingFunction from sol.