1

I know that this is a "beginner question", and, in fact, I am. I want to improve my code, as it takes really too much time to run. I have already read some other discussions, like here, but I am struggling in translating the simplifications with Table or Do into my example.

The For cycle I want to improve is the following:

zh = 0.4;
list = {};

For[i = 1, i < 10000, i++, 
    With[{zg = -0.6, dh = i10^-2}, 
         nsol = Block[{eps = $MachineEpsilon}, 
         NDSolve[{phi''[x] + 2phi'[x]/x + (2(zg + zhExp[-zh*x/dh])/x + 1)phi[x] == 0, phi[eps] == 1, phi'[eps] == -(zg + zh)}, phi, {x, eps, 20000}, 
                 WorkingPrecision->MachinePrecision, AccuracyGoal->15, PrecisionGoal->8, MaxSteps->Infinity]]];

    AppendTo[list, 1/Evaluate[(15000phi[15000])^2 + ((15000-Pi/2)*phi[15000-Pi/2])^2 /. nsol[[1]]]];]

Clearly, this code, written in this way, is highly inefficient. Also, I need to do more of these, with different values for zg inside With, and make some plots out of the lists.

Anyone that can help me with this noob question? Thanks a lot!

Lele
  • 428
  • 2
  • 10
  • Eliminate the AppendTo inside the loop and use Sow and Reap instead. – Daniel Huber Oct 21 '20 at 16:12
  • Can you provide me with an example? – Lele Oct 21 '20 at 16:17
  • You should take a look at this, if you haven't yet: https://mathematica.stackexchange.com/q/134609/12 The documentation of Sow and Reap has several relevant examples. – Szabolcs Oct 21 '20 at 16:43
  • 2
    Actually I don't see any need for Sow/Reap here. Just use Table. – Szabolcs Oct 21 '20 at 16:49
  • @Szabolcs I have already explicitly written that I had had a look into that link, but that I am not able to use Table in my case. – Lele Oct 21 '20 at 17:25
  • Here is a simple example: Reap[ Do[Sow[i], {i, 5}] ]. Note that Reap returns a list of 2 entries, the first one is the result of the expression inside Reap (here: Null), the second is a list of the elements Sow returned. – Daniel Huber Oct 21 '20 at 18:08
  • 1
    I don't see why Table wouldn't work—what exactly have you tried? Also, the code you posted here doesn't run (which is why I can't directly give you the Table translation). Please show a complete, minimal, working example. – Szabolcs Oct 21 '20 at 18:58
  • I mistakenly wrote i=0 as the starting point of the For loop in the question, but it is i=1, this is why it did not work for you. I've just corrected the error. Also, it appears that MMA complains if you define the differential equation before the For cycle, hence I put it dexplicitly inside. – Lele Oct 22 '20 at 12:44
  • Yes, With does literal replacement. a = b; With[{b=1}, a] does not work. You need Block for that, not With. – Szabolcs Oct 22 '20 at 12:51
  • You still didn't explain why this is not a straightforward translation to Table. – Szabolcs Oct 22 '20 at 15:44

1 Answers1

4

Here is an example that uses Table (well, ParallelTable) instead of For. I've also used ParametricNDSolveValue instead of With, mostly to simplify the Table.

phifunc = ParametricNDSolveValue[
  {\[Phi]''[x] + 2 \[Phi]'[x]/x + (2 (zg + zh Exp[-zh x/dh])/x + 1) \[Phi][x] == 0
   , \[Phi][$MachineEpsilon] == 1
   , \[Phi]'[$MachineEpsilon] == -(zg + zh)
   }
  , 1/((15000 \[Phi][15000])^2 + ((15000 - \[Pi]/2) \[Phi][
         15000 - \[Pi]/2])^2)
  , {x, $MachineEpsilon, 20000}
  , {dh \[Element] Reals, zg \[Element] Reals, zh \[Element] Reals}
]
list = ParallelTable[phifunc[i 10^-2, -0.6, 0.4], {i, 1, 10000}]

As far as timing goes, it still seems to take a very long time to run 10000 evaluations, so another answer might provide a faster method. I vaguely recall a way to functionalize calls to NDSolve in a way that stores calculations for faster repeated calls, but I can't find the link.

Josh Bishop
  • 741
  • 3
  • 8