3

I have been trying to plot those specific recursively defined functions p1 and q1 in Mathematica:

Δt = 10^(-1);
q1[j_] := q1[j] = q1[j - 1] + Δt*p1[j - 1];
p1[j_] := p1[j] = p1[j - 1] - Δt*Sin[q1[j - 1]];
q1[0] = π/12;
p1[0] = 0;
graph1q = ListPlot[Table[{j, q1[j]}, {j, 0, 20}], PlotStyle -> Green];
graph1p = ListPlot[Table[{j, p1[j]}, {j, 0, 20}], PlotStyle -> Green];

If you set jmax = 10, it's really fast, however the times raises too much when calculating more points. I have tried to put jmax = 50 and left for lunch. After an hour it didn't finish. Doing the same algorithm in Excel take less than a second!

Anybody got an idea why Mathematica has been taking so long to calculate the points?

Sektor
  • 3,320
  • 7
  • 27
  • 36
Geo
  • 461
  • 3
  • 12
  • 2
    It is quite fast if you do this with approximate numerical instead of exact evaluation. Try it with q1[0] = Pi/12.;p1[0] = 0.; – Daniel Lichtblau Nov 05 '14 at 16:41
  • @DanielLichtblau I did mention this, but was wrong to think of this as of a secondary reason, and to blame memoization. Changed my post now, apparently this is indeed the main and only reason for the slowness. – Leonid Shifrin Nov 05 '14 at 16:52

1 Answers1

7

The problem

The main problem seems that you have been using exact numbers for computations, since your delta is de facto an exact number. This led to generation of very large symbolic expressions, which also contributed to the slow-down. To solve this, one simply needs to use numerical values from the start, increasing precision if needed.

So, you can simply use N@Pi/2 and 0. as your starting values, and the code runs fast.

Iterative approach

One other possibility us to use an iterative approach, based on Nest (a similar case was discussed here). Initially, I posted this code since I blamed memoization. Apparently, the main problem is symbolic instead of numeric numbers in computations. Still, I'll keep this code just as an illustration.

pq = 
   NestList[
      {#1 + \[CapitalDelta]t*#2, #2 - \[CapitalDelta]t*Sin[#1]} & @@ # &, 
      N@{Pi/12, 0}, 
      20
   ];

The above generated both ps and qs, and takes no time. So now, 

graph1q = ListPlot[Table[{j, pq[[j + 1, 1]]}, {j, 0, 20}], PlotStyle -> Green]
graph1p = ListPlot[Table[{j, pq[[j + 1, 2]]}, {j, 0, 20}], PlotStyle -> Green]

If you go to larger values of j and use numerical approach, it may be a good idea to increase the precision at some point, because chances are that this procedure may lead to loss of precision at every iteration.

Community
  • 1
Leonid Shifrin
  • 114,335
  • 15
  • 329
  • 420
  • I think the primary reason for the slowness of Table is the large symbolic expressions generated. Using N early on is worth stressing. (You were still editing the answer while I was typing.) – Michael E2 Nov 05 '14 at 16:27
  • @MichaelE2 Will do, thanks. But have a look at my edit, I also fixed the recursive code. – Leonid Shifrin Nov 05 '14 at 16:28
  • 1
    Comment for the OP: It might be interesting to note that the memory used in the OP's definitions to store p1[j], q1[j] up to j == 20 given by ByteCount[{DownValues[p1], DownValues[q1]}] is about 0.5 GB; for Leonid's, it's less than 10K. – Michael E2 Nov 05 '14 at 16:42
  • @MichaelE2 Looks like my reasoning was partly wrong. Memoization seems to work also in the OP's case. Will change my answer. – Leonid Shifrin Nov 05 '14 at 16:46
  • Setting the condicitons numerical was such a simple problem that I would never find an answer. Thanks – Geo Nov 05 '14 at 18:40