Solving an Equation and then plotting a graph
2 次查看(过去 30 天)
显示 更早的评论
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]')
采纳的回答
Star Strider
2019-12-1
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.
25 个评论
Shimon Katzman
2019-12-1
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
2019-12-1
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:
![1Solving an Equation and then plotting a graph m- 2019 12 01.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/251785/1Solving%20an%20Equation%20and%20then%20plotting%20a%20graph%20m-%202019%2012%2001.png)
Shimon Katzman
2019-12-2
Wow, once again Thank you. you are the best!
Something is wrong with the graph, it should be somethink close to this:
![1.jpg](https://www.mathworks.com/matlabcentral/answers/uploaded_files/251803/1.jpeg)
I will try to figure out what is wrong and if i will struggle i will ask for a help again.
Thank you Star.
Shimon Katzman
2019-12-2
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:
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
2019-12-2
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).
(of course, plotting phi,m and not epscm,cdivd)
Star Strider
2019-12-2
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
2019-12-2
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:
![2.jpg](https://www.mathworks.com/matlabcentral/answers/uploaded_files/252019/2.jpeg)
Star Strider
2019-12-2
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.
Star Strider
2019-12-2
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
2019-12-3
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
2019-12-3
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 :( )
![3.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/252458/3.png)
Star Strider
2019-12-3
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
2019-12-3
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
2019-12-3
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
2019-12-4
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..
![4.jpg](https://www.mathworks.com/matlabcentral/answers/uploaded_files/252521/4.jpeg)
Shimon Katzman
2019-12-4
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
2019-12-4
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)