Code is running continuously, but never ends

h_s=[1 2 5 6];
h_d = [3 4 8 6];
MA = [3.5 1.5 2.5];
D_d = 0.0254*[1.25 1.75 2.25];
D_s = [0.03175 0.04445 0.05715];
D_p = 0.0254*[2.5 3.5 4.5];
L = [0.10 0.18 0.16];







for k1 = 1:numel(h_s)  % suction head m
    for k2 = 1:numel(h_d) % delivery head m
        for k3 = 1:numel(MA) % mechanical advantage
            for k4 = 1:numel(D_d) %delivery pipe diameter
                for k5 = 1:numel(D_s) %suction pipe diamter...0.0254*[1.25 1.75 2.25]
                    for k6 = 1:numel(D_p) %piston diameter m
                        for k7 = 1:numel(L)% height of cylinder stoke
                            
                            
                            

y=[0.15 0.094]

%yv = linspace(min(y), max(y), 10);
v0 = 0.01;                      
m =64.5; 
x = 0.7; 
b = 0.625;
a = 0.250; 
g = 9.8;
[A_d(k4,:)] = (D_d(k4)^2)*3.14*0.25; 
                            
                            for a1= 1:numel(A_d)
                                
                                [A_s(k5,:)] = (D_s(k5)^2)*3.14*.25; 
                                for a2= 1:numel(A_s)
                                    

l_d = 4.572; 
l_s = 6.096;
[A_p(k6,:)] = ((D_p(k6))^2)*3.14*0.25; 
                                    
                                    for a3= 1:numel(A_p)
                                        
m_p = 0.2; 
rho = 1000; 
uk = 0.06; 
uf = 0.4; 
up = 0.00089; 
h = 0.001;
t_w=  0.05;  
[A_w(k6,:)] = t_w*D_p(k6)*3.14; 
                                        for a4= 1:numel(A_w)
                                            
[z(k3,:)]=(0.2946*MA(k3)-0.1761);
                                            for a5= 1:numel(z)
                                                
M_tr= 2;
L_tr=0.5;
P=173.7;
R_pin=0.05;]
[M_T(k7,k3,a1,a3,a5,:)]=2*m_p*(A_d(a1)/A_p(a3))^2+(L(k7)+l_d+l_s)*rho*((A_d(a1))^2/A_p(a3))+(2*M_tr*L_tr^2/a^2)+(z(a5))^2*MA(k3)^2*m*(A_d(a1)/A_p(a3))^2;
                                                
                                                for a6= 1:numel(M_T)
                                                    
                                                    [c0(a3,a1,a6,:)]=((1-2*x)*A_p(a3)*P/(A_d(a1)*M_T(a6)));
                            
                            disp(c0);
                            
                                                    for  t0 = 1:numel(c0)
                                                        
                                

                                    [c1(k1,k2,k7,k3,a2,a1,a3,a6,:)]=(rho*g*((h_s(k1))*A_s(a2)+(h_d(k2))*A_d(a1)+(L(k7))*A_p(a3))+z(+1)*m*g*MA(k3))/M_T(a6);
                                                        
                                                        for  t1 = 1:numel(c1)
                                                            
                                    [c2(a3,a6,:)]=(2*rho*g*A_p(a3))/M_T(a6);
                                                            for  t2 = 1:numel(c2)
                                                                
                                                                
                                       [c3(a1,a6,a3,k7,:)]=(2*uf*b+25.12*up*L(k7)*h)*A_d(a1)/(h*M_T(a6)*A_p(a3));
                                        
                                                                for  t3 = 1:numel(c3)
                                                                    
                                           [c4(k4,k5,a1,a2,a6,:)]=(0.124*rho^0.75*up^0.25*[l_d*D_d(k4)^0.25+l_s*D_s(k5)^0.25]*(A_d(a1)/A_s(a2))^1.75)/M_T(a6);
                                            
                                                                    for  t4 = 1:numel(c4)
                                                                        
                                                                        
                                               [c5(a1,a2,a3,a6,:)]=rho*[((A_d(a1))^2/A_s(a2))*[0.5+0.5+0.2+0.4+0.5]+0.5*((A_d(a1))^2/A_p(a3))+A_d(a1)*(0.5+0.4+0.2+1)+0.5*((A_d(a1))^2/A_p(a3))]/(2*M_T(a6));%please check the constants.
                                                
                                                                        for  t5 = 1:numel(c5)
                                                                        
                                                    [y,v(t0,t1,t2,t3,t4,t5,:)] = ode45((@(y,v)((c0(t0)/v^2)-(c1(t1)/v)-(c2(t2)*y/v)-c3(t3)-c4(t4)*v^0.75-c5(t5)*v)) ,y ,v0);
                                                    %disp(v,y) % displays y and x_velocity
                                                
                                                                        end
                                                                    end
                                                                end
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end



