Random number in system of equations, solving using NDSolve & WhenEvent

Question

I need to analyse the effect of random forcing on the system of coupled equations. For example, I have shown the equations below.

$$\begin{pmatrix} x1''\\x2''\\x3''\\x4''\\x5'' \end{pmatrix} + M_{(5 \times 5)} \begin{pmatrix} x1\\x2\\x3\\x4\\x5 \end{pmatrix} + V2_{(5 \times 1)} xb = \begin{Bmatrix} f1\\ f2\\ f3\\ f4\\ f5 \end{Bmatrix}, \\ \\ \text{and } \;\;\; xb=f(x5,t)$$

Here, F={f1,f2,f3,f4,f5} is the random force vector. I have used WhenEvent which induces random number at ceratain dt interval for random force. For two simplified coupled equations, I have mentioned the code below according to this answer.

tf=100; dt=0.001;  
sol2=NDSolve[{x''[t]+x'[t]+2*x[t]-y[t]==r1[t],y''[t]+y'[t]+2*y[t]-x[t]==r2[t],WhenEvent[Mod[t,dt]==0,
{r1[t],r2[t]}->RandomReal[{-0.01,0.01},2]],{x[0]==y[0]==x'[0]==y'[0]==r1[0]==r2[0]==0}},{x,y},{t,0,tf},
DiscreteVariables->{r1,r2},Method->{"StiffnessSwitching",Method->{"ExplicitRungeKutta",Automatic}},
MaxSteps->Infinity,AccuracyGoal->7,PrecisionGoal->7][[1]];

Above code is working, but, I need to carryout the analysis for a large system of equations, say more than 20. So, for that I have written down the following code.

tf = 10; dt = 0.01;  N1 = 10;  \Omega=2;
Y1[t_] := Table[Subscript[x, i][t], {i, 2, N1}]; F[t_] := Table[Subscript[f, i][t], {i, 2, N1}];   
a = ConstantArray[0, N1 - 1]; a[[N1-1]]=1;
  xb[t] = Subscript[x, N1][t] \Cos[\Omega t];
eqns = Thread[D[Y1[t], t, t] + D[Y1[t], t] + Y1[t] ==a*xb[t]+ F[t]];   
X0 = ConstantArray[0, N1 - 1]; Xd0 = ConstantArray[1, N1 - 1]; 
sol = NDSolve[{eqns, WhenEvent[Mod[t, dt] == 0, F[t] -> RandomReal[{-0.01, 0.01}, 2]], 
    Y1[0] == F[0] == X0, Y1'[0] == Xd0}, {Y1[t], F[t]}, {t, 0, tf}, DiscreteVariables -> F, 
   Method -> {"StiffnessSwitching",  Method -> {"ExplicitRungeKutta", Automatic}}, 
   MaxSteps -> Infinity, AccuracyGoal -> 7, PrecisionGoal -> 7][[1]]

Plot[Evaluate[{ Subscript[x, N1][t]} /. First@sol], {t, 0, tf}, 
 PlotRange -> Full] 

It is working well in the absence of random number generator and WhenEvent approach. While consideration of whenevent for random number generation, it gives the following error.

NDSolve::underdet: There are more dependent variables, {Subscript[f, 2][t],Subscript[f, 3][t],Subscript[x, 2][t],Subscript[x, 3][t],{Subscript[f, 2][t],Subscript[f, 3][t]}}, than equations, so the system is underdetermined.

Can someone please help me, to solve the issue?

Answer 1

Updated to Accommodate Change in $a$ Definition

The OP changed the definition of the symbol, $a$ from my original answer. The code is essentially the same as before:

tf = 10; dt = 0.1; N1 = 10; Ω = 2;
a = ConstantArray[0, N1 - 1]; a[[N1 - 1]] = 1;
zeros = ConstantArray[0, N1 - 1];
Y1[t_] := Table[Subscript[x, i][t], {i, 2, N1}]; 
F[t_] := Table[Subscript[f, i][t], {i, 2, N1}];
dvbls = Head[#] & /@ F[t];
xb[t] = Subscript[x, N1][t] Cos[Ω t];
eqns = Thread[D[Y1[t], t, t] + D[Y1[t], t] + Y1[t] == a*xb[t] + F[t]];
ics = Thread[Y1[0] == zeros]~Join~Thread[Y1'[0] == zeros + 1];
istates = Thread[F[0] == zeros];
events = WhenEvent[Mod[t, dt] == 0, 
     Subscript[f, #][t] -> RandomReal[{-0.2, 0.2}]] & /@ Range[2, N1];
(* Join System of Equations to be Solved *)
system = Join[eqns, ics, istates, events];
sol = NDSolve[system, Join[Y1[t], F[t]], {t, 0, tf}, 
   DiscreteVariables -> dvbls, 
   Method -> {"StiffnessSwitching", 
     Method -> {"ExplicitRungeKutta", Automatic}}, 
   MaxSteps -> Infinity, AccuracyGoal -> 7, PrecisionGoal -> 7];
Plot[Evaluate[Table[Subscript[x, i][t], {i, 2, N1}] /. sol], {t, 0, 
  tf}, PlotRange -> Full]
Plot[Evaluate[Table[Subscript[f, i][t], {i, 2, N1}] /. sol], {t, 0, 
  tf}, PlotRange -> Full]

Original Answer

I am not sure if your equations sense (e.g., in your case, the Cos term cancels because $a$ is a zero array) or just a simplified version, but you will need to Thread your initial conditions and states like you did for eqns. Also, you will need to create a list of events. Below is an example (I adjusted some parameters for visualization):

tf = 10; dt = 0.1; N1 = 10; Ω = 2;
a = ConstantArray[0, N1 - 1];
zeros = ConstantArray[0, N1 - 1];
Y1[t_] := Table[Subscript[x, i][t], {i, 2, N1}]; 
F[t_] := Table[Subscript[f, i][t], {i, 2, N1}];
dvbls = Head[#] & /@ F[t];
xb[t] = Subscript[x, N1][t] Cos[Ω t];
eqns = Thread[D[Y1[t], t, t] + D[Y1[t], t] + Y1[t] == a*xb[t] + F[t]];
ics = Thread[Y1[0] == zeros]~Join~Thread[Y1'[0] == zeros + 1];
istates = Thread[F[0] == zeros];
events = WhenEvent[Mod[t, dt] == 0, 
     Subscript[f, #][t] -> RandomReal[{-0.2, 0.2}]] & /@ Range[2, N1];
(* Join System of Equations to be Solved *)
system = Join[eqns, ics, istates, events];
sol = NDSolve[system, Join[Y1[t], F[t]], {t, 0, tf}, 
   DiscreteVariables -> dvbls, 
   Method -> {"StiffnessSwitching", 
     Method -> {"ExplicitRungeKutta", Automatic}}, 
   MaxSteps -> Infinity, AccuracyGoal -> 7, PrecisionGoal -> 7];
Plot[Evaluate[Table[Subscript[x, i][t], {i, 2, N1}] /. sol], {t, 0, 
  tf}, PlotRange -> Full]
Plot[Evaluate[Table[Subscript[f, i][t], {i, 2, N1}] /. sol], {t, 0, 
  tf}, PlotRange -> Full]