Solving an Equation and then plotting a graph

2 次查看（过去 30 天）
Shimon Katzman
Shimon Katzman 2019-12-1
评论： Star Strider 2019-12-4
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
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 个评论
显示 23更早的评论
Shimon Katzman
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
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 中查找有关 Graphics Object Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by