Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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. >>
share|improve this question
1  
Then set $RecursionLimit to a bigger value, or even Infinity. You will have to be careful if you do $RecursionLimit = Infinity, though. –  J. M. May 18 '13 at 17:14
1  
@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. May 18 '13 at 17:26
    
Thanks. What is the correct command for blocking? Block[{$RecursionLimit = Infinity}, rez] in my example? –  Gasper May 18 '13 at 17:29
2  
@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
show 5 more comments

1 Answer

up vote 3 down vote accepted

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.

share|improve this answer
    
Works as charm, thanks! –  Gasper May 18 '13 at 18:22
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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