Recursion depth exceeded

Question

Below is my code for numerical solving of PDE with Crank Nicolson scheme.

 ClearAll[L, M, \[CapitalDelta]t, \[CapitalDelta]x]
    L = Pi // N; M = 25;
    \[CapitalDelta]x = L/M;
    \[CapitalDelta]t = \[CapitalDelta]x^2; k = 1;
    U0[x_] = x (Pi - x);
    U[0, n_] := 0.;
    U[M, n_] := 0.;
    sol[0] = Table[U[j, 0] -> U0[j \[CapitalDelta]x], {j, 1, M - 1}];
    sol[n_] := 
     Module[{vars, eqns}, vars = Table[U[j, n], {j, 1, M - 1}]; 
      eqns = Table[
         U[j, n] - U[j, n - 1] == 
          1/2 ((k \[CapitalDelta]t )/\[CapitalDelta]x^2 (U[j + 1, n] - 
               2 U[j, n] + U[j - 1, n] + U[j + 1, n - 1] - 2 U[j, n - 1] +
                U[j - 1, n - 1])), {j, 1, M - 1}] /. sol[n - 1];
      Solve[eqns, vars][[1]]
     ]

The problem is, that I have to eveluate value of scheme for (Pi/4,2). That means i have to put 127 for n in code below. The problem is, that the program doesn't want to eveluate for such big numbers.

rez = Table[
  Table[{\[CapitalDelta]x j, \[CapitalDelta]t n, U[j, n]}, {j, 0, 
     M}] /. sol[n], {n, 127}]

But if you go bigger for n, like 100, 400, and more it gives:

$RecursionLimit::reclim: Recursion depth of 256 exceeded. >>

Then set $RecursionLimit to a bigger value, or even Infinity. You will have to be careful if you do $RecursionLimit = Infinity, though. — J. M.'s missing motivation, May 18 '13 at 17:14
@J.M. Actually, I think that setting $RecursionLimit = Infinity is never appropriate. I realize that I may be biased by my software development (rather than problem - solving) perspective, but it is all too easy to lose unsaved data by using this setting. — Leonid Shifrin, May 18 '13 at 17:22
@Leonid, of course, I usually do that as Block[{$RecursionLimit = Infinity}, (* stuff *)] instead of setting it globally. I did say that OP has to be careful; in particular, he should be sure that the recursion he has is guaranteed to terminate, but usually people are unable or unwilling to do such proofs. — J. M.'s missing motivation, May 18 '13 at 17:26
Thanks. What is the correct command for blocking? Block[{$RecursionLimit = Infinity}, rez] in my example? — energyMax, May 18 '13 at 17:29
@J.M. Sure, I just wanted to reinforce your comment. Re- Block - this is better, but IMO not enough. Actually, presence of absence of Block is irrelevant for my main argument. My point is that even with Block, setting $RecursionLimit = Infinity can easily overflow the stack which leads to a kernel crash. Happened to me more than once :-). — Leonid Shifrin, May 18 '13 at 17:30
@J.M. It also depends on how recursion is used. If one uses it for algorithm which is deeply recursive (high recursion depth), then the right way to do this in Mathematica is to perform some form of tail call optimization. Note that tail call - optimized functions in Mathematica require modifications of $IterationLimit rather than $RecursionLimit, and the former are safe (it is safe to set $IterationLimit to Infinity). — Leonid Shifrin, May 18 '13 at 17:34
@Leonid, I had forgotten about tail-call optimization. You are most certainly right. ;) Now, I wonder if OP's function can be reimplemented in that manner... then again, if memory serves, there's a nonrecursive implementation of Crank-Nicolson; I wonder why this isn't the route taken. — J. M.'s missing motivation, May 18 '13 at 17:37
@Gasper, if you're going to gamble with that setting, you'll want to do it in the definition of sol[]; that is, sol[n_Integer] := Block[{$RecursionLimit = Infinity}, (* stuff *)] — J. M.'s missing motivation, May 18 '13 at 17:42
I put n=1200 and after 10min it's still not evaluated. Any reccomendations for fastening the calculation process? — energyMax, May 18 '13 at 17:42
@J.M. In this case, conversion to iterative version is made trivial via memoization, see my answer. — Leonid Shifrin, May 18 '13 at 18:04

score 3 · Accepted Answer · edited Apr 13 '17 at 12:55

Ok, your case is the one which can be easily helped by memoization. Basically, it will serve to compute the set of values in sol in reverse order, de facto making it just a convenient device to turn your algorithm into an iterative one without changing it much.

All you have to do is to change your definition of sol as:

sol[n_] := sol[n] = Module[your-old-code]

and then, if you want say n=1200, first execute

Scan[sol,Range[1200]]

This will automatically add definitions for specific values of sol for all n from 1 to 1200, and then you simply call

sol[1200]

(* 
   {U[1,1200]->1.92724*10^-9,U[2,1200]->3.82408*10^-9,
    U[3,1200]->5.66062*10^-9,U[4,1200]->7.40789*10^-9,
    U[5,1200]->9.03832*10^-9,U[6,1200]->1.05262*10^-8,
    U[7,1200]->1.18481*10^-8,U[8,1200]->1.29832*10^-8,
    U[9,1200]->1.39134*10^-8,U[10,1200]->1.46243*10^-8,
    U[11,1200]->1.51045*10^-8,U[12,1200]->1.53466*10^-8,
    U[13,1200]->1.53466*10^-8,U[14,1200]->1.51045*10^-8,
    U[15,1200]->1.46243*10^-8,U[16,1200]->1.39134*10^-8,
    U[17,1200]->1.29832*10^-8,U[18,1200]->1.18481*10^-8,
    U[19,1200]->1.05262*10^-8,U[20,1200]->9.03832*10^-9,
    U[21,1200]->7.40789*10^-9,U[22,1200]->5.66062*10^-9,
    U[23,1200]->3.82408*10^-9,U[24,1200]->1.92724*10^-9}
*)

and you get this under a couple of seconds.

Recursion depth exceeded

1 Answers1