Solving an Equation and then plotting a graph

1 ビュー (過去 30 日間)
Shimon Katzman 2019 年 12 月 1 日
Commented: 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]')

0 件のコメント

サインイン to comment.

採用された回答

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 月 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.
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..
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.
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.

サインイン to comment.