%% I am not able to get output and the code is not showing any error.
%% It is going to run for long time.
%% Please, upload this into your matlab tool, then you  can have idea of problem. Please, tell me how to proceed.

[Walter Roberson]

The line

R_pin=0.05;]

is invalid.

[Adam]

How long is 'continuously'? You have 19 nested for loops, unless I miscounted! That's not going to be especially performant!

You can stick a breakpoint in while code is running to see where it currently is. If you do this initially in the deepest for loop, after some amount of time you will get an idea how long it is taking. Or you could use

doc waitbar

around the outer loop to also get an idea of how long it is taking.

[SUSHANT CHAVAN]

I am not getting, where do I add following command? @Adam

Will it before last 'end' line?

doc waitbar

[KALYAN ACHARJYA]

Adam suggested to create a waitbar

https://in.mathworks.com/help/matlab/ref/waitbar.html

[Walter Roberson]

Go to the MATLAB command line and give the command

doc waitbar

and read what shows up.

[SUSHANT CHAVAN]

@Adam

I have studied it step by step...It is not coming out of last loop..How to tackle?

Answer 1

[Guillaume]

As others have said, you need to debug the code you've written, using the various tools that matlab offer. waitbar is one, you could also use the debugger to step through your problem line by line.

I would suspect that the code is not even doing what you meant it to do. As Adam said you have an extremely large number of nested loops and some of them don't make much sense. For example, we have

%inside the k1, k2, k3, ..., k6, k7 loops 
[A_d(k4,:)] = (D_d(k4)^2)*3.14*0.25; 
for a1= 1:numel(A_d)
    %...

Note that this code is inside the k5, k6 and k7 loops, yet doesn't depend at all on these indices. So for a start the whole calculation will be repeated max(k5)*max(k6)*max(k7) times each time yielding the exact same result. Certainly a waste of time, but not the biggest problem.

The right hand side of the equation is always going to be scalar, so the : on the left hand side is misleading. It's equivalent to 1. We also have some unnecessary [] brackets. None of this is a problem, but it makes me wonder if something was intended all along.

The biggest problem is that A_d grows at each step of the k4 loop. When k4 is 1, A_d is only one element, then when k4 is 2, you add a new element as A_d(2), etc. So, the number of steps that the a1 loop is going to take depends on the current k4. I'm fairly certain that was not the intent.

[SUSHANT CHAVAN]

Actually, what I tried to do in this code is, I am storing the values of constants A_d, A_s.......etc. in the matrix form and calling a particular value in the equation which includes that particular constant.

I wanted to call each value from the created vectors by using 'for ' loops one by one and find the value of 'v', which is the solution of differential equation.

I have created such for loops to store the values of each individual constants in matrix form, then I could call it in differential equation.

As 'Guillaume' suggested, A_d grows at each step, but I am taking the value of D_d from the matrix given at the start one by one and storing it in A_d matrix.

Then whenever A_d comes in expressions used later, A_d would be an single scaler value. I want to call a individual A_d(Area) for corressponding D_d(diameter). That was the approach.

Please, tell your thought about this?

[Guillaume]

It's not clear to me how many independent variables you have. If the only independent variables are just: h_s, h_d, MA, D_d, D_s, D_p and L (so 7 variables, then you should have just 7 loops. (or just one loop with a different design). The A_d, etc., if dependent should not trigger another loop, you use the same index as used for the variable that depends to.

In addition, there's no reason to calculate A_d inside the loop, you can calculate all values at once before any loop, it's simply:

A_d = D_d .^ 2 * pi / 4;  %use pi instead of 3.14, it's more accurate. Note the use of .^ since D_d is a vector
%same for A_S and probably others I haven't looked at.

If I understood, correctly what you're trying to do, the code should be:

%independent variables
h_s = [1, 2, 5, 6];
h_d = [3, 4, 8, 6];
MA = [3.5, 1.5, 2.5];
%etc. all the variables to loop over

%constants, they shouldn't be in the loops since they're constant.
v0 = 0.01;                      
m = 64.5; 
x = 0.7; 
%etc. they're peppered all over your code, it's a mess.

%dependent variables
A_d = D_d .^ 2 * pi / 4;  %will be indexed by k4
A_s = D_s .^ 2 * pi / 4;  %will be indexed by k5
%etc.

%preallocation of results, 7-D arrays
M_T = zeros(numel(h_s), numel(h_d), numel(MA), ...etc for all 7 variables);
c0 = M_T;  %or repeat the zeros lines above
%etc. all the output variables

