1
$\begingroup$

My code is the following:

Clear["Global`*"]
Off[FindRoot::lstol]

Vn = (-G*Mn)/Sqrt[x^2 + y^2 + cn^2];
Vd = (-G*Md)/Sqrt[x^2 + y^2 + (s + h)^2];
Vb = (G*Mb)/(2*a)*(ArcSinh[(x - a)*(y^2 + c^2)^(-1/2)] - 
                   ArcSinh[(x + a)*(y^2 + c^2)^(-1/2)]);
Vh = (-G*Mh)/Sqrt[x^2 + y^2 + ch^2];
Vt = Vn + Vd + Vb + Vh;

G = 1; Mn = 400; cn = 0.25;
Md = 7000; s = 3; h = 0.175;
Mb = 3500; a = 2.5; c = 1;
Mh = 20000; ch = 20;
Ωb = 4.5;
E0 = -3245;

The initial data list

data0 = Table[{ϕ, L}, {ϕ, -π, π, 0.01}, {L, -110, 190, 1}];
ICs = Flatten[data0, 1];
norb = Length[ICs];
Print["N = ", norb]

Now I need to filter them using a For loop in order to obtain the valid initial conditions

data = {};
For[j = 1, j <= norb, j++, 
    ϕ0 = ICs[[j, 1]]; 
    L0 = ICs[[j, 2]];
    vt = L0/R0;

    x0 = R0*Cos[ϕ0];
    y0 = R0*Sin[ϕ0];
    px0 = -vt*Sin[ϕ0];
    py0 = vt*Cos[ϕ0];  
    H = 1/2*(px^2 + py^2) + Vt - Ωb*(x*py - y*px);
    H0 = H /. {x -> x0, y -> y0, px -> px0, py -> py0};

    sol = FindRoot[H0 == E0, {R0, 5}];
    R00 = R0 /. sol[[1]];
    E00 = H0 /. {R0 -> R00};
    Which[E00 == E0, AppendTo[data, {R00, ϕ0, L0}]]
];

The correct grid on the (ϕ, L)-plane is the following

enter image description here

As you can see, not all ICs are allowed since there are forbidden (white) regions on the grid. For all allowed initial conditions I know the triplet {ϕ, L, R}.

However the For loop requires a substantial amount of time in order to check all the initial conditions of the initial list.

Is there a way to speed up the procedure?

$\endgroup$

