# An efficient way to translate a triple for-loop from Matlab to Mathematica

This is my first post on StackExchange, so if I violate any etiquettes, I apologize in advance. I am a senior engineering student who has been using Mathematica for a year now. For one of my classes, called reservoir simulation, I am trying to obtain the pressure distribution matrix. My instructor exclusively uses Matlab to model the problems we do in class. I am not very good at Matlab, so I decided to use Mathematica.

The instructor gave us his code in Matlab, and I want to translate it into Mathematica. There are a lot of constants (All of them are a 12 x 9 matrix) in the problem I managed to get all of them into my Mathematica notebook, so they are all defined. What I will not share are the constant but assume they are already in the notebook. Below is the Matlab code my instructor used which I am explicitly trying to translate.

First, he initialized many matrices like this:

P_old = ones(9,12) * 7750;

Pwf = 5600 ; % psia

Pe = zeros(9,12);

E = zeros(9,12);

Ge = zeros(9,12);

SGe = zeros(9,12);

Pw = zeros(9,12);

W = zeros(9,12);

Gw = zeros(9,12);

SGw = zeros(9,12);

Pn = zeros(9,12);

N = zeros(9,12);

Gn = zeros(9,12);

SGn = zeros(9,12);

Ps = zeros(9,12);

S = zeros(9,12);

Gs = zeros(9,12);

SGs = zeros(9,12);

omega = zeros(9,12); % productivity index

omega_history = zeros(9,12,61);

Store = zeros(9,12,61);

q = zeros(9,12);

q(3,9) = -650;

dt = 1;

A = zeros(108,108);

B = zeros(108,1);

gamma = zeros(9,12);

nx = 12;

ny = 9;


Then he made a very long triple for-loop. The indices for this for-loop are problem specific which is what I am trying to replicate.

for t = 1:61

Store(:,:,t) = P_old;

for j = 1:9

for i = 1:11

Pe(j,i) = (Vb(j,i+1) * P_old(j,i+1) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j,i+1));

E(j,i) = (2 * kx1(j,i+1) * Ax(j,i+1) * kx1(j,i) * Ax(j,i) ) / ((kx1(j,i+1) * Ax(j,i+1) * x(j,i)) + (kx1(j,i) * Ax(j,i) * x(j,i+1))) * (1 + 9.0 * 10^(-6) * (Pe(j,i) - 14.7)) * 1/(  a * (Pe(j,i))^(3) + b * (Pe(j,i))^(2) + c * Pe(j,i) + d  );

Ge(j,i) = G(j,i+1) - G(j,i);

SGe(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Pe(j,i) - 14.7))) );

end
end

Pe(isnan(Pe)) = 0;

E(isnan(E)) = 0;

SGe(isnan(SGe)) = 0;
for j = 1:9

for i = 2:12

Pw(j,i) = (Vb(j,i-1) * P_old(j,i-1) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j,i-1));

W(j,i) = (2 * kx1(j,i-1) * Ax(j,i-1) * kx1(j,i) * Ax(j,i) ) / ((kx1(j,i-1) * Ax(j,i-1) * x(j,i)) + (kx1(j,i) * Ax(j,i) * x(j,i-1))) *  (1 + 9.0 * 10^(-6) * (Pw(j,i) - 14.7)) * 1/(  a * (Pw(j,i))^(3) + b * (Pw(j,i))^(2) + c * Pw(j,i) + d  );

Gw(j,i) = G(j,i-1) - G(j,i);

SGw(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Pw(j,i) - 14.7))) );
end
end
Pw(isnan(Pw)) = 0;

W(isnan(W)) = 0;

SGw(isnan(SGe)) = 0;
for j = 2:9

for i = 1:12

Pn(j,i) = (Vb(j-1,i) * P_old(j-1,i) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j-1,i));

N(j,i) = (2 * ky1(j-1,i) * Ay(j-1,i) * ky1(j,i) * Ay(j,i) ) / ((ky1(j-1,i) * Ay(j-1,i) * y(j,i)) + (ky1(j,i) * Ay(j,i) * y(j-1,i))) * (1 + 9.0 * 10^(-6) * (Pn(j,i) - 14.7)) * 1/(  a * (Pn(j,i))^(3) + b * (Pn(j,i))^(2) + c * Pn(j,i) + d  );

