# NDSolve with arrays and Tables of Equations and WhenEvent

I've been toying around with NDSolve for a while, and read through the website. By doing so I discovered that I could use it for vectors and arrays with much pleasure. So I wanted to write a simple simulator with multiple balls or particles interacting with Lennard Jones potentials. When I write this way I have no problem

partsim[p0_, v0_, mass_] :=
Module[{rp},
rp[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
rp[t] /. NDSolve[
Join[
Table[
{
(*Initial COnditions*)
rp[[i]] == p0[[i]],
rp'[[i]]  == v0[[i]],

(*Newton Equations*)
rp''[t][[i]] == (1/mass[[i]])*
Sum[If[j != i, ljForce[i, j, rp[t]], 0], {j, 1,
nParticles}]
(*Third newton law in the future*)

}
, {i, nParticles}](*End of Table*)
](*End of Join*)
, rp[t]
, {t, 0, ftime}] // First
](*End of Function*);


The real problem arises when I'm trying to also add a WhenEvent for each of the particles present in the system. Since it's 3d I would like them to stay inside a bounding box and invert the sign of their velocity every time they touch the walls or at least come withing tolerance.

partsim[p0_, v0_, mass_] :=
Module[{rp, vel},
rp[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
vel[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
rp[t] /. NDSolve[
Join[
Table[
{
(*Initial COnditions*)
rp[[i]] == p0[[i]],
rp'[[i]] == vel[[i]] == v0[[i]],
(*Newton Equations*)
rp'[t][[i]] == vel[t][[i]],
rp''[t][[i]] == (1/mass[[i]])*
Sum[If[j != i, pairForce[i, j, rp[t]], 0], {j, 1,
nParticles}]

}
, {i, nParticles}](*End of Table*)

(*Discrete Events*)
(*Bouncing on walls*)

, Table[{
WhenEvent[
{Norm[rp[t][[u]] - wallcoord] < .01  ,
Norm[rp[t][[u]]] <  .01}
, vel[t][[u]] -> -vel[t][[u]]
](*End of WhenEvent*)
}, {u, nParticles}]

](*End of Join*)

, rp[t]
, {t, 0, ftime}, DiscreteVariables -> {vel}] // First
](*End of Function*);


The error that appears seems to be related to the fact that u remains unevaluated, as attempting to Table it outside of the program will reveal. I have added DiscreteVariables -> {vel} because I thought that since we are subjecting velocity to a sudden change it would be necessary. Is there a way of adding a list of working WhenEvent, one for each particle, to the system of equations? Many thanks

• Please post a minimally working example, with everything required to use it - any data, parameters, arguments to use, etc.. As it stands, it is an exercise in mind-reading. – ciao Mar 6 '14 at 11:50

In lieu of a working example in the current OP to work with, yes, this can be done. An example of two bouncing balls with differing velocity retention on bounce:

eventtab =
Table[With[{n = n},
WhenEvent[y[n][t] == 0, y[n]'[t] -> (-0.95 + n/100) y[n]'[t]]], {n,
2}]

NDSolve[{y''[t] == -9.81, y == 5, y' == 0,
y''[t] == -4, y == 3, y' == 0, eventtab}, {y,
y}, {t, 0, 10}]

Plot[Evaluate[{y[t], y[t]} /. %], {t, 0, 10}] • This was a very useful answer, I did not think about using With. I knew that WhenEvent had the Hold attribute, but I did not know I could bypass it. – BurgerAutomata Mar 6 '14 at 13:32