2
$\begingroup$

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. >>
$\endgroup$
  • 1
    $\begingroup$ Then set $RecursionLimit to a bigger value, or even Infinity. You will have to be careful if you do $RecursionLimit = Infinity, though. $\endgroup$ – J. M. will be back soon May 18 '13 at 17:14
  • 1
    $\begingroup$ @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. $\endgroup$ – Leonid Shifrin May 18 '13 at 17:22
  • $\begingroup$ @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. $\endgroup$ – J. M. will be back soon May 18 '13 at 17:26
  • $\begingroup$ Thanks. What is the correct command for blocking? Block[{$RecursionLimit = Infinity}, rez] in my example? $\endgroup$ – energyMax May 18 '13 at 17:29
  • 2
    $\begingroup$ @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). $\endgroup$ – Leonid Shifrin May 18 '13 at 17:34
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Works as charm, thanks! $\endgroup$ – energyMax May 18 '13 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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