Problem with for loop within a loop of a loop

The following code calculates the excess pore water pressure u after consolidation, and generates a m x n matrix, where the thickness of the soil layer i divided into m equal parts and u is calculated for n time steps.

In this simplified exampel, beta is a constant, but I want to impliment specifik values of beta for different equal sized depths indicated by m.

So say I have a beta-vector containing 5 values and a soil layer divided into m=50 equal parts. beta(1) will run from m(1):m(10), beta(2) will run from m(11):m(20) and so forth.

Does anyone have a suggestion on how to impliment the beta-vector?

beta=0.2; n= 10; m= 5
u_1=[60 54 41 29 19 15]';
u(:,1)=u_1;
for j=(1:n)
  for i=(2:m) 
    u(i,j+1)=u(i,j)+beta*(u(i-1,j)+u(i+1,j)-2*u(i,j));
    for i=(m+1) % As we look at a half-closed layer
      u(i,j+1)=u(i,j)+beta(end)*(2*u(i-1,j)-2*u(i,j))
    end
  end
end

(Ultimatly I will end up with lenght(beta)=30, m= 300 and n=>130 and even with beta as a constant it takes some time to run.)

Much appriciated!

[Walter Roberson]

for i=(2:m) 
    u(i,j+1)=u(i,j)+beta*(u(i-1,j)+u(i+1,j)-2*u(i,j));
for i=(m+1) % As we look at a half-closed layer
    u(i,j+1)=u(i,j)+beta(end)*(2*u(i-1,j)-2*u(i,j))   %HERE
end
end

On the line I marked with HERE, does i refer to for i=(m+1) or does it refer to the enclosing for i=(2:m) ?

[Mette Marie Bondegaard Petersen]

Ohh yes - actually the matrix ends up with the dimensions (m+1) x (n+1) for the lower boundary layer and initial conditions at t=0 respectively.

So it is within the loop of i=(2:m) that I need to impliment the beta-vector

Does that make sense?

[Walter Roberson]

Perhaps

for i=(2:m) 
    u(i,j+1)=u(i,j)+beta*(u(i-1,j)+u(i+1,j)-2*u(i,j));
end
for i=(m+1) % As we look at a half-closed layer
    u(i,j+1)=u(i,j)+beta(end)*(2*u(i-1,j)-2*u(i,j));
end

?

[Mette Marie Bondegaard Petersen]

Well as I see it, this gives me the same result as before.

I still have trouble implimenting beta as a vector, with specified values for different equal sized depths indicated by m.

Am I missing something or is my initial question unclear?

..and thank you for taking the time!

EDIT:

Ok. So now I get it. Your fix made it run much much better, so thank you for that. One down, one to go :)

[Walter Roberson]

Well as I see it, this gives me the same result as before.

Yes, but faster! ;-)

[Mette Marie Bondegaard Petersen]

Very helpful! :)

[Guillaume]

Note that

for i = (m+1)  %why the () ?

is simply

i = m+1;

[Mette Marie Bondegaard Petersen]

Understood :)

This is a result of me trying to figure out how to run the loop for different intervals of m, corresponding to different values of beta, as the actual beta is a vector.

Answer 1

[Guillaume]

Probably, the biggest source of slow down is the lack of preallocation of u. As a result, it grows one column at a time, necessating reallocation and copying on each j step.

u = zeros(m+1, n+1);  %preallocation
u(:, 1) = u1;
for j = 1:n   %why the () ?
    for i = 2:m   %so u(1, j) is 0 for all but j = 1? 
        u(i,j+1)=u(i,j)+beta*(u(i-1,j)+u(i+1,j)-2*u(i,j));
    end
    u(m+1,j+1)=u(m+1,j)+beta*(2*u(m,j)-2*u(m+1,j));  %not sure why there was a beta(end) here.
end

[Mette Marie Bondegaard Petersen]

I have been playing around with the code in order to implement beta as a vector.

The preallocation makes sense – thank you.

Yes, u(1,j) is 0 for all but j=1, with the assumption that the initial pressure instantaneously becomes zero on the upper permeable boundary.

As for the last loop i=m+1, this calculates the pressure, u on the lower closed boundary. As I am trying to implement beta as a vector, the lower boundary is calculated from the last value of beta, i.e. beta(end).

[Guillaume]

I completely missed your question about beta.

The simplest thing is to repelem your different beta across the rows. e.g.

beta = [0.2, 0.3, 0.5, 0.7, 0.9];
ubeta = repelem(beta(:), m/numel(beta));

and then you use ubeta(i) in your equation.

Also, I've just noticed that you don't need the i loop:

for j = 1:n
    u(2:end-1, j+1) = u(2:end, j) + ubeta(2:end) .* (u(1:end-2, j) + u(3:end, j) - 2*u(2:end-1, j));
    u(end, j+1) = u(end, j) + 2*ubeta(end) *(u(end-1, j) - u(end, j));
end

[Mette Marie Bondegaard Petersen]

I see what you did there - nice, that will work.

Thank you so much!