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