%the loops
for k1 = 1:numel(h_s)  % suction head m
    for k2 = 1:numel(h_d) % delivery head m
        for k3 = 1:numel(MA) % mechanical advantage
            for k4 = 1:numel(D_d) %delivery pipe diameter
                for k5 = 1:numel(D_s) %suction pipe diamter...0.0254*[1.25 1.75 2.25]
                    for k6 = 1:numel(D_p) %piston diameter m
                        for k7 = 1:numel(L)% height of cylinder stoke
                            M_T(k1, k2, k3, k4, k5, k6, k7) = 2*m_p * (A_d(k4)/A_p(k6))^2 + (L(k7)+l_d+l_s)*rho*((A_d(k4))^2/A_p(k6))+(2*M_tr*L_tr^2/a^2)+(z(k3))^2*MA(k3)^2*m*(A_d(k4)/A_p(k6))^2;
                            %etc. everything uses k1 to k7
                        end
                    end
                end
            end
        end
    end
end

Another option is to ndgrid all your inputs and just use a linear index:

%create 7-d arrays:
[h_s, h_d, MA, D_d, D_s, D_p, L] = ndgrid([1, 2, 4, 5],    ...suction head m
                                          [3, 4, 8, 6],    ...delivery head m
                                          [3.5, 1.5, 2.5], ...mechanical advantage
                                          0.0254*[1.25, 1.75, 2.25], ...delivery pipe diameter
                                          [0.03175, 0.04445, 0.05715], ...suction pipe diameter
                                          0.0254*[2.5, 3.5, 4.5], ...piston diameter
                                          [0.10, 0.18, 0.16]); %stroke

%constants as before
v0 = 0.01;                      
m =64.5; 
x = 0.7; 
%etc.

%dependent variables as before:
A_d = D_d.^2 * pi/4;
%etc.

%no loop needed for M_T and other variables. Calculate everything at once as 7_d array
%need to use dotted operators when dealing with matrices
M_T =  2*m_p*(A_d./A_p).^2 + (L+l_d+l_s)*rho.*(A_d.^2./A_p) + 2*M_tr*L_tr^2/a^2 +z.^2.*MA.^2*m.*(A_d./A_p).^2;
%etc. for the c variables

%preallocate v
v = zeros(size(h_s)); %as a 7_d variable
%loop over each element using linear indexing
for k = 1:numel(v)
    [~, v(k)] = ode45(@(y,v) c0(k)/v^2 - c1(k)/v - c2(k)*y/v -c3(k) - c4(k)*v^0.75 - c5(k)*v, y , v0);
end

[SUSHANT CHAVAN]

Thank you so much for giving your time for this. I will try by above method and get back to you. Thanks.

[SUSHANT CHAVAN]

h_s=[1 2 5 6];
h_d = [3 4 8 6];
MA = [3.5 1.5 2.5];
D_d = 0.0254*[1.25 1.75 2.25];
D_s = [0.03175 0.04445 0.05715];
D_p = 0.0254*[2.5 3.5 4.5];
L = [0.10 0.18 0.16];
%%
y=[0 L]
%I want to select different interval for each loop. It means, for first run
%y shoild be [0 0.10], for second run y should be [0 0.18] and so on. How
%do I do it?
yv = linspace(min(y), max(y), 20);
v0 = 0.01;                         
m =64.5; 
x = 0.7; 
b = 0.625; 
a = 0.250; 
g = 9.8; 
l_d = 4.572; 
l_s = 6.096; 
m_p = 0.2; 
rho = 1000; 
uk = 0.06; 
uf = 0.4; 
up = 0.00089; 
h = 0.001; 
t_w=  0.05;   
M_tr= 2;
L_tr=0.5;
P=173.7;
R_pin=0.05;

%%
A_d = D_d.^2* pi / 4 ;
A_s = D_s.^2 * pi / 4 ;
A_p = D_p.^2 * pi / 4;
A_w = t_w * D_p.* pi ;
                                        
%%
M_T = zeros(numel(MA), numel(A_d), numel(A_p), numel(L) );
c0 = zeros(numel(MA), numel(A_d), numel(A_p), numel(L));
c1 = zeros(numel(MA), numel(A_d), numel(A_p), numel(L), numel(A_s), numel(h_s), numel(h_d));
c2 = zeros(numel(MA), numel(A_d), numel(A_p), numel(L));
c3 = zeros(numel(MA), numel(A_d), numel(A_p), numel(L));
c4 = zeros(numel(MA), numel(A_d), numel(A_p), numel(L), numel(A_s));
c5 = zeros(numel(MA), numel(A_d), numel(A_p), numel(L), numel(A_s));
%%

