fsolve求解12阶行列式方程(超越方程),遇到错误。
Show older comments
E2=76e9;
rho2=10500;
%%%%
E1=160e9;
rho1=2300;
%%%%
h1=2*1e-6;
h2=5*1e-6;
h0=h2-h1;
a=150*1e-6;
b=350*1e-6;
l=500*1e-6;
z=10*1e-6;
%%%%
co=1; %缩小系数
syms omiga;
beta1=(10.15179*(co*omiga)^0.5) %beta1、2、3分别为以下行列式中的参数
beta2=(10.1533*(co*omiga)^0.5)
beta3=(10.15179*(co*omiga)^0.5)
for i=1:2000 %求解超越方程特征根,以下12阶行列式中omiga为待求解的未知数
root(i)=fsolve(@(omiga) det([0,1,0,1,0,0,0,0,0,0,0,0;
beta1,0,beta1,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,sin(beta3*l),cos(beta3*l),sinh(beta3*l),cosh(beta3*l);
0,0,0,0,0,0,0,0,beta3*cos(beta3*l),-beta3*sin(beta3*l),beta3*cosh(beta3*l),beta3*sinh(beta3*l);
sin(beta1*a),cos(beta1*a),sinh(beta1*a),cosh(beta1*a),-sin(beta2*a),-cos(beta2*a),-sinh(beta2*a),-cosh(beta2*a),0,0,0,0;
0,0,0,0,sin(beta2*b),cos(beta2*b),sinh(beta2*b),cosh(beta2*b),-sin(beta3*b),-cos(beta3*b),-sinh(beta3*b),-cosh(beta3*b);
beta1*cos(beta1*a),-beta1*sin(beta1*a),beta1*cosh(beta1*a),beta1*sinh(beta1*a),-beta2*cos(beta2*a),beta2*sin(beta2*a),-beta2*cosh(beta2*a),-beta2*sinh(beta2*a),0,0,0,0;
0,0,0,0,beta2*cos(beta2*b),-beta2*sin(beta2*b),beta2*cosh(beta2*b),beta2*sinh(beta2*b),-beta3*cos(beta3*b),beta3*sin(beta3*b),-beta3*cosh(beta3*b),-beta3*sinh(beta3*b);
-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^2)*sin(beta1*a),-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^2)*cos(beta1*a),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^2)*sinh(beta1*a),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^2)*cosh(beta1*a),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*sin(beta2*a),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*cos(beta2*a),-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*sinh(beta2*a),-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*cosh(beta2*a),0,0,0,0;
0,0,0,0,-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*sin(beta2*b),-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*cos(beta2*b),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*sinh(beta2*b),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^2)*cosh(beta2*b),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^2)*sin(beta3*b),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^2)*cos(beta3*b),-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^2)*sinh(beta3*b),-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^2)*cosh(beta3*b);
-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^3)*cos(beta1*a),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^3)*sin(beta1*a),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^3)*cosh(beta1*a),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta1)^3)*sinh(beta1*a),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*cos(beta2*a),-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*sin(beta2*a),-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*cosh(beta2*a),-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*sinh(beta2*a),0,0,0,0;
0,0,0,0,-(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*cos(beta2*b),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*sin(beta2*b),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*cosh(beta2*b),(2*E2*z*((h2^3-h1^3)/3))*((beta2)^3)*sinh(beta2*b),((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^3)*cos(beta3*b),-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^3)*sin(beta3*b),-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^3)*cosh(beta3*b),-((2*z/3)*(E1*(h1)^3+E2*((h2)^3-(h1)^3)))*((beta3)^3)*sinh(beta3*b)]),i,optimset('display','off'));
end
ro1=root;
j_tem=1;
k_tem=1;
for i_tem=1:1999 %去重根
b_tem(i_tem)=ro1(i_tem+1);
end
b_tem(2000)=ro1(2000);
c_tem=b_tem-ro1;
c_tem_find_first= find( c_tem>0.1);
r_tem(j_tem)=ro1(c_tem_find_first(1));
r_tem(j_tem+1)=b_tem(c_tem_find_first(1));
j_tem=j_tem+2;
for i_c=c_tem_find_first(2):2000
if c_tem(i_c)>0.1
r_tem(j_tem)=b_tem(i_c);
j_tem=j_tem+1;
end
end
r_tem;
for i_re=1:j_tem-2
root_return(i_re)=r_tem(i_re+1);
end
root_return;
ri=root_return/co


Accepted Answer
More Answers (0)
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!