How can I reduce the computation time for this MATLAB code? It always shows busy without giving any output

%Input variables 
clear all
close all
l=2;%metre
E=210*10^9;%Pa
do=0.127;%5in=0.127m
di=0.1086;%4.276in=0.1086m
I=(pi/64)*(do^4-di^4);
A=(pi/4)*(do^2-di^2);
wnom=29.0192;%19.5 lb/ft = 29.012 kg/m Nominal_Weight
Mass =(wnom)*l;
Volume=A*l;
Area=do*l;
rho=Mass/Volume;
P=0;
cd=0;
%ksp=0;
ksp=2.0201*10^9;
mu=8.9*10^(-4);
rhowater=1000;%kg/m3
weightofair=0.3587;%N
sol=[5.2036 4.5721;1.8105 4.6666;6.5143 4.3060;6.4432 8.0821;3.3297 5.0217;3.0705 5.0641];
response=0;
response1r=0;
response1i=0;
j=1;
%i=1;
t=1000;
h=0.1;
for i=1
%for x=0:l/4:l
for h=0.1:0.2:0.5
%     response=0;
%     response1r=0;
%     response1i=0;
%     figure
%for i=1:4
%for P=0:-5000:-10000   
    temp=sol(i,1);
    s(j,1)=temp;
    temp1=sol(i,2);
    s(j,2)=temp1;
    alpha=(temp)/l;
    s(j,3)=alpha;
    beta=(temp1)/l;
    s(j,4)=beta;
%     if h == 0
%         w2=(weightofair)/l;  
%     else
%         Masswater=rhowater*Area*h*0.0254;
%         w2=(Masswater*9.81)/l;   
%     end
%    delstat=(1/(E*I))*(-((wnom*(l-x)^4)/(24))-((wnom*l^3*x)/(6))+((wnom*l^4)/24));%meter
    w1=h;
    s(j,5)=h;
    if h==0
        cd=0;
    else
        cd=-(mu*Area)/(h*0.0254);
    end
    s(j,6)=cd;
    s(j,7)=ksp;
    s(j,8)=P;
    R1=(rho*A);
    R2=(1i*cd);
    R3=(-E*I*alpha^4)+(P*alpha^2)+ksp;
    finalomega=[R1;R2;R3];
    finalomega1=roots(finalomega);
    s(j,9)=finalomega1(1,1);
    z1=finalomega1(1,1);
    o1=isreal(finalomega1(1,1));
    if o1==0
        z1real=real(finalomega1(1,1));
        z1imag=imag(finalomega1(1,1));
    end    
    s(j,10)=finalomega1(2,1);
    z2=finalomega1(2,1);
    o2=isreal(finalomega1(2,1));
    if o2==0
        z2real=real(finalomega1(2,1));
        z2imag=imag(finalomega1(2,1));
    end 
    if ksp==2.0201*10^9
        omega=z2imag;
    else
        omega=z2;
    end
    k1=((cosh(temp)+cos(temp))/(sinh(temp)+sin(temp)));
    s(j,11)=k1;
    int1=1+(((sin(2*beta*l)/4)+(sinh(2*beta*l)/4)-(cos(beta*l)*sinh(beta*l))-(cosh(beta*l)*sin(beta*l)))/beta);
    int2=((-(k1^2)/(2*beta))*((cos(beta*l)*sin(beta*l))-(2*cos(beta*l)*sinh(beta*l))+(2*cosh(beta*l)*sin(beta*l))-(cosh(beta*l)*sinh(beta*l))));
    int3=((k1*(sin(beta*l)-sinh(beta*l))^2)/beta);
    int=int1+int2-int3;
    const=1/sqrt(rho*A*int);
    s(j,12)=const;
    x=0:l/4:l;
    w=const*((cosh(alpha*x)-cos(beta*x))-(k1*(sinh(alpha*x)-((alpha/beta)*sin(beta*x)))));    