Gn(j,i) = G(j-1,i) - G(j,i);

SGn(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Pn(j,i) - 14.7))) );

end
end

Pn(isnan(Pn)) = 0;
N(isnan(N)) = 0;
SGn(isnan(SGn)) = 0;

for j = 1:8

for i = 1:12

Ps(j,i) = (Vb(j+1,i) * P_old(j+1,i) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j+1,i));

S(j,i) = (2 * ky1(j+1,i) * Ay(j+1,i) * ky1(j,i) * Ay(j,i) ) / ((ky1(j+1,i) * Ay(j+1,i) * y(j,i)) + (ky1(j,i) * Ay(j,i) * y(j+1,i))) * (1 + 9.0 * 10^(-6) * (Ps(j,i) - 14.7)) * 1/(  a * (Ps(j,i))^(3) + b * (Ps(j,i))^(2) + c * Ps(j,i) + d  );

Gs(j,i) = G(j+1,i) - G(j,i);

SGs(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Ps(j,i) - 14.7))) );

end
end

Ps(isnan(Ps)) = 0;
S(isnan(S)) = 0;
SGs(isnan(SGs)) = 0;

if P_old(4,4) <= Pwf

omega(4,4) = 0;
else

omega(4,4) = ( (2 * pi .* sqrt(kx1(4,4) * ky1(4,4)) .* thickness(4,4) * (1 + 9.0 * 10^(-6) * (P_old(4,4) - 14.7))) / ( (  a * (P_old(4,4))^(3) + b * (P_old(4,4))^(2) + c * P_old(4,4) + d  ) * ( log(re/0.25) ) ) ) ;

end

omega_history(:, :, t) = omega;

for j = 1:9

for i = 1:12

C(j,i) = - (E(j,i) + W(j,i) + N(j,i) + S(j,i) + gamma(j,i) + omega(j,i));

C(C==0) = 1;

end
end

for j = 1:9

for i = 1:12

Q(j,i) = (-1 .* omega(j,i) .* Pwf) - (gamma(j,i) .* P_old(j,i)) + ( E(j,i) .* SGe(j,i) .* Ge(j,i) )  + ( W(j,i) .* SGw(j,i) .* Gw(j,i) )  + ( N(j,i) .* SGn(j,i) .* Gn(j,i) )  + ( S(j,i) .* SGs(j,i) .* Gs(j,i) )  - q(j,i);

end
end

Q(isnan(Q)) = 0;

for j = 1:9

for i = 1:12

r = i + nx * (j-1);

A(r,r) = C(j,i);

if r >1

A(r,r-1) = W(j,i);

end

if r < nx * ny

A(r,r+1) = E(j,i);
end

if r + nx <= nx * ny

A(r,r+nx) = S(j,i);

end

if r > nx

A(r,r-nx) = N(j,i);

end

B(r) = Q(j,i);

end
end

P_new = (A\B);

P_new = reshape(P_new,[12,9]);
P_new = P_new';

P_old = P_new;

end


As you can see above, it is a very involved messy code. I attempted to do this in Mathematica using the "Module" function and the "Do" function, but I cannot replicate the results. My futile attempt is the following:

Module[{i, j}, Do[pe = (volume[[j, i + 1]] pi[[j, i + 1]] +
volume[[j, i]]  pi[[j, i]])/(volume[[j, i]] +
volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet, {i, 1, 8, 1}, {j, 1, 10, 1}]]


This is just one small for-loop of the bigger for-loop, and I do not get what I want.

The end (P_new) result should be something like this (The "*" represents an indeterminate): • Can you provide a smaller yet essence-preserving example? – Αλέξανδρος Ζεγγ Nov 11 '18 at 6:56
• hey @H.Alanzi welcome! Are you sure "(All of them are a 12 x 9 matrix)" because most of zeros and ones are called with arguments (9, 12) which-if I'm not mistaken-means that they are the transpose of what the quoted claim above states (please check the documentation)... – user42582 Nov 11 '18 at 8:22
• If you just want to make a naive translation, I don't think there'll be any barrier because the MATLAB code only involves basic algebraic calculation, you must have made simple mistake if you don't get what you want. If you want to implement the algorithm in a better way (make the code more elegant, efficient, etc. ) then you should show us the algorithm in traditional math notation, because it's hard to figure it out directly from the code. – xzczd Nov 11 '18 at 8:43

Okay, let's focus on this snippet:

Module[{i, j},
Do[
pe = (volume[[j, i + 1]] pi[[j, i + 1]] +
volume[[j, i]]  pi[[j, i]])/(volume[[j, i]] +
volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet,
{i, 1, 8, 1}, {j, 1, 10, 1}]]


Creating some dummy data:

volume = RandomReal[{-1, 1}, {10, 9}];
pi = RandomReal[{-1, 1}, {10, 9}];


You use Do here and not For. That's indeed great for many reasons and shows that you already know about the problems with For. Moreover, Do scopes its iterators, hence we don't need Module.

So let's strip it off. Moreover, by default, Do increases its iterators by 1

Do[
pe = (volume[[j, i + 1]] pi[[j, i + 1]] +
volume[[j, i]] pi[[j, i]])/(volume[[j, i]] +
volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet,
{i, 1, 8}, {j, 1, 10}]


This reassigns pe over and over again. But I got the impression that you mean pe to be a matrix. The "Matlab-way" of constructing it would be to allocate a 0-matrix and to fill it afterwards (actually, a good Matlab programmer would do it differently):

pe = ConstantArray[0, {10, 8}];
Do[pe[[j, i]] = (volume[[j, i + 1]] pi[[j, i + 1]] +
volume[[j, i]] pi[[j, i]])/(volume[[j, i]] +
volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet,
{i, 1, 8, 1}, {j, 1, 10, 1}]


In Mathematica, we can use Table instead.

pe = Table[(volume[[j, i + 1]] pi[[j, i + 1]] +
volume[[j, i]] pi[[j, i]])/(volume[[j, i]] +
volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet,
{j, 1, 10, 1}, {i, 1, 8, 1}]


Notice that I had to reorder the iterators in order to produce the same matrix as before and not its transpose.

Both Matlab and Mathematica are interpreted languages (in contrast to compiled languages like C or fortran). And in interpreted languages, it is usually way more performant to avoid loop constructions if possible and to replace them by vectorized code. Your instructor must have aimed at showing you how not to program in Matlab. Indeed, the code is overly convoluted and almost all interior loops could be replaced by vectorized operations (at the first glance).

So, in this example, the following produces the same result, and is one to two orders of magnitude faster:

pe = With[{a = volume pi},
(a[[;; , ;; -2]] + a[[;; , 2 ;;]])/(volume[[;; , ;; -2]] + volume[[;; , 2 ;;]])
]/. Indeterminate -> "*" // Quiet


This should also lead to much more concise and more legible code.

Another example (it is so striking that I believe that the instructor did that intentionally):

for j = 1:9

for i = 1:12

C(j,i) = - (E(j,i) + W(j,i) + N(j,i) + S(j,i) + gamma(j,i) + omega(j,i));

C(C==0) = 1;

end
end


would be simply (in Matlab syntax):

C = - (E + W + N + S + gamma + omega);
C(C==0) = 1;


Moreover, putting C(C==0) = 1; into the body of the loop is just plain crazy; it has to be performed just once.

# Remarks

For the understanding of ;; see Span; it is essentially the analogue to Matlab's :. (Notice the three-argument versions of these functions differ in the ordering of the arguments.) See also this post for more analogies between Matlab and Mathematica.

• “instructor must have aimed at showing you how not to program in Matlab” – I like your optimism, but going by what atrocities I've seen in quite a few other Matlab scripts, I wouldn't be so sure... – leftaroundabout Nov 11 '18 at 23:59
• Is there some packages that we can convert matlab script to mma script, it would be great helpful for many MMA users. – ABCDEMMM Nov 17 '18 at 23:01
• @ABCDEMMM Unfortunately, I don't think that such a package exists. It would be really hard to develop. But it is possible to run Matlab code from Mathematica with MatLink. – Henrik Schumacher Nov 17 '18 at 23:10
• @HenrikSchumacher, thanks a lot, I have some code about topo, still I want to use it in MMA. okay, thanks a lot again! – ABCDEMMM Nov 18 '18 at 0:09
• @ABCDEMMM You're welcome. Glad to be of help. – Henrik Schumacher Nov 18 '18 at 8:15