Solving an Equation and then plotting a graph

Hi everyone,

Trying to solve an equation and plot a graph but its going wrong and I am not good enought at MATLAB to fix it.

The procedure:

1) epscm running from 0 to 100 (*1E-3).

2) the connection between epss to epscm depends on the varaible of c.

3)compression = f1(c)

4)tension = f2(c) =sigmasteel(c)*As1/1000

5) find the c out of f1(c)=f2(c) for every epscm running from 0 to 100 (*1E-3).

6) with the c i found for evey epscm i calculate phi(i) and M(i)

7) plotting a graph of phi,M

Thank you very much.

b=400;  %mm
d=500;  %mm
fck1=30;  %Mpa
Ecshah=57000/145*(fck1*145)^0.5;  %Mpa
Es=200000;  %Mpa
Esh=8500; %Mpa
As1=3000;  %mm^2
fy=500;  %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck1*145+0.00195;
kshah=0.025*fck1*10^3;
A=Ecshah*eps0/fck1;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;

epscmv = linspace(0.1, 100, 5000)*1E-3;
for i=1:numel(epscmv);
    epscm = epscmv(i);
    epss=@(c) (d-c)/c*epscm;
    funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);    
    compression=@(c) b*fck1*c/epscm*integral(funCshah,0,epscm)/1000;
    sigmaSteel=@(c) Es*epss .* (epss<=epsy) + fy .* (epss>epsy & epss<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss)./(epssu-epssh)).^(1/P)) .* (epss>epssh & epss<=epssu) + 0 .* (epss>epssu);
    tension=@(c) sigmaSteel.*As1/1000;
    c(i)=fsolve(@(c) compression(c)-tension(c),1);
    funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0); 
    M(i)=b*fck1*c(i)/epscm*integral(funM,0,epscm)/1000000;
    phi(i)=epscm/c(i);
end

[Mmax,idx]=max(M) %[kNm]
phiAtMmax=phi(idx)  %[1/mm]
epsAtMmax=epscmv(idx)*1000  %[promil]


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

Answer 1

[Star Strider]

In these two lines:

sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
tension=@(c) sigmaSteel(c).*As1/1000;

note that ‘epss’ and ‘sigmaSteel’ respectively are functions, so it is necessary to evaluate them as functions in their respective anonymous functions. (I corrected those lines, so simply copy them and paste them to your code in place of the present syntax.)

This ran without error with these changes, and the plot worked.

I am not exactly sure what you are doing. However if you are creating the anonymous functions in each iteration of the loop simpply to use the current value of ‘espcm’, a more efficient solution would be to create them before the loop, add ‘espcm’ to the argument list of each function that uses it, then evaluate the existing, pre-defined functions inside the loop.

P.S. — Thank you again for the email alerting me to your Question.

[Shimon Katzman]

Hi Star,

1) Can you please attach here what did it plot to you? because it didnt work for me .. :(

b=400;  %mm
d=500;  %mm
fck1=30;  %Mpa
Ecshah=57000/145*(fck1*145)^0.5;  %Mpa
Es=200000;  %Mpa
Esh=8500; %Mpa
As1=3000;  %mm^2
fy=500;  %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck1*145+0.00195;
kshah=0.025*fck1*10^3;
A=Ecshah*eps0/fck1;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;

epscmv = linspace(0.1, 20, 5000)*1E-3;
for i=1:numel(epscmv);
    epscm = epscmv(i);
    epss=@(c) (d-c)/c*epscm;
    funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);    
    compression=@(c) b*fck1*c/epscm*integral(funCshah,0,epscm)/1000;
    sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
    tension=@(c) sigmaSteel(c).*As1/1000;
    c(i)=fsolve(@(c) compression(c)-tension(c),1);
    funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0); 
    M(i)=b*fck1*c(i)/epscm*integral(funM,0,epscm)/1000000;
    phi(i)=epscm/c(i);
end

[Mmax,idx]=max(M) %[kNm]

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

2) What is the best way to create a single plot with 4 graphs with different fck,As parameters:

for example graph1: fck=30,As=3000 graph 2: fck=30,As=5000, graph 3: fck=90,As=3000, graph 4: fck=90, As=5000.

Thank you very much.

[Star Strider]

1) This runs for me without errors:

