# how to solve a non linear system of ODEs

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, x2 == 0, x1' == 0.1,
x2' == 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[]^2 +
k/2 (sol[]^2 + (sol[] - sol[])^2) + α (sol[[
1]]^3 + (sol[] - sol[])^3) + β (sol[[
1]]^4 + (sol[] - sol[])^4), {t, 0, i}], {i, n}]];
uMedia1=Accumulate[u1];


It takes too much to plot sol[] (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;

• Could you give parameter values? Also, nota bene for others: EulerEquations requires Needs["VariationalMethods"]. Nov 5, 2017 at 18:48
• Could you give us your value for n and clarify what you are doing with u1 and uMedia1? If you are just trying to calculate and plot the average potential energy you need only do something like u1avg[tf_]:=NIntegrate[PotEnergy1[t],{t,0,tf}]/tf where PotEnergy1 is the argument of your Nintegrate' statement in u1. You can then Plot this. If you are trying to do something different please tell us.
– qbit
Nov 6, 2017 at 5:43

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[]^2 +
k/2 (sol[]^2 + (sol[] -
sol[])^2) + \[Alpha] (sol[]^3 + (sol[] -
sol[])^3) + \[Beta] (sol[]^4 + (sol[] -
sol[])^4), en == 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.

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, x2 == 0, x1' == 0.1, x2' == 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}
], {x1, x2, en}, {t, 0, 10000}][]


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.

• I appreciate your help. I'm following the last suggestion (solving for x1, x2, en1 and en2 at the same time) and it look very fast. Nov 7, 2017 at 1:37
• I appreciate your help. I'm following the last suggestion (solving for x1, x2, en1 and en2 at the same time) and it look very fast. But now I'm not able to calculate the Variance of (en1,en2) as a function of time. Could you please help? Which is the correct syntax? I give the details for en2 editing my primary post. Nov 7, 2017 at 1:50