How to speed up ReplaceAll in Table?

Question

I am writing a code which solve a list of equations and then create a Table with resulst of equation solving. In applications where I use this code, I create a bunch of equations and my computations can be very long (more than 3-4 days).

I investigate that Table creation with ReplaceAll function is my narrow point, which calculate much more longer than other parts. Here is an example:

  In[171]:= ClearAll["Global`*"]

In[172]:= t1 = AbsoluteTime[];

In[173]:= m = 150;
p[1, k_] = 0;
p[m, k_] = 1;
p[j_, 1] = 0;
p[j_, m] = 0;
dx = dy = 0.01;

In[179]:= 
eqns = Flatten[
   Table[1/dx^2 (p[j + 1, k] - 2 p[j, k] + p[j - 1, k]) + 
      1/dy^2 (p[j, k + 1] - 2 p[j, k] + p[j, k - 1]) == 0, {j, 2, 
     m - 1}, {k, 2, m - 1}]];

In[180]:= vars = Flatten[Table[p[j, k], {j, 2, m - 1}, {k, 2, m - 1}]];

In[181]:= sol = Solve[eqns, vars][[1]];

In[182]:= AbsoluteTime[] - t1

Out[182]= 0.8736015

In[183]:= 
pTable = Table[{j dx, k dy, p[j, k] /. sol}, {j, 1, m}, {k, 1, m}];

In[184]:= AbsoluteTime[] - t1

Out[184]= 38.4112677

Is there any way to improve this operation? Last time when I just move "ReplaceAll" command near p[j, k] it speed up calculation around 40%.

Maybe there is another hint, which can speed up replacement in this example?

I would be grateful for any help.

Answer 1

What you have laid out is inefficient in a number of ways. But, I will tackle your question, as is, first, then I will lay out how I would approach it.

Since your bottleneck is

Table[{j dx, k dy, p[j, k] /. sol}, {j, 1, m}, {k, 1, m}]

you need to use a more efficient layout for your replacement rules: Dispatch. Replace

sol = Solve[eqns, vars][[1]];

with

sol = Dispatch@Solve[eqns, vars][[1]];

and you will notice an enormous speed up.

As to how I would approach this problem, I would recast it as a matrix problem. The simplest method for doing that is to use CoefficientArrays which returns something of the form {b, m}. Then, this problem becomes

{b, m} = CoefficientArrays[eqns, vars];
solution = LinearSolve[m, b];

where solution is your p[j,k] array flattened out.