MATLAB: Solving an Equation and then plotting a graph

functionhelpMATLAB

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

Best Answer

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.
Related Question