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[0][[i]] == p0[[i]],
rp'[0][[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[0][[i]] == p0[[i]],
rp'[0][[i]] == vel[0][[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