我的解微分方程最后一步运行不了
程序如下:Clear allformat short
symsw
L=1;b=0.08;h=0.002;h1=0.1e3;E=7.8e10;E1=7.86e10;rou=7.8e3;rou1=7.5e3;
e31=-9.29;g=9.8;
L1=0.3;L2=0.5;P1=1000;P2=1000; %几何物理参数
v=100;
F=-2*b*e31*v;
Ma=P1/g;
P=P1+P2;
I=b*h*h*h/12; %梁截面惯性矩
Ia=Ma*L*L*L/3; %质量块转动惯量
lm=((P1+P2)/(E*I)); %--λ
u=(rou*b*h*w*w/(E*I)); %--μ
EI0=E*I+E1*h*b*h1*(h/2+h1); %--(EI)*
lm1=(((P1+P2)-F)/EI0); %--λ1
u1=(w*w*(rou*b*h+2*rou1*b*h1)/EI0); %--μ1
k1=((u+lm*lm/4).^(1/2)+lm/2).^(1/2);
k2=((u+lm*lm/4).^(1/2)-lm/2).^(1/2);
k11=((u1+lm1*lm1/4).^(1/2)+lm1/2).^(1/2);
k22=((u1+lm1*lm1/4).^(1/2)-lm1/2).^(1/2);
r1r1=1;r1r3=1;r2r2=k1;r2r4=k2; %P>=F时矩阵元素
r3r9=k1*k1*k1*sin(k1*L)+P1*k1*sin(k1*L)/E*I;
r3r10=-k1*k1*k1*cos(k1*L) - P1*k1*cos(k1*L)/E*I;
r3r11=k2*k2*k2*sinh(k2*L) - P1*k1*sinh(k1*L)/E*I;
r3r12=k2*k2*k2*cosh(k2*L) - P1*k1*cosh(k1*L)/E*I;
r4r9=-k1*k1*cos(k1*L)-w*w*Ia*k1*sin(k1*L)/(E*I);
r4r10=-k1*k1 *sin(k1*L)+ w*w*Ia*k1*cos(k1*L)/(E*I);
r4r11=k2*k2* cosh(k1*L)+w*w*Ia* k2*sinh(k2*L)/(E*I);
r4r12=k2*k2* sinh(k2*L)+w*w*Ia* k2*cosh(k2*L)/(E*I);
r5r1=cos(k1*L1);r5r2=sin(k1*L1);
r5r3=cosh(k2*L1);r5r4=sinh(k2*L1);
r5r5=-cos(k11*L1);r5r6=-sin(k11*L1);
r5r7=-cosh(k22*L1);r5r8=-sinh(k22*L1);
r6r1=-k1*sin(k1*L1);r6r2=k1*cos(k1*L1);
r6r3=k2*sinh(k2*L1);r6r4=k2*cosh(k2*L1);
r6r5=k11*sin(k11*L1);r6r6=-k11*cos(k11*L1);
r6r7=-k22*sinh(k22*L1);r6r8=-k22*cosh(k22*L1);
r7r1=-k1*k1*cos(k1*L1);r7r2=-k1*k1*sin(k1*L1);
r7r3=k2*k2*cosh(k2*L1);r7r4=k2*k2*sinh(k2*L1);
r7r5=k11*k11*cos(k11*L1);r7r6=k11*k11*sin(k11*L1);
r7r7=-k22*k22*cosh(k22*L1);r7r8=-k22*k22*sinh(k22*L1);
r8r1=k1*k1*k1*sin(k1*L1);r8r2=-k1*k1*k1*cos(k1*L1);
r8r3=k2*k2*k2*sinh(k2*L1);r8r4=k2*k2*k2*cosh(k2*L1);
r8r5=-k11*k11*k11*sin(k11*L1);r8r6=k11*k11*k11*cos(k11*L1);
r8r7=-k22*k22*k22*sinh(k22*L1);r8r8=-k22*k22*k22*cosh(k22*L1);
r9r5=cos(k11*L2);r9r6=sin(k11*L2);
r9r7=cosh(k22*L2);r9r8=sinh(k22*L2);
r9r9=-cos(k1*L2);r9r10=-sin(k1*L2);
r9r11=-cosh(k2*L2);r9r12=-sinh(k2*L2);
r10r5=-k11*sin(k11*L2);r10r6=k11*cos(k11*L2);
r10r7=k22*sinh(k22*L2);r10r8=k22*cosh(k22*L2);
r10r9=k1*sin(k1*L2);r10r10=-k1*cos(k1*L2);
r10r11=-k2*sinh(k2*L2);r10r12=-k2*cosh(k2*L2);
r11r5=-k11*k11*cos(k11*L2);r11r6=-k11*k11*sin(k11*L2);
r11r7=k22*k22*cosh(k22*L2);r11r8=k22*k22*sinh(k22*L2);
r11r9=k1*k1*cos(k1*L2);r11r10=k1*k1*sin(k1*L2);
r11r11=-k2*k2*cosh(k2*L2);r11r12=-k2*k2*sinh(k2*L2);
r12r5=k11*k11*k11*sin(k11*L2);r12r6=-k11*k11*k11*cos(k11*L2);
r12r7=k22*k22*k22*sinh(k22*L2);r12r8=k22*k22*k22*cosh(k22*L2);
r12r9=-k1*k1*k1*sin(k1*L2);r12r10=k1*k1*k1*cos(k1*L2);
r12r11=-k2*k2*k2*sinh(k2*L2);r12r12=-k2*k2*k2*cosh(k2*L2);
r5r5r=-cosh(k11*L1);r5r6r=-sinh(k11*L1); %P<F时部分矩阵元素
r5r7r=-cos(k22*L1);r5r8r=-sin(k22*L1);
r6r5r=-k11*sinh(k11*L1);r6r6r=-k11*cosh(k11*L1);
r6r7r=k22*sin(k22*L1);r6r8r=-k22*cos(k22*L1);
r7r5r=-k11*k11*cosh(k11*L1);r7r6r=-k11*k11*sinh(k11*L1);
r7r7r=k22*k22*cos(k22*L1);r7r8r=k22*k22*sin(k22*L1);
r8r5r=-k11*k11*k11*sinh(k11*L1);r8r6r=-k11*k11*k11*cosh(k11*L1);
r8r7r=-k22*k22*k22*sin(k22*L1);r8r8r=k22*k22*k22*cos(k22*L1);
r9r5r=cosh(k11*L2);r9r6r=sinh(k11*L2);
r9r7r=cos(k22*L2);r9r8r=sin(k22*L2);
r10r5r=-k11*sinh(k11*L2);r10r6r=k11*cosh(k11*L2);
r10r7r=-k22*sin(k22*L2);r10r8r=k22*cos(k22*L2);
r11r5r=k11*k11*cosh(k11*L2);r11r6r=k11*k11*sinh(k11*L2);
r11r7r=-k22*k22*cos(k22*L2);r11r8r=-k22*k22*sin(k22*L2);
r12r5r=k11*k11*k11*sinh(k11*L2);r12r6r=k11*k11*k11*cosh(k11*L2);
r12r7r=k22*k22*k22*sin(k22*L2);r12r8r=-k22*k22*k22*cos(k22*L2);
R1=[r1r1,0,r1r3,0,0,0,0,0,0,0,0,0; %P>=F矩阵形式
0,r2r2,0,r2r4,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,r3r9,r3r10,r3r11,r3r12;
0,0,0,0,0,0,0,0,r4r9,r4r10,r4r11,r4r12;
r5r1,r5r2,r5r3,r5r4,r5r5,r5r6,r5r7,r5r8,0,0,0,0;
r6r1,r6r2,r6r3,r6r4,r6r5,r6r6,r6r7,r6r8,0,0,0,0;
r7r1,r7r2,r7r3,r7r4,r7r5,r7r6,r7r7,r7r8,0,0,0,0;
r8r1,r8r2,r8r3,r8r4,r8r5,r8r6,r8r7,r8r8,0,0,0,0;
0,0,0,0,r9r5,r9r6,r9r7,r9r8,r9r9,r9r10,r9r11,r9r12;
0,0,0,0,r10r5,r10r6,r10r7,r10r8,r10r9,r10r10,r10r11,r10r12;
0,0,0,0,r11r5,r11r6,r11r7,r11r8,r11r9,r11r10,r11r11,r11r12;
0,0,0,0,r12r5,r12r6,r12r7,r12r8,r12r9,r12r10,r12r11,r12r12;];
R2=[r1r1,0,r1r3,0,0,0,0,0,0,0,0,0; %P<F矩阵形式
0,r2r2,0,r2r4,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,r3r9,r3r10,r3r11,r3r12;
0,0,0,0,0,0,0,0,r4r9,r4r10,r4r11,r4r12;
r5r1,r5r2,r5r3,r5r4,r5r5r,r5r6r,r5r7r,r5r8r,0,0,0,0;
r6r1,r6r2,r6r3,r6r4,r6r5r,r6r6r,r6r7r,r6r8r,0,0,0,0;
r7r1,r7r2,r7r3,r7r4,r7r5r,r7r6r,r7r7r,r7r8r,0,0,0,0;
r8r1,r8r2,r8r3,r8r4,r8r5r,r8r6r,r8r7r,r8r8r,0,0,0,0;
0,0,0,0,r9r5r,r9r6r,r9r7r,r9r8r,r9r9,r9r10,r9r11,r9r12;
0,0,0,0,r10r5r,r10r6r,r10r7r,r10r8r,r10r9,r10r10,r10r11,r10r12;
0,0,0,0,r11r5r,r11r6r,r11r7r,r11r8r,r11r9,r11r10,r11r11,r11r12;
0,0,0,0,r12r5r,r12r6r,r12r7r,r12r8r,r12r9,r12r10,r12r11,r12r12;];
ifP>=F
R=R1
else
R=R2
end
M=det(R) %从此开始显示目标太大???
solve(M)
页:
[1]