Hello. when I run the code below, i get the error NDSolveValue::reinitfail: Unable to reinitialize the system at t = 11.692137665979317` within specified tolerances. which i think is due to numerical instability in the calculations. any help on how to solve this issue would be greatly appreciated.
ClearAll["Global`*"]
m[n_] := Piecewise[{{0.1, n == 40}}, 0.0001];
l[n_] := 0.02;
\[Zeta][t_] := 0.1 Cos[7 t]
r[n_, t_] := \[Zeta][t] UnitVector[3, 3] +
Sum[UnitVector[3, 1]*l[i] Sin[\[Theta][i, t]] -
UnitVector[3, 3]*l[i] Cos[\[Theta][i, t]], {i, n}];
rdot[n_, t_] := D[r[n, t], t];
tt[nn_] := 1/2 Sum[m[n]*Total[rdot[n, t]^2], {n, nn}];
vv[nn_] := g*Sum[m[n]*r[n, t][[3]], {n, nn}];
ll[nn_] := tt[nn] - vv[nn];
rr[nn_] := 1/2 cd*Sum[Total[rdot[n, t]^2], {n, nn}];
eqLag[k_, nn_] :=
D[D[ll[nn], D[\[Theta][k, t], t]], t] -
D[ll[nn], \[Theta][k, t]] == -D[rr[nn], D[\[Theta][k, t], t]];
Block[{nn = 40,(*m=Function[1],l=Function[1],\[Zeta]=Function[1],*)
cd = 0.00005, g = 9.81},
eqs = Table[eqLag[n, nn], {n, nn}](*//Simplify*);
initialData =
Flatten[Table[{\[Theta][k, 0] == 0.05,
Derivative[0, 1][\[Theta]][k, 0] == 0}, {k, nn}]];
sols = NDSolveValue[
Flatten@{eqs, initialData,
WhenEvent[
0.02 (Sin[\[Theta][1, t]] + Sin[\[Theta][2, t]] +
Sin[\[Theta][3, t]] + Sin[\[Theta][4, t]] +
Sin[\[Theta][5, t]] + Sin[\[Theta][6, t]] +
Sin[\[Theta][7, t]] + Sin[\[Theta][8, t]] +
Sin[\[Theta][9, t]] + Sin[\[Theta][10, t]] +
Sin[\[Theta][11, t]] + Sin[\[Theta][12, t]] +
Sin[\[Theta][13, t]] + Sin[\[Theta][14, t]] +
Sin[\[Theta][15, t]] + Sin[\[Theta][16, t]] +
Sin[\[Theta][17, t]] + Sin[\[Theta][18, t]] +
Sin[\[Theta][19, t]] + Sin[\[Theta][20, t]] +
Sin[\[Theta][21, t]] + Sin[\[Theta][22, t]] +
Sin[\[Theta][23, t]] + Sin[\[Theta][24, t]] +
Sin[\[Theta][25, t]] + Sin[\[Theta][26, t]] +
Sin[\[Theta][27, t]] + Sin[\[Theta][28, t]] +
Sin[\[Theta][29, t]] + Sin[\[Theta][30, t]] +
Sin[\[Theta][31, t]] + Sin[\[Theta][32, t]] +
Sin[\[Theta][33, t]] + Sin[\[Theta][34, t]] +
Sin[\[Theta][35, t]] + Sin[\[Theta][36, t]] +
Sin[\[Theta][37, t]] + Sin[\[Theta][38, t]] +
Sin[\[Theta][39, t]] + Sin[\[Theta][40, t]])
<= 0,
Derivative[0, 1][\[Theta]][40, t] -> -0.8*
Derivative[0, 1][\[Theta]][40, t]]},
Table[\[Theta][k, t], {k, nn}], {t, 0, 20},
Method -> {"EquationSimplification" -> "Residual"}]]