Marking a specific point in a graph made out of a loop

Hi everyone.

Trying to mark 2 points at the 4 curves in the graph:

1) Marking the peak in each graph (the Mmax point)

2) Marking the point where the value of M belongs to epss=epssh in each graph.

Thank You very much.

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);
    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+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
hold off
grid on
xlabel('phi  [1/mm]')
ylabel('Moment  [kNm]')
legend(lgdstr)

[Shimon Katzman]

Can anyone help please please?

[Star Strider]

Working on it. Windows 10 gave me a BSOD a few minutes ago for the second time in as many days (I have no kind works for Micro$oft just now), and it took a while to put the pieces of everything back together.

Answer 1

[Star Strider]

Try this:

b=400;  %mm
d=500;  %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30  30  90  90]; %Mpa
tic
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);
    end
    idx(k) = find(diff(M(:,k)) < 0, 1, 'first')%[kNm]
    Mmax(k) = M(idx(k),k) %[kNm]
    phiAtMmax(k)=phi(idx(k))  %[1/mm]
    epsAtMmax(k)=epscmv(idx(k))*1000  %[promil]
end
toc
% [Mmax,idx]=max(M) %[kNm]
% phiAtMmax=phi(idx)  %[1/mm]
% epsAtMmax=epscmv(idx)*1000  %[promil]

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
hold all
for k = 1:numel(fckv)
    hp(k) = plot(phi(1:idx+50,k), M(1:idx+50,k));
%     plot([1 1]*epss_c, ylim, ':k')
    lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
    plot(phi(idx(k),k),M(idx(k),k),'r*')  % marking the Mmax
end
hold off
grid on
xlabel('phi  [1/mm]')
ylabel('Moment  [kNm]')
hl = legend(hp,lgdstr);
hl.FontSize = 7;
% text(epss_c, 400, sprintf('\\leftarrow epss_c = %.1f', epss_c), 'HorizontalAlignment','left')

This correctly plots the maxima of every curve (although the green curve maximum is beyond the green curve plotted limit (30 Mpa, 3000 mm^3), and I corrected the legend problem as well, so that it only refers to the lines, not the markers.

I do not understand about ‘epss_c’ with respect to every curve. It appears to be beyond every curve. In any event, I have no idea how to calculate it with respect to plotting it as a function of ‘phi’.

The Windows 10 crash (that I did not immediately recognise) put me abouit two hours behind.

[Shimon Katzman]

Hi Star and thank you very very much for your help.

I cant understand why the peak of the yellow and the purple graph is beyond the graph.

The code suppose to stop the graph when it gets to the maximum of M..

phi values are wrong too:

phiAtMmax =

   1.0e-04 *

    0.1586    0.1928    0.2454    0.2233
    
    

Do you know why?

Thank you once again

[Star Strider]

The find call should return the same value as the maximum. I changed that to see if I could deal with the oscillation. I have no idea why that value and others are at the end or far from the the curve ends.

Since ‘phi’ is a 2D matrix, it is necessary to provide two subscripts to it. (The alternative is to use ‘linear indexing’ that will work, although it may not return the value you want.) I have not plotted ‘phi’, however this may provide the correct results for ‘phiAtMax’:

[Mmax(k),idx(k)]=max(M(:,k)) %[kNm]
phiAtMmax(k)=phi(idx(k),k)  %[1/mm]

Experiment with the ‘phi’ indexing to get the result you want.

In the plot loop, ‘idx’ needs to be subscripted as ‘idx(k)’:

hp(k) = plot(phi(1:idx(k)+50,k), M(1:idx(k)+50,k));

EDIT — (5 Dec 2019 at 13:15)

With these changes:

phiAtMmax =
    15.8590e-006    11.5660e-006    61.5382e-006    35.0650e-006

epsAtMmax =
     2.7779    3.2775     4.2367     3.7970

I have not looked at Walter’s code.

EDIT — (5 Dec 2019 at 13:32)

[Shimon Katzman]

Now it is great!

I am so thankfull to you, can is ask another question please?

Trying to plot a graph of epscmv,c but cant figure out why again it errors me:

"Index exceeds matrix dimensions."

b=400;  %mm
d=500;  %mm
Asv = [3000 5000 3000 5000]; %mm^2

fckv = [30  30  90  90]; %Mpa
num_fckv = numel(fckv);

epscmv = linspace(0.1, 100, 5000)*1E-3;
num_epscvm = numel(epscmv);

fsolve_opts = optimoptions('fsolve','Display','off');

M = zeros(num_epscvm, num_fckv);
phi = zeros(num_epscvm, num_fckv);
idx = zeros(1, num_fckv);
Mmax = zeros(1, num_fckv);
phiAtMmax = zeros(1, num_fckv);
epsAtMmax = zeros(1, num_fckv);

tic
for k = 1:num_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;
    As = Asv(k);

    c0 = 1000;
    c = zeros(1, num_epscvm);
    for i=1:num_epscvm
        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), c0, fsolve_opts);
        c0 = c(i);
        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(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);
    end
    idx(k) = find(diff(M(:,k)) < 0, 1, 'first')
    [Mmax(k),idx(k)]=max(M(:,k)) %[kNm]
    phiAtMmax(k)=phi(idx(k),k)  %[1/mm]
    epsAtMmax(k)=epscmv(idx(k))*1000 %[promil]
end
toc


figure
hold all
hp = zeros(1, num_fckv);
lgdstr = cell(1, num_fckv);
for k = 1:num_fckv
    hp(k) = plot(epscmv(1:idx(k)), c(1:idx(k),k));
%     plot([1 1]*epss_c, ylim, ':k')
    lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
    plot(epscmv(idx(k)),c(idx(k),k),'r*')  % marking the Value when Mmax
end
hold off
grid on
xlabel('epsilon cm')
ylabel('c [mm]')
hl = legend(hp,lgdstr);
hl.FontSize = 7;
% % text(epss_c, 400, sprintf('\\leftarrow epss_c = %.1f', epss_c), 'HorizontalAlignment','left')

Thank you very much.

[Star Strider]

As always, my pleasure!

In your code, ‘c’ is a vector, not a matrix. To create a matrrix for it, the easiest way would be to add:

c_mtx(i,k) = c(i);

just before the end of the ‘i’ loop, then change the plot call to:

hp(k) = plot(epscmv(1:idx(k)), c_mtx(1:idx(k),k));

EDIT — (5 Dec 2019 at 14:03)

That appears to plot correctly. You will need to change the related calls to show the correct xlabel and ylabel strings, and remove or change the second plot call.

[Shimon Katzman]

Thank you very much!!!!

[Star Strider]

My pleasure!

Answer 2

[Walter Roberson]

Repaired code is attached. Some optimizations were made as well.

[Shimon Katzman]

Wow Thank you very much.

Thank you.