b=400;  %mm
d=500;  %mm
fck1=30;  %Mpa
Ecshah=57000/145*(fck1*145)^0.5;  %Mpa
Es=200000;  %Mpa
Esh=8500; %Mpa
As1=3000;  %mm^2
fy=500;  %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck1*145+0.00195;
kshah=0.025*fck1*10^3;
A=Ecshah*eps0/fck1;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
epscmv = linspace(0.1, 100, 5000)*1E-3;
for i=1:numel(epscmv);
    epscm = epscmv(i);
    epss=@(c) (d-c)/c*epscm;
    funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);    
    compression=@(c) b*fck1*c/epscm*integral(funCshah,0,epscm)/1000;
    sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
    tension=@(c) sigmaSteel(c).*As1/1000;
    c(i)=fsolve(@(c) compression(c)-tension(c),1);
    funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0); 
    M(i)=b*fck1*c(i)/epscm*integral(funM,0,epscm)/1000000;
    phi(i)=epscm/c(i);
end
[Mmax,idx]=max(M) %[kNm]
phiAtMmax=phi(idx)  %[1/mm]
epsAtMmax=epscmv(idx)*1000  %[promil]
plot(phi(1:idx), M(1:idx))
grid on
xlabel('phi  [1/mm]')
ylabel('Moment  [kNm]')

It just takes a while.

2) This slight modification of the other code also runs for me without errors, and produces four plots on one set of axes:

b=400;  %mm
d=500;  %mm
% fck1=30;  %Mpa
Asv = [3000 5000 3000 5000];
fckv = [30  30  90  90];
for k = 1:numel(fckv)
    fck1 = fckv(k)
    Ecshah=57000/145*(fck1*145)^0.5;  %Mpa
    Es=200000;  %Mpa
    Esh=8500; %Mpa
    % As1=3000;  %mm^2
    fy=500;  %Mpa
    fsu=750; %Mpa
    epssh=0.009;
    epssu=0.075;
    eps0=1.027*10^-7*fck1*145+0.00195;
    kshah=0.025*fck1*10^3;
    A=Ecshah*eps0/fck1;
    P=Esh*((epssu-epssh)/(fsu-fy));
    epsy=fy/Es;
    epscmv = linspace(0.1, 100, 5000)*1E-3;
    As1 = Asv(k);
    for i=1:numel(epscmv);
        epscm = epscmv(i);
        epss=@(c) (d-c)/c*epscm;
        funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
        compression=@(c) b*fck1*c/epscm*integral(funCshah,0,epscm)/1000;
        sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
        tension=@(c) sigmaSteel(c).*As1/1000;
        c(i)=fsolve(@(c) compression(c)-tension(c),1);
        funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
        M(i,k)=b*fck1*c(i)/epscm*integral(funM,0,epscm)/1000000;
        phi(i,k)=epscm/c(i);
    end
end
[Mmax,idx]=max(M) %[kNm]
phiAtMmax=phi(idx)  %[1/mm]
epsAtMmax=epscmv(idx)*1000  %[promil]

figure
hold all
for k = 1:numel(fckv)
    plot(phi(1:idx,k), M(1:idx,k))
    lgdstr{k} = sprintf('fck = %2d, As = %4d',fckv(k), Asv(k));
end
hold off
grid on
xlabel('phi  [1/mm]')
ylabel('Moment  [kNm]')
legend(lgdstr)

It takes even longer because of the loops, and eventually produces this plot:

[Shimon Katzman]

Wow, once again Thank you. you are the best!

Something is wrong with the graph, it should be somethink close to this:

I will try to figure out what is wrong and if i will struggle i will ask for a help again.

Thank you Star.

[Star Strider]

I very much appreciate your compliment!

As always, my pleasure!

[Shimon Katzman]

Hi Star,

Can't figure out what is wrong.

isn't it wierd that i get the Mmax at idx=3?

the graph of the epcm,sigmaSteel is wrong too.. if you remember it should be similiar to the graph here:

https://www.mathworks.com/matlabcentral/answers/494042-error-nonscalar-arrays-of-function-handles-are-not-allowed-use-cell-arrays-instead?s_tid=prof_contriblnk

