Message Boards Message Boards

Numerical instability in simulation

Posted 1 month ago

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"}]]
POSTED BY: Ee Kin Chan
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract