# Update initial conditions recursively/sequentially within ParallelTable using NDSolve

I am modelling 3D motion of particles inside a force-field using second-order coupled PDEs, which I solve from time t=t0 to t=t1 with initial conditions x0,y0,z0. I have a total number of "particles" n with subscript i and for each of these particles I model j collisions.

It is relatively straightforward to obtain the position vector r[t1]=={x[t1],y[t1],z[t1]} and the velocity vector r'[t1] for a single particle p_i by solving the PDEs once, but what I want to do is keep track of these outputs and use them as inputs in the subsequent calculation. That is to say, I want to then solve my PDEs from t=t1 to t=t2 using initial conditions r[t1] and r'[t1], and so on until I'm done . My code so far is:

f[nn_] :=
Module[{n = nn},
temp = Quiet@
ParallelTable[(NDSolve[{x''[t] == -x'[t] - x[t],
y''[t] == -y'[t] - y[t], z''[t] == -z'[t] - z[t], x[0] == 1,
x'[0] == 1, y[0] == 1, y'[0] == 1, z[0] == 1,
z'[0] == 1}, {x, y, z}, {t, tjs[[i, j]], tjs[[i, j + 1]]}]) //
First, {i, n}, {j, Length[tjs[[i]]] - 1}]]
nn = 20;
et = 2.1 10^7;
p = 1;
pdft[t_, p_] = et p 10^-6 E^(-(et p 10^-6) t);
trndm[p_] :=
RandomVariate[ProbabilityDistribution[pdft[t, p], {t, 0, Infinity}]]
tj := Prepend[(sum = 0;
Reap[While[sum < 0.05, sum = Sow[sum + trndm[p]]]][[2, 1]]), 0]
tjs = Table[tj, nn];
Table[(tjs[[i, Length[tjs[[i]]]]] = 0.05), {i, nn}];
f[nn]


So, each time I run the module f[nn] it properly solves j collisions for each particle i, but it does so using the same initial conditions every time.

I've tried to provide a MWE so if something is missing or needs clarification please let me know.

• Your code does not work as posted (NDSolve complains that the system is underdetermined). 1) copy and past the code you posted, see if you can fix it so it runs; 2) for now, do away with the complication of parallel execution and concentrated on getting the code to run. – MarcoB Mar 6 at 22:15
• I intended to provide the algorithm without going into details, but I'll add a version of the code that actually executes now. – Milad Dakka Mar 7 at 1:45
• @MarcoB, it is done! – Milad Dakka Mar 7 at 1:52
• @MiladDakka What do you want to get in the result? – Alex Trounev Mar 7 at 11:56
• @AlexTrounev, I have answered my own question! What I personally want is to perform a statistical analysis of the result of millions of collisions for a particular research question I am investigating. For this purpose, I needed to "follow" the trajectory of atoms for a certain period of time, a trajectory that includes several "collisions", which I modelled as velocity changes. See the answer below for more info. – Milad Dakka Mar 19 at 11:02

I have found an answer to my own question using advice from a Mathematica guru and I had to update this question as a result.

So, it turns out the WhenEvent function in Mathematica is the most wonderful answer to this question, as it allows one to specify an event (in my case a collision occurring at pre-determined times) and what one wants to happen as a result (in my case, the atom experiences a velocity change due to the collision) neatly within the NDSolve function.

So my example code, which includes "collisions" at times t=0.01, t=0.05 and t=0.08, now looks something like this:

NDSolveValue[
{
x''[t] == -x'[t] - x[t],
y''[t] == -y'[t] - y[t],
z''[t] == -z'[t] - z[t],
x[0] == 1, x'[0] == 1,
y[0] == 1, y'[0] == 1,
z[0] == 1, z'[0] == 1,
WhenEvent[
{t==0.01,t==0.05,t==0.08},
{x'[t]->RandomReal[{0,1}],
y'[t]->RandomReal[{0,1}],
z'[t]->RandomReal[{0,1}]}
]
},
{x,y,z},
{t,0,0.1}
]


Here is a very helpful link to Events and Discontinuities in DEs, with examples for many different cases that could be of use to anyone who comes across this kind of issue.

• Are you sure this code describes a collision of molecules? – Alex Trounev Mar 19 at 11:36
• Hi Alex, this code isn't meant to represent real collisions, only the structure of the code that could do so. The forces in this case are not very realistic and the collisions are assumed to result in random velocities. But yes, the essential idea is useful for describing collisions (or events). – Milad Dakka Mar 24 at 10:06
• There are several posts on this topic, for example mathematica.stackexchange.com/questions/38687/… – Alex Trounev Mar 24 at 11:22