b=400;  %mm
d=500;  %mm
fck1=30;  %Mpa
Ecshah=57000/145*(fck1*145)^0.5;  %Mpa
Es=200000;  %Mpa
Esh=8500; %Mpa
As1=3000;  %mm^2
fy=500;  %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck1*145+0.00195;
kshah=0.025*fck1*10^3;
A=Ecshah*eps0/fck1;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
epscmv = linspace(0.1, 100, 5000)*1E-3;
for i=1:numel(epscmv);
    epscm = epscmv(i);
    epss=@(c) (d-c)/c*epscm;
    funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);    
    compression=@(c) b*fck1*c/epscm*integral(funCshah,0,epscm)/1000;
    sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
    tension=@(c) sigmaSteel(c).*As1/1000;
    c(i)=fsolve(@(c) compression(c)-tension(c),1);
    funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0); 
    M(i)=b*fck1*c(i)/epscm*integral(funM,0,epscm)/1000000;
    phi(i)=epscm/c(i);
    epss1(i)= (d-c(i))/c(i)*epscm;
    sigmaSteel1(i)= Es*epss1(i) .* (epss1(i)<=epsy) + fy .* (epss1(i)>epsy & epss1(i)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss1(i))./(epssu-epssh)).^(1/P)) .* (epss1(i)>epssh & epss1(i)<=epssu) + 0 .* (epss1(i)>epssu);
end

[Mmax,idx]=max(M) %[kNm]
phiAtMmax=phi(idx)  %[1/mm]
epsAtMmax=epscmv(idx)*1000  %[promil]


plot(epscmv(1:idx), sigmaSteel1(1:idx))
grid on
xlabel('epsilon cm')
ylabel('SigmaSteel [Mpa]')

Thank you very much.

[Shimon Katzman]

also, the graph of (phi,M) does not suppose to be so different from the graph from this post:

(the only "critical" thing that changed is the "tension" function).

https://www.mathworks.com/matlabcentral/answers/492893-error-vectors-must-be-the-same-length?s_tid=prof_contriblnk

(of course, plotting phi,m and not epscm,cdivd)

[Star Strider]

As always, my pleasure!

It took a few minutes to experiment with this. One problem is that the fsolve call is highly dependent on the initial parameter estimate (as are all nonlinear solvers). Giving it a much larger initial estimate (I experimented with 100 and 1000) begins to produce the result you want:

c(i)=fsolve(@(c) compression(c)-tension(c),1000);

with ‘idx’ being 135 with this (1000) initial estimate.

I also experimented with the plot calls to see what the complete ‘sigmaSteel1’ result looked like:

figure
plot(epscmv(1:idx), sigmaSteel1(1:idx))
grid on
xlabel('epsilon cm')
ylabel('SigmaSteel [Mpa]')

figure
plot(epscmv, sigmaSteel1)
grid on
xlabel('epsilon cm')
ylabel('SigmaSteel [Mpa]')

figure
semilogy(epscmv, sigmaSteel1)
% plot(epscmv(1:idx), sigmaSteel1(1:idx))
grid on
xlabel('epsilon cm')
ylabel('SigmaSteel [Mpa]')

figure
loglog(epscmv, sigmaSteel1)
% plot(epscmv(1:idx), sigmaSteel1(1:idx))
grid on
xlabel('epsilon cm')
ylabel('SigmaSteel [Mpa]')

I also looked at your code again (specifically ’P’ and ‘SigmaSteel1’) and found no errors, comparing it to ‘1.jpg’. The plot still does not look exactly like the one in ‘1.jpg’, however it is much closer.

[Shimon Katzman]

Star YOU ARE THE B E S T!

Thank you so much.

just a small question, how do i continute to plot the graph just a litlle bit after the max, like the picture:

[Star Strider]

I very much appreciate your compliment!

The easiest way is to add some arbitrary number to ‘idx’ in the plot calls, for example:

plot(phi(1:idx+10,k), M(1:idx+10,k))

I chose 10 here, however anything less that sums to less than or equal to the column size of the variables will work.

[Shimon Katzman]

Thank you very much!

[Star Strider]

As always, my pleasure!

[Shimon Katzman]

Is it possible to continue the graph untill 95% of Mmax (after the peak)?

[Star Strider]

It took a while to reply because I need to test everything I post and the code takes a while to run.

It is possible, however it does not appear to make much difference in the plot. Posting only the relevant portion of the code:

[Mmax,idx]=max(M) %[kNm]
idx95 = find(M(idx:end) < 0.95*Mmax, 1, 'first')        % Find Index Beyond ‘idx’ Where Curve Is < 95% Of Max
idxext = idx + idx95                                    % New (Extended) Index
phiAtMmax=phi(idx)  %[1/mm]
epsAtMmax=epscmv(idx)*1000  %[promil]