for k1 = 1:numel(h_s)  
    for k2 = 1:numel(h_d) 
        for k3 = 1:numel(MA) 
            for k4 = 1:numel(D_d) 
                for k5 = 1:numel(D_s) 
                    for k6 = 1:numel(D_p) 
                        for k7 = 1:numel(L)
                            z(k3) = 0.2946 * MA(k3) - 0.1761 ;                        
                            M_T(k4,k6,k7,k3) = 2*m_p * (A_d(k4)/A_p(k6))^2 + (L(k7)+l_d+l_s)*rho*((A_d(k4))^2/A_p(k6))+(2*M_tr*L_tr^2/a^2)+(z(k3))^2*MA(k3)^2*m*(A_d(k4)/A_p(k6))^2;                                               
                            c0(k6,k4,k7,k3)=((1-2*x)*A_p(k6)*P)/(A_d(k4)*M_T(k4,k6,k7,k3));
                            if y < L(k7)/2
                                c1(k1,k5,k2,k4,k7,k6,k3)=(rho*g*((h_s(k1))*A_s(k5)+(h_d(k2))*A_d(k4)+(L(k7))*A_p(k6))- z(1)*m*g*MA(k3))/M_T(k4,k6,k7,k3) ;
                            end
                            if y > L(k7)/2
                                c1(k1,k5,k2,k4,k7,k6,k3)=(rho*g*((h_s(k1))*A_s(k5)+(h_d(k2))*A_d(k4)+(L(k7))*A_p(k6))+ z(1)*m*g*MA(k3))/M_T(k4,k6,k7,k3) ;
                            end
                            %I want to take value of c1 according to
                            %value of l, it means if y is less than L/2
                            %I will use first value given above
                            %otherwise second conditions apply nd these
                            %particular value of c1 will be used in
                            %diff equation
                            %c1(k1,k5,k2,k4,k7,k6,k3)=(rho*g*((h_s(k1))*A_s(k5)+(h_d(k2))*A_d(k4)+(L(k7))*A_p(k6))+ z(1)*m*g*MA(k3))/M_T(k4,k6,k7,k3) ;
                            c2(k6,k4,k7,k3)=(2*rho*g*A_p(k6))/M_T(k4,k6,k7,k3);
                            c3(k7,k4,k6,k3)=(2*uf*b+25.12*up*L(k7)*h)*A_d(k4)/(h*M_T(k4,k6,k7,k3)*A_p(k6));
                            c4(k4,k5,k6,k7,k3)=(0.124*rho^0.75*up^0.25*[l_d*D_d(k4)^0.25+l_s*D_s(k5)^0.25]*(A_d(k4)/A_s(k5))^1.75)/M_T(k4,k6,k7,k3);
                            c5(k4,k5,k6,k7,k3)=rho*[((A_d(k4))^2/A_s(k5))*[0.5+0.5+0.2+0.4+0.5]+0.5*((A_d(k4))^2/A_p(k6))+A_d(k4)*(0.5+0.4+0.2+1)+0.5*((A_d(k4))^2/A_p(k6))]/(2*M_T(k4,k6,k7,k3)) ;%please check the constants.
                            [y ,v] = ode45((@(y,v)((c0(k6,k4,k7,k3)/v^2)-(c1(k1,k5,k2,k4,k7,k6,k3)/v)-(c2(k6,k4,k7,k3)*y/v)- c3(k7,k4,k6,k3)- c4(k4,k5,k6,k7,k3)*v^0.75-c5(k4,k5,k6,k7,k3)*v)) ,yv ,v0);
                            %Whatever value of 'v' I get from the solution of
                            %this diff equation is just 20*1 matrix, but I want
                            %it to 20*no. of combinations according to values
                            %of independent variables like h_s, h_d, D_d...etc.
                            %In short, I don't just want last value of 'v',
                            %rather I want the value of 'v' for each
                            %combination. Please, suggest solution?
                        end
                    end
                end
            end
        end
    end
end

@Guillaume, As per your suggestion, I have applied the first method(not ndgrid, before that)

Please, see the the comments, I incorporated in the code. You will get about my concerns. Thanks.

[Guillaume]

"I want to select different interval for each loop"

For y? If so, simply define y inside the loops as

y = linspace(0, L(k7), 20);

However if y is a vector, then your if y < ... and if y > ... (and what happens when y = ... ?) don't make any sense.

But more importantly, I don't understand what the purpose of y is. Your ode is supposed to solve but the y you've created is not used by the ode (it's a completely different y variable internal to the solver). If the value of c1 is supposed to change within the ode depending on the y internal to the solver then you'll have to completely change the ode design.

[SUSHANT CHAVAN]

Actually, y is the range/ domain over which ode would supposed to run. I divided y into 20 parts and calculate 'v ' for each part. L is the length of the piston. I have to calculate instantaneous velocity of the piston, for this I divided L into 20 parts.

If I remove the 'y' from the code, it will not going to run.

and for "if y < ... and if y > ..." this case, I want to take value of y . Means, I have cut L into 20 parts, each one constitute as y...If a y< L/2...at a particular stage I will have some particular value of c1..and if y>L/2 I will have some other value of c1. As our y is proceeding up into 20 parts so its relation with L/2 also going to change, somewhere it will be greater than L/2 and vice versa.