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) + \[Alpha]α (x1[
t]^3 + (x2[t] - x1[t])^3 + (x2[t] - x1[t])^3 + (L -
x2[t])^3) + \[Beta]β (x1[
t]^4 + (x2[t] - x1[t])^4 + (x2[t] - x1[t])^4 + (L -
x2[t])^4));
where all coefficients are constant. I ask MathematicaMathematica to calculate the Euler equations as follows
ee = EulerEquations[Lag, {x1[t], x2[t]}, t]
Hence I ask MMathematica to solve the corresponding non linear system of two second order ODEs, using NDSolveNDSolve
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]]
NDSolveNDSolve 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) + \[Alpha]α (sol[[
1]]^3 + (sol[[2]] - sol[[1]])^3) + \[Beta]β (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; \[Alpha]α = 0.5; \[Beta]β = 0; \[Lambda]λ = 0.1 ;1;
