if true
syms P
l=100; a=50; b=20; h0=50; E=200; h1=50;
for n=1:0.02:4
I0=b*h0^3/12; I1=b*h1^3/12;
L0=sqrt(P/(E*I0-P*n));
L1=sqrt(P/(E*I1-P*n));
v=0.3;
s=0.5;
f=1.93-3.07*s+14.53*s^2-25.11*s^3+25.8*s^4;
k=(E*I1)/(6*pi*h1*f*(1-v^2));
a11=cos(L1*l) - cos(L1*a) - L1*a*sin(L1*l) + L1*l*sin(L1*l);
a12=sin(L1*l) - sin(L1*a) + L1*a*cos(L1*l) - L1*l*cos(L1*l);
a13=cos(L0*a) - 1;
a14=sin(L0*a) - L0*a;
a21=-(L1*sin(L1*a) - L1*sin(L1*l)) - L1^2*cos(L1*a)*((P*n - E*I1)/k);
a22=(L1*cos(L1*a) - L1*cos(L1*l)) - L1^2*sin(L1*a)*((P*n - E*I1)/k);
a23=L0*sin(L0*a);
a24=(L0 - L0*cos(L0*a));
a31=-L1^2*cos(L1*a)*(P*n - E*I1);
a32=-L1^2*sin(L1*a)*(P*n - E*I1);
a33=L0^2*cos(L0*a)*(P*n - E*I0);
a34=L0^2*sin(L0*a)*(P*n - E*I0);
a41=L1^3*sin(L1*a)*(P*n - E*I1);
a42=-L1^3*cos(L1*a)*(P*n - E*I1);
a43=-L0^3*sin(L0*a)*(P*n - E*I0);
a44=L0^3*cos(L0*a)*(P*n - E*I0);
A=[a11 a12 a13 a14;a21 a22 a23 a24;a31 a32 a33 a34;a41 a42 a43 a44];
delta1=det(A);
delta2=subs(delta1);
delta=simplify(delta2);
F=matlabFunction(delta);
P1=[];
for j = 0.01:10
PP = fzero(F, j);
P1 = [P1; PP];
end
P_cr = min(P1(P1 > 0));
D=plot(n,P_cr,'k.','linewidth',1.5);
DD=double(D);
hold on
end
end
Best Answer