WhenEvent, Storing Sow data, Loops

Question

My problem is essentially to do with manipulating and storing Sow data in a nested loop structure.

I have a set of differential equations describing time-evolution of a bunch of variables. Eventually I want the initial conditions i0, j0, k0 to be parameters. I have a range for each of these that I want to scan over. Eventually, I want to find out the t value at which point i[t]==1.

So far I can set up an NDSolveValue taking parameteres, and Reap for a specific set of parameters.

test[i0_, j0_, k0_] := NDSolveValue[{coupled nonlinear first order ODEs,
i[0] == i0, j[0] == j0, k[0] == k0, WhenEvent[i[t] == 1, Sow@{i0, j0, k0};
Sow[t]]}, {i, j, k}, {t, 0, 15}]; 

test[0.4, 0.001, 0.03]; // Reap

The output makes sense, giving the t-value as 11.2293 for this case which seems correct:

{Null, {{{0.4, 0.001, 0.03}, 11.2293}}}

I have been struggling with scanning across a specified i0, j0, k0 range. Here is something that kind of works (StopIntegration is to avoid a pole that happens at some late t):

Table[Table[Table[
NDSolveValue[{the ODEs, i[0] == i0, j[0] == j0, k[0] == k0, 
WhenEvent[i[t] == 1, Print[t,i0,j0,k0]; 
Sow[t, {i0,j0,k0}]; "StopIntegration"]}, {i, j, k}, {t, 0, 25}], 
{i0,0.3, 0.4, 0.1}], {j0, 0.0001, 0.005, 0.001}], {k0, 0.1, 0.5, 0.1}]

And this gives me:

22.1212....0.3....0.0001....0.1

11.2014....0.4....0.0001....0.1

6.79366....0.5....0.0001....0.1

22.1209....0.3....0.0021....0.1

etc. (ellipses added by me)

for each combination of values, which also seem correct. However, I am struggling with saving these results in a proper data (.dat) file. I want to ultimately make a table with columns i0, j0, j0, t which can be exported as a .dat file. Ideally, I would also like to append to this table the columns f(i0,j0,k0), g(i0,j0,k0), which are certain functions of the values. How do I create a professional data table out of this data (I've been using 'Print' because I don't know how to use Sow/Reap data properly)?

Edit: If I add Reap@ to the start of the expression it gives me a list of t-values for i0=i1 then for i0=i2 then i0=i3, that is, for each constant i from the i range, with the others varying; followed by a list of t-values for j0=j1,....

At this stage my query is mainly to do with formatting the data, i.e. generating a neat exportable table. How do I convert the Reap output into a data table?

Any help would be appreciated!

Answer 1

There are always at least a dozen different ways of doing anything in Mathematica, and that doesn't even count turning function names into punctuation. Here are three ways.

Assume that stuff does your calculations, ends up assigning a value to t and that multiple steps are separated by semicolons.

Table with multiple iterators will create nested lists, but Flatten with an appropriate 2nd argument will reduce that to a simple matrix.

nBy4Matrix = Flatten[Table[stuff; {t, i, j, k}, {i,1,3}, {j,1,3}, {k,1,3}], 2]

Sow and Reap will let you accumulate results from any sequence of calculations.

nBy4Matrix = Reap[Do[stuff; Sow[{t, i, j, k}], {i,1,3}, {j,1,3}, {k,1,3}]][[2,1]]

Map will let you perform a calculation on each item from a list, each of which happens to be a triple in your case, and return a list of the results to you.

nBy4Matrix = Map[({i,j,k}=#; stuff; {t, i, j, k})&, {{1,1,1}, {2,3,1}, {3,1,3}, {2,1,1}}]

Pick one way that you will be able to remember in the future and be able to use without making any mistakes. That is the real goal of all this.