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

Question

Shimon Katzman 2019-12-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494721-marking-a-specific-point-in-a-graph-made-out-of-a-loop

关闭： MATLAB Answer Bot 2021-8-20

在 MATLAB Online 中打开

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)

2 个评论
显示无隐藏无

Shimon Katzman 2019-12-4

Can anyone help please please?

Star Strider 2019-12-4

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 2019-12-4

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494721-marking-a-specific-point-in-a-graph-made-out-of-a-loop#answer_404804

在 MATLAB Online 中打开

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2019-12-5

编辑：Star Strider 2019-12-5

在 MATLAB Online 中打开

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)

1Marking a specific point in a graph made out of a loop - 2019 12 05.png

Shimon Katzman 2019-12-5

在 MATLAB Online 中打开

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')