figure
plot(epscmv(1:idxext), sigmaSteel1(1:idxext))
grid on
xlabel('epsilon cm')
ylabel('SigmaSteel [Mpa]')

figure
plot(epscmv, M, '-b',    epscmv(idxext), M(idxext), 'r+')
grid

The second figure (here) plotting ‘M’, demonstrates that the code does actually test for and return the correct value. (The ‘idxext’ value is 256.)

[Shimon Katzman]

Wow. Thank you Star, like always, you helped me so much.

[Star Strider]

As always, my pleasure!

[Shimon Katzman]

Hi Star,

How can i find the value of M,c,epscmv when epss=epssh (=0.009) and to mark it on the graph?

Thank You.

[Shimon Katzman]

Also i need to mark the Mmax.

i tryed this:

for k = 1:numel(fckv)
    plot(phi(1:idx+50,k), M(1:idx+50,k))
    lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
    plot(phi(idx,k),M(idx,k),'r*')  % marking the Mmax
end

But it plotted this:

(The legend is wrong too :( )

[Star Strider]

As always, my pleasure.

I am not quite sure what you want.

Try this (at the end of the present code):

epss_c = fsolve(@(c) epss(c)-epssh, 100)
c_idx = find(c >= epss_c, 1, 'first')
epscm_exact = interp1(c([-5 5]+c_idx), epscmv([-5 5]+c_idx), epss_c)

figure
plot(epscmv, sigmaSteel1, '-k',  epscmv, M, '-b',    epscmv, c, '-r')
hold on
plot(epscmv(c_idx), sigmaSteel1(c_idx), '+k',  epscmv(c_idx), M(c_idx), '+b',    epscmv(c_idx), c(c_idx), '+r')
hold off
grid
legend('sigmaSteel1', 'M', 'c', 'epss = 0.009', 'Location','E')

figure
plot(epscmv, sigmaSteel1, '-k',  epscmv, M, '-b',    epscmv, c, '-r', epscmv, epss_c*ones(size(epscmv)), '-g')
grid
legend('sigmaSteel1', 'M', 'c', 'epss = 0.009', 'Location','E')

The fsolve call calculates the value of ‘epss’ at ‘epssh’, the find call returns the closest index, and ‘epscm_exact’ is the interpolated value of ‘epscm’ at ‘epss_c’.

Edit the plot calls to get the result you want.

I also timed it, out of curiosity. It takes 2½ minutes to run on my desktop:

Elapsed time is 150.782281 seconds.

[Star Strider]

Shimon Katzman’s Answer moved here —

Thank you very much, i will try it.

Another question.

i am trying to plot a epscmv,sigmaSteel1 graph but it Errors me :((((

b=400;  %mm
d=500;  %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30  30  90  90]; %Mpa

for k = 1:numel(fckv)
    fck = fckv(k);
    Ecshah=57000/145*(fck*145)^0.5;  %Mpa
    Es=200000;  %Mpa
    Esh=8500; %Mpa
    fy=500;  %Mpa
    fsu=750; %Mpa
    epssh=0.009;
    epssu=0.075;
    eps0=1.027*10^-7*fck*145+0.00195;
    kshah=0.025*fck*10^3;
    A=Ecshah*eps0/fck;
    P=Esh*((epssu-epssh)/(fsu-fy));
    epsy=fy/Es;
    epscmv = linspace(0.1, 100, 5000)*1E-3;
    As = Asv(k);
    for i=1:numel(epscmv);
        epscm = epscmv(i);
        epss=@(c) (d-c)/c*epscm;
        funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
        compression=@(c) b*fck*c/epscm*integral(funCshah,0,epscm)/1000;
        sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
        tension=@(c) sigmaSteel(c).*As/1000;
        c(i)=fsolve(@(c) compression(c)-tension(c),1000);
        funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
        M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
        phi(i,k)=epscm/c(i);
        epss1(i)= (d-c(i))/c(i)*epscm;
        sigmaSteel1(i)= Es*epss1(i) .* (epss1(i)<=epsy) + fy .* (epss1(i)>epsy & epss1(i)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss1(i))./(epssu-epssh)).^(1/P)) .* (epss1(i)>epssh & epss1(i)<=epssu) + 0 .* (epss1(i)>epssu);
    end
end
[Mmax,idx]=max(M) %[kNm]
phiAtMmax=phi(idx)  %[1/mm]
epsAtMmax=epscmv(idx)*1000  %[promil]

figure
hold all
for k = 1:numel(fckv)
    plot(epscmv(1:idx,k), sigmaSteel1(1:idx,k))
    lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));

