Problem with for loop within a loop of a loop

1 view (last 30 days)
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!
  8 Comments
Guillaume
Guillaume on 7 Sep 2019
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.

Sign in to comment.

Accepted Answer

Guillaume
Guillaume on 7 Sep 2019
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
  3 Comments
Guillaume
Guillaume on 8 Sep 2019
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!

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Tags

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!