Solving an Equation and then plotting a graph

조회 수: 5 (최근 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
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 CenterFile Exchange에서 Logical에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by