end
hold off
grid on
xlabel('epsilon cm')
ylabel('Sigma Steel [Mpa]')
legend(lgdstr)

returns:

Index exceeds matrix dimensions.

Error in advencedconcrete4NEW (line 76)
    plot(epscmv(1:idx,k), sigmaSteel1(1:idx,k))

[Star Strider]

When I time the main part of this code, the result is:

Elapsed time is 669.804737 seconds.

So no immediate miracles!

That aside, tthis runs without error:

tic
b=400;  %mm
d=500;  %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30  30  90  90]; %Mpa
for k = 1:numel(fckv)
    fck = fckv(k);
    Ecshah=57000/145*(fck*145)^0.5;  %Mpa
    Es=200000;  %Mpa
    Esh=8500; %Mpa
    fy=500;  %Mpa
    fsu=750; %Mpa
    epssh=0.009;
    epssu=0.075;
    eps0=1.027*10^-7*fck*145+0.00195;
    kshah=0.025*fck*10^3;
    A=Ecshah*eps0/fck;
    P=Esh*((epssu-epssh)/(fsu-fy));
    epsy=fy/Es;
    epscmv = linspace(0.1, 100, 5000)*1E-3;
    As = Asv(k);
    for i=1:numel(epscmv);
        epscm = epscmv(i);
        epss=@(c) (d-c)/c*epscm;
        funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
        compression=@(c) b*fck*c/epscm*integral(funCshah,0,epscm)/1000;
        sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
        tension=@(c) sigmaSteel(c).*As/1000;
        c(i)=fsolve(@(c) compression(c)-tension(c),1000);
        funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
        M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
        phi(i,k)=epscm/c(i);
        epss1(k,i)= (d-c(i))/c(i)*epscm;
        sigmaSteel1(k,i)= Es*epss1(k,i) .* (epss1(i)<=epsy) + fy .* (epss1(k,i)>epsy & epss1(k,i)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss1(k,i))./(epssu-epssh)).^(1/P)) .* (epss1(k,i)>epssh & epss1(k,i)<=epssu) + 0 .* (epss1(k,i)>epssu);
    end
    [Mmax(k),idx(k)]=max(M(:,k)) %[kNm]
    phiAtMmax(k)=phi(idx(k))  %[1/mm]
    epsAtMmax(k)=epscmv(idx(k))*1000  %[promil]
end
toc

figure
hold all
for k = 1:numel(fckv)
    plot(epscmv(1:idx(k)), sigmaSteel1(k,1:idx(k)))
    lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
end
hold off
grid on
xlabel('epsilon cm')
ylabel('Sigma Steel [Mpa]')
legend(lgdstr, 'Location','NW')

Here, ‘epss1’ and ‘sigmaSteel1’ are now both (4x5000) matrices, and ‘M’ a (5000x4) matrix. It would likely be worthwhile to preallocate them, however I would not expect a significant increase in speed, since the fsolve call takes most of the time. The matrix formulations may be useful if you want to emulate the code for determining the same indices as in my previous Comment.

I have no idea where the oscillations come from in the third iteration (90, 3000). Using:

idx(k) = find(diff(M(:,k) < 0, 1, 'first')
Mmax(k) = M(idx,k)

might be an acceptable alternative to the max function if you want to define the ‘Mmax’ value as the last value before the oscillations begin. This will also find ‘Mmax’ correctly for all the values of ‘M’ that do not oscillate (assuming the oscillation is due to oscillations in ‘M’, since I did not explore that).

[Shimon Katzman]

Star thank you so mucj, i will try it :)

[Star Strider]

As always, my pleasure!

[Shimon Katzman]

