2

I noticed that for large matrix, the rule /. is very slow. In the example below, the combination of memorization and the replacement "/.sol[n-1]" inside the definition of sol[n] turns the code extremely slow. Just to calculate sol[1], it takes around 20 sec. So for sol[90] it will take 30 min! Is there any way to use a faster "replacement" technique than /. ?

ClearAll["Global`*"]
sp = 100; ϵ = 0.2; β = 0.03; ρ = 0.5; α = 0.8; γ = 0.8;
V0 = 100; 
J = 50; K = 50;
T = 100; Δt = 0.5;
Xmax = 300; Xmin = -200;
Vmax = 300; Vmin = 0;
ΔX = (Xmax - Xmin)/K;
ΔV = (Vmax - Vmin)/J;
F0[Xk_, Visc_] := Visc - sp;
F[k_, 0, n_] := 0;
F[k_, J, n_] := 0;
F[0, j_, n_] := 0;
F[K, j_, n_] := 0;
a[j_] := 0.5 β j Δt;
b[j_] := 1 + ( ϵ^2 j^2 + β) Δt;
c[j_] := -0.5 β j Δt;
d[j_] := -(0.25 ϵ ρ j/ ΔX) Δt;
sol[0] = Flatten[
  Table[F[k, j, 0] -> F0[Xmin + k ΔX, j ΔV], {k, 1, K - 1}, {j, 1, J - 1}]];
sol[n_] := 
  sol[n] = Module[{vars, eqns}, 
    vars = Flatten[Table[F[k, j, n], {k, 1, K - 1}, {j, 1, J - 1}]]; 
    (*THIS IS THE   PROBLEMATIC STEP!*)
    eqns = Flatten[
      Table[a[j] F[k, j - 1, n] + b[j] F[k, j, n] + c[j] F[k, j + 1, n] + 
        d[j] (F[k + 1, j + 1, n] + F[k - 1, j - 1, n] - 
        F[k - 1, j + 1, n] - F[k + 1, j - 1, n]) == 
        F[k, j, n - 1], {k, 1, K - 1}, {j, 1, J - 1}]] /. sol[n - 1];
    temp = Solve[eqns, vars][[1]];
    FF = ArrayReshape[Values@temp, {K - 1, J - 1}];
    Fc[k_, j_] := FF[[k, j]];
    Flatten[Table[F[k, j, n] -> Fc[k , j ], {k, 1, K - 1}, {j, 1, J - 1}]]]
sol[1] // Timing
sol[900] // Timing (*This will take forever!*)
m_goldberg
  • 107,779
  • 16
  • 103
  • 257
duartedias
  • 21
  • 2
  • Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Feb 21 '15 at 04:36
  • Dispatch is your friend... that said, the generation can be done more directly/efficently also. – ciao Feb 21 '15 at 06:05
  • Thanks, rasher! The difference on speed is amazing. I wonder why isn't dispatch the default built in /. Do you know how can access the content of a Dispatch output without /. ? I've notice that the output is composed by "Length" and "Rules", but I would like to access "Rules" without creating another Table to use /. – duartedias Feb 21 '15 at 12:01
  • (1) It costs time to build the dispatch table. (2) Related: http://mathematica.stackexchange.com/questions/31613/how-to-speed-up-calculations-with-large-number-of-replacement-rules – Michael E2 Feb 21 '15 at 14:32

0 Answers0