Stopping a loop on the peak

Hi

I want to stop the loop when i get to the peak of the M

b=300;  %mm
d=400;  %mm
fc=40;  %Mpa
Ecm=22*(fc/10)^0.3*10^3;  %Mpa
Es=200000;  %Mpa
As=2400;  %mm^2
fy=400;  %Mpa
epsc1=min(2.8/1000,0.7*fc^0.31/1000);
epscu=3.5/1000;
k=1.05*Ecm*epsc1/fc;

epscmv = linspace(0.05, 3.5, 50)*1E-3;
for i=1:50(epscmv);
    epscm = epscmv(i);
    funC=@(epsc) (k*epsc/epsc1-(epsc/epsc1).^2)./(1+(k-2)*epsc/epsc1);
    compression=@(c) b*fc*c/epscm*integral(funC,0,epscm)/1000;
    tension=@(c) min(Es*(d-c)/c*epscm*As/1000,fy*As/1000);
    c(i)=fsolve(@(c) compression(c)-tension(c),1);
    funM=@(epsc)(k*epsc./epsc1-(epsc./epsc1).^2)./(1+(k-2)*epsc./epsc1).*(d-c(i)+(c(i)./epscm).*epsc); 
    M(i)=b*fc*c(i)/epscm*integral(funM,0,epscm)/1000000;
    phi(i)=epscm/c(i);
    if M(i)<M(i-1)
        break;
    end
end

[Mmax,idx]=max(M)
phiAtMmax=phi(idx)

plot(phi, M)
grid on
xlabel('phi  [1/mm]')
ylabel('Moment  [kNm]')

A kind forum member helped me but now I get the Error:

Attempted to access M(0); index must be a positive integer or logical.

Error in advencedconcrete32a (line 22)
    if M(i)<M(i-1)

Answer 1

[Star Strider]

Change the if block to:

    if (i > 1) & (M(i)<M(i-1))
        break;
    end

and your code appears to give the desired result.

This avoids the ‘zero index’ problem by comparing only the second and later values of ‘M’.

[Shimon Katzman]

thank you star, but it doesnt stop the loop at the maximum like you can see in the picture i attached:

the maximum should be at 2.5 phi. why is that?

the code:

b=300;  %mm
d=400;  %mm
fc=40;  %Mpa
Ecm=22*(fc/10)^0.3*10^3;  %Mpa
Es=200000;  %Mpa
As=2400;  %mm^2
fy=400;  %Mpa
epsc1=min(2.8/1000,0.7*fc^0.31/1000);
epscu=3.5/1000;
k=1.05*Ecm*epsc1/fc;

epscmv = linspace(0.05, 3.5, 50)*1E-3;
for i=1:50(epscmv);
    epscm = epscmv(i);
    funC=@(epsc) (k*epsc/epsc1-(epsc/epsc1).^2)./(1+(k-2)*epsc/epsc1);
    compression=@(c) b*fc*c/epscm*integral(funC,0,epscm)/1000;
    tension=@(c) min(Es*(d-c)/c*epscm*As/1000,fy*As/1000);
    c(i)=fsolve(@(c) compression(c)-tension(c),1);
    funM=@(epsc)(k*epsc./epsc1-(epsc./epsc1).^2)./(1+(k-2)*epsc./epsc1).*(d-c(i)+(c(i)./epscm).*epsc); 
    M(i)=b*fc*c(i)/epscm*integral(funM,0,epscm)/1000000;
    phi(i)=epscm/c(i);
    if (i > 1) & (M(i)<M(i-1))
        break;
    end
end

[Mmax,idx]=max(M)
phiAtMmax=phi(idx)

plot(phi, M)
grid on
xlabel('phi  [1/mm]')
ylabel('Moment  [kNm]')

[Star Strider]

As always, my pleasure!

One problem is that this:

for i=1:50(epscmv);

should be replaced by:

for i=1:numel(epscmv)

With that change and increasing the resolution of ‘epscmv’ (by changing the third argument to linspace):

epscmv = linspace(0.05, 3.5, 250)*1E-3;

gives these results:

Mmax =
        342.5

idx =
   197

phiAtMmax =
   2.5061e-05

Using the original definition of ‘epscmv’ (with only 50 elements):

epscmv = linspace(0.05, 3.5, 50)*1E-3;

the results are:

Mmax =
       342.49

idx =
    40

phiAtMmax =
   2.5409e-05

Increasing the number of elements of linspace to 500 produces:

phiAtMmax =
   2.5078e-05

with the value of ‘M’ essentially unchanged from the other estimates.

It will become even more accurate (up to the limits of floating-point calculation precision) with increased resolution.

[Shimon Katzman]

That is GREAT!!!

THANK YOU SO MUCH!

[Star Strider]

As always, my pleasure!