Differential equations substituting list back in solution

I'm trying to solve some Lagrange differential equations of motion, where $$L = T - U$$.

How can I use the solution for each x1, x2, x3, ..., xn?

Clear["Global*"]
Needs["VariationalMethods"]
n = 10;
t0 = 10;
Subscript[x, 0][t_] := 0;
Subscript[x, n + 1][t_] := 0;
m = 1;
v0 = 0.1;
Table[Subscript[k, j] = 1, {j, 0, n + 1}];
ic = Table[{Subscript[x, j][0] == 0, Derivative[1][Subscript[x, j]][0] ==RandomReal[{-v0, v0}]}, {j, 1, n}];

ue[x_, t_, k_, n_] :=
(1/2)*Sum[Subscript[k, j]*(Subscript[x, j][t] - Subscript[x, j - 1][t])^2,
{j, 1, n + 1}];
te[x_, t_, n_] := (1/2)*m*Sum[Derivative[1][Subscript[x, j]][t]^2, {j, 1, n}];

lg[x_, t_, k_, n_] := te[x, t, n] - ue[x, t, k, n];

eqm[j_] := D[lg[x, t, k, n], {Subscript[x, j][t], 1}] -
D[D[lg[x, t, k, n], {Derivative[1][Subscript[x, j]][t], 1}], {t, 1}];

• 1) There was a question about frequencies, not exact solutions. 2) There are two types of boundary conditions with $k_0=0; 5$. Mar 10, 2019 at 18:41
• Please look at the instructions here to copy readable code from MMA: How to copy code from MMA so it looks good on this site Mar 10, 2019 at 20:58

Numerical solution

Clear["Global*"]
Needs["VariationalMethods"]
n = 10; t0 = 10; Subscript[x, 0][t_] := 0; Subscript[x, n + 1][t_] := 0; m = 1; v0 = 0.1;
Table[Subscript[k, j] = 1, {j, 0, n + 1}];
ic = Table[{Subscript[x, j][0] == 0, Derivative[1][Subscript[x, j]][0] ==
RandomReal[{-v0, v0}]}, {j, 1, n}];
ue[x_, t_, k_, n_] :=
(1/2)*Sum[Subscript[k, j]*(Subscript[x, j][t] - Subscript[x, j - 1][t])^2,
{j, 1, n + 1}];
te[x_, t_, n_] := (1/2)*m*Sum[Derivative[1][Subscript[x, j]][t]^2, {j, 1, n}];
lg[x_, t_, k_, n_] := te[x, t, n] - ue[x, t, k, n];
eqm[j_] := D[lg[x, t, k, n], {Subscript[x, j][t], 1}] -
D[D[lg[x, t, k, n], {Derivative[1][Subscript[x, j]][t], 1}], {t, 1}];
sol = NDSolveValue[Flatten[{Table[eqm[i] == 0, {i, 1, n}], ic}],
Table[Subscript[x, i][t], {i, 1, n}], {t, 0, t0}]
Plot[sol, {t, 0, t0}, PlotLegends -> Automatic]


Calculation of frequencies for a system of 100 particles with springs of different stiffness

    ks1[i_] := 3 + 2*(-1)^i
ks2[i_] := 3 - 2*(-1)^i
m = 1;
omega[k_, n_] :=
Block[{k0 = k, P = n},
eq1 = -m*(x[1]'')[t] + k0*(x[0][t] - x[1][t]) -
5 (x[1][t] - x[2][t]);
eqP = -m*(x[P]'')[t] + k0*(x[P + 1] - x[P][t]) -
5 (-x[P - 1][t] + x[P][t]);
eqn = Table[-m*x[i]''[t] - ks1[i]*(x[i][t] - x[i - 1][t]) -
ks2[i]*(x[i][t] - x[i + 1][t]), {i, 2, P - 1}];
Eq1 = Join[{eq1, eqn, eqP}] /. {x[0][t] -> 0, x[P + 1][t] -> 0};
op = Table[x[i_][t_] := A[i]*Exp[I*om*t], {i, 1, P}];
Eq2 = Flatten[Eq1 /. t -> 0];
Eq = Table[Eq2[[i]] == 0, {i, 1, P}];
matr = CoefficientArrays[Eq, Table[A[i], {i, 1, P}]];
M = Last[Normal[matr]];
f = Det[M];
s = Solve[f == 0, om];
lst = Table[N[om /. s[[i]], 10], {i, 1, Length[s]}];
Drop[Sort[lst], P]]


An example of using the function omega[]

lst0 = omega[0, 100];
lst5 = omega[5, 100];
{ListPlot[Abs[lst0 - lst5], PlotRange -> All],
ListPlot[{lst0, lst5},
PlotLegends -> {Row[{"k0 = ", 0}], Row[{"k0 = ", 5}]}],
ListPlot[lst0, PlotLabel -> Row[{"k0 = ", 0}]],
ListPlot[lst5, PlotLabel -> Row[{"k0 = ", 5}]]}


• @Alrubaie I fixed the code. Mar 12, 2019 at 21:11