%     sgm1=cos(alpha*l*1i);
%     sgm2=sin(alpha*l*1i);
%     sgm3=cosh(l)+sinh(l);
%     term1=-(24*alpha^8*wnom*sgm1*1i)-(24*beta^8*w*cos(beta*l)*1i)+(24*alpha^2*beta^6*wnom*1i)+(24*alpha^6*beta^2*wnom*1i)+(6*alpha^6*beta^6*l^4*wnom*1i)+(24*alpha^8*wnom*sgm1*cos(beta*l)*1i)+(24*beta^8*wnom*sgm1*cos(beta*l)*1i);
%     term2=-(24*alpha^2*beta^6*wnom*sgm1*1i)-(24*alpha^6*beta^2*wnom*cos(beta*l)*1i)-(24*alpha*beta^7*wnom*sgm2*sin(beta*l))+(24*alpha^7*beta*wnom*sgm2*sin(beta*l))+(24*alpha^3*beta^6*l*wnom*sgm2)-(24*alpha^7*beta^2*l*wnom*sgm2);
%     term3=(24*alpha^2*beta^7*l*wnom*sin(beta*l)*1i)-(24*alpha^6*beta^3*l*wnom*sin(beta*l)*1i)-(12*alpha^4*beta^6*l^2*wnom*sgm1*1i)+(12*alpha^8*beta^2*l^2*wnom*sgm1*1i)-(12*alpha^2*beta^8*l^2*wnom*cos(beta*l)*1i)+(12*alpha^6*beta^4*l^2*wnom*cos(beta*l)*1i);
%     term4=(3*alpha^4*beta^8*l^4*wnom*sgm1*cos(beta*l)*1i)+(3*alpha^8*beta^4*l^4*wnom*sgm1*cos(beta*l)*1i)-(4*alpha^3*beta^8*l^3*wnom*cos(beta*l)*sgm2)+(4*alpha^4*beta^7*l^3*wnom*sgm1*sin(beta*l)*1i)+(4*alpha^7*beta^4*l^3*wnom*cos(beta*l)*sgm2);
%     term5=-(4*alpha^8*beta^3*l^3*wnom*sgm1*sin(beta*l)*1i)-(3*alpha^5*beta^7*l^4*wnom*sgm2*sin(beta*l))+(3*alpha^7*beta^5*l^4*wnom*sgm2*sin(beta*l));
%     term6=-(24*alpha^6*beta^7*sin(beta*l)*sgm3)+(24*alpha^7*beta^6*sgm2*sgm3*1i);
%     etai0=(rho*A*const)*((term1+term2+term3+term4+term5)/(term6));
     etai0=(rho*A)*(w1*const)*(((sinh(alpha*l)-(k1*cosh(alpha*l))+k1)/alpha)-(((beta*sin(beta*l))+(k1*alpha*cos(beta*l))-(alpha*k1))/beta^2));%
    %    w=const*((cosh(beta*x)-cos(beta*x))-(k1*(sinh(beta*x)-sin(beta*x))));
%    etai0=(rho*A)*(w1*const)*(((sinh(beta*l)-sin(beta*l))/beta)-(k1*((cosh(beta*l)+cos(beta*l)-2)/beta)));%Solved Manually
    s(j,13)=etai0;    
if cd==0
        etait=etai0*cos(omega*t);
        response1=(w*etait);
     %   response=response+response1;
        force=ksp*response1;
        for u=1:5
            while response1(1,u) > (h*0.0254) || response1(1,u) < -(h*0.0254)
                ksp=ksp+(10^9);
                response1(1,u)=force(1,u)/ksp;
                as(u,1)=response1(1,u);
                as(u,2)=ksp;
            end
        end
else
    etaitreal=etai0*cos(omega*t);
    etaitimag=etai0.*exp((cd*t)/(rho*A));
    response1real=real(w*etaitreal);
%     response1r=response1r+response1real;
    forcereal=ksp*response1real;  
    for u2=1:5
        while response1real(1,u2) > (h*0.0254)|| response1real(1,u2) < -(h*0.0254)
            ksp=ksp+(10^9);
            response1real(1,u2)=forcereal(1,u2)/ksp;
        end
            asr(u2,1)=response1real(1,u2);
            asr(u2,2)=ksp;
    end     
    response1imag=(w*etaitimag);
    forceimag=ksp*response1imag;
    for u3=1:5
        while response1imag(1,u3) > (h*0.0254)|| response1imag(1,u3) < -(h*0.0254)
            ksp=ksp+(10^9);
            response1imag(1,u3)=forceimag(1,u3)/ksp;
        end
            asi(u3,1)=response1imag(1,u3);
            asi(u3,2)=ksp;    
    end       
end
    if cd==0
        plot(x,response1)
    else
        plot(x,response1imag)
    end
%    title('Dynamic Response of Euler beam (ksp = 0)');
    xlabel('Time (sec)')
    ylabel('Deflection (m)')
    hold on
    j=j+1;
end
end
% figure
%end
r=array2table(s,'variableNames',{'temp','temp1','alpha','beta','delstat','cd','ksp','P','z1','z2','k1','const','etai0'});
r
% response1real
%as
%asr
asi

Answer 1

[per isakson]

I put your code in a function and profiled it. See profile, Profile execution time for functions. I let it run a couple of minutes and clicked Pause and Quit Debugging.

I observed two issues.

1. Some variables needs Preallocation. (Hover over the orange bars in two the right. They have something to say!)

.

2. The while loop converge very slowly - if at all.

[Abhijeet Chodankar]

It didn't fix the issue. I have implemented the pre-allocation principle, still, it takes a lot of time with no output. Is there a way to make this program execute faster? Did the program execute completely for you? How can I fix the while loop issue?

[per isakson]

No, I clicked Pause and Quit Debugging (see above). It is trapped in the while-loop.

Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
Inf, 0.007620
for ever

• See Debug a MATLAB Program
• See Examine Values While Debugging

[Abhijeet Chodankar]

Thank you very much!!!! This is a very important lead.