Hi Star i changed the max function, but now it errors me :(

tic
b=400;  %mm
d=500;  %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30  30  90  90]; %Mpa
for k = 1:numel(fckv)
    fck = fckv(k);
    Ecshah=57000/145*(fck*145)^0.5;  %Mpa
    Es=200000;  %Mpa
    Esh=8500; %Mpa
    fy=500;  %Mpa
    fsu=750; %Mpa
    epssh=0.009;
    epssu=0.075;
    eps0=1.027*10^-7*fck*145+0.00195;
    kshah=0.025*fck*10^3;
    A=Ecshah*eps0/fck;
    P=Esh*((epssu-epssh)/(fsu-fy));
    epsy=fy/Es;
    epscmv = linspace(0.1, 100, 5000)*1E-3;
    As = Asv(k);
    for i=1:numel(epscmv);
        epscm = epscmv(i);
        epss=@(c) (d-c)/c*epscm;
        funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
        compression=@(c) b*fck*c/epscm*integral(funCshah,0,epscm)/1000;
        sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
        tension=@(c) sigmaSteel(c).*As/1000;
        c(i)=fsolve(@(c) compression(c)-tension(c),1000);
        funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
        M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
        phi(i,k)=epscm/c(i);
        epss1(k,i)= (d-c(i))/c(i)*epscm;
        sigmaSteel1(k,i)= Es*epss1(k,i) .* (epss1(i)<=epsy) + fy .* (epss1(k,i)>epsy & epss1(k,i)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss1(k,i))./(epssu-epssh)).^(1/P)) .* (epss1(k,i)>epssh & epss1(k,i)<=epssu) + 0 .* (epss1(k,i)>epssu);
    end
    idx(k) = find(diff(M(:,k) < 0, 1, 'first'))%[kNm]
    Mmax(k) = M(idx,k) %[kNm]
    phiAtMmax(k)=phi(idx(k))  %[1/mm]
    epsAtMmax(k)=epscmv(idx(k))*1000  %[promil]
end
toc

figure
hold all
for k = 1:numel(fckv)
    plot(epscmv(1:idx(k)), sigmaSteel1(k,1:idx(k)))
    lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
end
hold off
grid on
xlabel('epsilon cm')
ylabel('Sigma Steel [Mpa]')
legend(lgdstr, 'Location','NW')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Error using diff
Dimension argument must be a positive integer scalar within indexing range.

Error in advencedconcrete4SigmaSteel (line 36)
    idx(k) = find(diff(M(:,k) < 0, 1, 'first'))%[kNm]
    
    

Also, the graphs should look similar to this, and before the oscillation its very different..

[Shimon Katzman]

Uh, i understood what was wrong.

the k was missing in (epss1(k,i)<=epsy) :

sigmaSteel1(k,i)= Es*epss1(k,i) .* (epss1(k,i)<=epsy) + fy .* (epss1(k,i)>epsy & epss1(k,i)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss1(k,i))./(epssu-epssh)).^(1/P)) .* (epss1(k,i)>epssh & epss1(k,i)<=epssu) + 0 .* (epss1(k,i)>epssu);

Thank You.

[Star Strider]

As always, my pleasure!

I am running your code now (it takes a while), and have drafted this Comment that I was waiting to post after I tested my latest update:

—————

I misplaced a parenthesis. (I copied this from some test code to be sure it worked, and forgot about ‘M’ needing subscripts.)

Those assignments should be:

idx(k) = find(diff(M(:,k)) < 0, 1, 'first')%[kNm]
Mmax(k) = M(idx(k),k) %[kNm]

With those changes the code runs:

Elapsed time is 645.655820 seconds.

The only problem is that ‘M’ apparently does not oscillate (or is not the source of the oscillations), so the changed code would have to be used on ‘sigmaSteel12’ to prevent the oscillations, although this could be done after the loop.

To detect the oscillations and end with the maximum value that is not after the oscillations begin, run this after the main loops:

for k = 1:numel(fckv)
    sS1_idx(k) = find(diff(sigmaSteel1(k,:)) <0, 1, 'first');
end

figure
hold all
for k = 1:numel(fckv)
    plot(epscmv(1:sS1_idx(k)), sigmaSteel1(k,1:idx(k)))
    lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
end
hold off
grid on
xlabel('epsilon cm')
ylabel('Sigma Steel [Mpa]')
legend(lgdstr, 'Location','eastoutside')

—————

However, if you have found the problem, those steps may not be necessary.

Your code just finished:

Elapsed time is 652.802339 seconds

and with those changes, runs without error.

I have not included the new code you posted. I will run it with that tomorrow, so I can see the result.