closed as unclear what you're asking by MarcoB, user9660, dr.blochwave, m_goldberg, Öskå Nov 9 '15 at 18:04

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Can't you just generate valid initial conditions to begin with? Generating a lot of pairs and then filtering them will have much worse complexity than a code that only generates valid initial conditions to begin with. $\endgroup$ – Marius Ladegård Meyer Nov 8 '15 at 15:38
  • $\begingroup$ @MariusLadegårdMeyer But I don't know beforehand which are the valid initial conditions. The filtering reveals the valid initial conditions. $\endgroup$ – Vaggelis_Z Nov 8 '15 at 15:43
  • $\begingroup$ All of these assignments are being done in the For loop which is adding substantial overhead. I would symbolically compute H0 then do a replacement like {phi0 -> #1, L0 -> #2} and apply the result to ICs $\endgroup$ – LLlAMnYP Nov 8 '15 at 16:32
  • $\begingroup$ @LLlAMnYP I don't exactly understand what you mean. Could you explain it further or post an answer? $\endgroup$ – Vaggelis_Z Nov 8 '15 at 16:36
  • 1
    $\begingroup$ Vaggelis_Z Could you explain in words what your code is trying to accomplish? I can't quite follow the logic, and I would still think that @Marius 's suggestion (generate the "right" numbers at the outset) might be a more effective way to go. Also, if you already have FORTRAN code that does what you want in 2 seconds, why don't you construct an interface to that, rather than trying to reinvent the wheel in Mathematica? $\endgroup$ – MarcoB Nov 8 '15 at 19:09
4
$\begingroup$

UPDATE 09.11.15

As @Marius has shown, it is quite easy to get an evenly spaced grid on ϕ, R0 for values of L0. First the equation H0 == E0 needs to be solved analytically for L0. Introduce initial conditions:

Vn = (-G*Mn)/Sqrt[x^2 + y^2 + cn^2];
Vd = (-G*Md)/Sqrt[x^2 + y^2 + (s + h)^2];
Vb = (G*Mb)/(2*a)*(ArcSinh[(x - a)*(y^2 + c^2)^(-1/2)] - 
                   ArcSinh[(x + a)*(y^2 + c^2)^(-1/2)]);
Vh = (-G*Mh)/Sqrt[x^2 + y^2 + ch^2];
Vt = Vn + Vd + Vb + Vh;

G = 1; Mn = 400; cn = 0.25;
Md = 7000; s = 3; h = 0.175;
Mb = 3500; a = 2.5; c = 1;
Mh = 20000; ch = 20;
Ωb = 4.5;
E0 = -3245;

vt = L0/R0;

x0 = R0*Cos[ϕ0];
y0 = R0*Sin[ϕ0];
px0 = -vt*Sin[ϕ0];
py0 = vt*Cos[ϕ0];  
H = 1/2*(px^2 + py^2) + Vt - Ωb*(x*py - y*px);
H0 = H /. {x -> x0, y -> y0, px -> px0, py -> py0};

Now solve:

sol = Solve[H0 == E0, L0] // FullSimplify;

In order to speed things up, let's get this into compiled form:

comp1 = With[{sol1 = L0 /. First@sol /. {ϕ0 -> #, R0 -> #2}}, 
   Compile[{ϕ0, R0}, Evaluate[sol1 & @@ {ϕ0, R0}]]];
comp2 = With[{sol2 = L0 /. Last@sol /. {ϕ0 -> #, R0 -> #2}}, 
   Compile[{ϕ0, R0}, Evaluate[sol2 & @@ {ϕ0, R0}]]];

Note the injection of parameters inside the compiled function (otherwise ϕ0, R0 in the compiled function are not the same symbols as \Phi[0], R0 in the solution. To avoid symbolic evaluation of compiled functions and also to handle the removable singularity at R0 == 0 augment this with a pattern matching only numeric arguments:

expr1[ϕ0_?NumericQ, R0_?NumericQ /; R0 > 0] := comp1[ϕ0, R0];
expr2[ϕ0_?NumericQ, R0_?NumericQ /; R0 > 0] := comp2[ϕ0, R0];
expr1[ϕ0_, 0.] = 0; expr2[ϕ0_, 0.] = 0;

Let's plot these solutions for some arbitrary ϕ0, for example 0.5:

ListPlot[{{expr1[.5, #], #} & /@ Range[0., 5., .01], {expr2[.5, #], #} & /@ Range[0., 5., .01]}]

L0(R0)

Blue is one branch of the solutions for L0, orange is another. And the points are evenly spaced on R0. But we want this the other way round. In fact, if we are thinking of R0 as a function of L0 and not the other way round, we see, that for several values of L0 there are two possible R0. Let's pick the lower branch.

ListPlot[{
  With[{tab = {expr1[.5, #], #} & /@ Range[0., 5., .01]},
   Pick[tab, {True}~
     Join~(Composition[TrueQ, (# < 0 &), First] /@ Differences[tab])]],
  With[{tab = {expr2[.5, #], #} & /@ Range[0.01, 5., .01]}, 
   Pick[tab, {True}~
     Join~(Composition[TrueQ, (# > 0 &), First] /@ Differences[tab])]]}]

Bottom branch

What I'm doing here, is picking only those elements of the first branch where for increasing R0 L0 decreases (as soon as it starts increasing we are in the upper branch of solutions) and similarly, increasing R0 -> increasing L0 in the second branch of L0(R0). The following expression will now simply yield a list of values for the lower branch in order of increasing L0:

With[{ϕ0 = #}, 
   With[{tab = {expr1[ϕ0, #], #} & /@ Range[0., 5., .01]},
     Pick[
       tab, {True}~
        Join~(Composition[TrueQ, (# < 0 &), First] /@ Differences[tab])] // 
      Reverse]~Join~
    With[{tab = {expr2[ϕ0, #], #} & /@ Range[0.01, 5., .01]},
     Pick[
      tab, {True}~
       Join~(Composition[TrueQ, (# > 0 &), First] /@ Differences[tab])]]] &[.5]

I've used an outer With to make this into a pure function which takes the value of ϕ0 as an argument and avoids conflict between inner and outer slots along the way.

Let's save this table for each value of ϕ0 from Range[-3.14, 3.14, 0.01]:

With[{ϕ0 = Round[#, 0.01]},
   r0tab[ϕ0] = 
    With[{tab = {expr1[ϕ0, #], #} & /@ Range[0., 5., .01]},
      Pick[tab, {True}~Join~(Composition[TrueQ, (# < 0 &), First] /@ Differences[tab])]
        // Reverse]
    ~Join~
    With[{tab = {expr2[ϕ0, #], #} & /@ Range[0.01, 5., .01]},
      Pick[tab, {True}~Join~(Composition[TrueQ, (# > 0 &), First] /@ Differences[tab])]]
   ] & /@ Range[-3.14, 3.14, 0.01];

This took roughly 10 seconds, already not too bad, (though I'm curious, how can this be done in a mere 2 seconds even with FORTRAN).

Now the table for each value of \[Phi0] needs to be resampled to get an even grid on L0:

With[{ϕ0 = Round[#, 0.01]}, 
     With[{min = Ceiling@First@First@r0tab[ϕ0], 
       max = Floor@First@Last@r0tab[ϕ0], 
       func = Interpolation@r0tab[ϕ0]},
      Table[{ϕ0, L0, func[L0]}, {L0, min, max, 3.}]]] & /@ 
   Range[-3.14, 3.14, 0.05] // Flatten[#, 1] & // ListPointPlot3D

Result

Adjust the step size in the Table function to change the sampling density for L0, but keep in mind that it is also dependent on the sampling density for R0 in the previous commands. Adjust the step size in Range[-3.14, 3.14, 0.05] for the sampling density in Phi but the values must be the same as those saved in r0tab.

$\endgroup$
  • $\begingroup$ When I evaluate the two With I get the following error: Syntax::sntxf: "TrueQ@" cannot be followed by "*(#<0&)@*First/@Differences[tab]". $\endgroup$ – Vaggelis_Z Nov 9 '15 at 9:46
  • $\begingroup$ @Vaggelis_Z what version are you using? The fix is probably a pair of parentheses that need inserting, but in older versions of MMA the shorthand operator for Composition was not present. $\endgroup$ – LLlAMnYP Nov 9 '15 at 9:47
  • $\begingroup$ I use v9.0 of MMA. $\endgroup$ – Vaggelis_Z Nov 9 '15 at 9:48
  • $\begingroup$ I see. The solution is to replace TrueQ@*(...)@*First with Composition[TrueQ, (# < 0 &), First]. Will update in a second. $\endgroup$ – LLlAMnYP Nov 9 '15 at 9:51
  • $\begingroup$ @Vaggelis_Z fixed to work with v9 syntax (infix notation for Composition was introduced in v10). $\endgroup$ – LLlAMnYP Nov 9 '15 at 9:53
4
$\begingroup$

You can solve the problem symbolically for either L0 or R0 once and for all, and then create a table using these solutions. It's most easy to solve for L0 since it does not appear in square roots or trigonometric functions.

Clear[L0,R0];
vt = L0/R0;

x0 = R0*Cos[ϕ0];
y0 = R0*Sin[ϕ0];
px0 = -vt*Sin[ϕ0];
py0 = vt*Cos[ϕ0];  
H = 1/2*(px^2 + py^2) + Vt - Ωb*(x*py - y*px);
H0 = H /. {x -> x0, y -> y0, px -> px0, py -> py0};
sol = Solve[H0 == E0, L0];

Since this is a quadratic equation in L0 you will get two solutions. Choosing the first one for no particular reason, you can then generate the table of allowed values using

Table[{R0, ϕ0, L0 /. sol[[1]]}, {R0, 0.1, 1, 1/621}, {ϕ0, -Pi, Pi, 0.01}];

You will have to choose physically reasonable values for R0 yourself. You can even inspect the solution and create a compiled version for a possible speed boost...

$\endgroup$
  • $\begingroup$ NO. R0 should be computed from a given pair of phi0 and L0. I don't know how to explain in detail but this has to do with maps and flows in nonlinear dynamics. It is the transformations form (\phi, L, R) to (x,y,px,py). Your solution gives completely wrong output. $\endgroup$ – Vaggelis_Z Nov 8 '15 at 19:45
  • $\begingroup$ @Vaggelis_Z "I don't know how to explain in detail" doesn't help much to others trying to help you. Perhaps you should try harder to "explain in detail" $\endgroup$ – Dr. belisarius Nov 8 '15 at 20:18
  • $\begingroup$ @belisariushassettled H0 is a function with three parameters (phi0,L0,R0). The first two (phi0,L0) are given at the beginning of the For loop. Then we use FindRoot to obtain the R0 value for a given par of (\phi0,L0). If R0 does not exist FindRoot gives a wrong answer, so at the end I use Which to obtain only the triplets (R0,phi0,L0) which give the desired energy level E0. $\endgroup$ – Vaggelis_Z Nov 8 '15 at 20:26
  • $\begingroup$ @Vaggelis_Z That's a description of what you've done, but it doesn't explain what you think is wrong with Marius' answer $\endgroup$ – Dr. belisarius Nov 8 '15 at 20:32
  • $\begingroup$ @Vaggelis_Z so you have some implicitly defined surface in the space of {\[Phi], L, R}, such that H[\[Phi], L, R] == E0. What difference does it make, which order to sample the grid in? $\endgroup$ – LLlAMnYP Nov 8 '15 at 20:32

Not the answer you're looking for? Browse other questions tagged or ask your own question.