设置的符号变量如何参与下面的数值计算，并求出这个变量的具体值，下面的sym f 就是我要求出来的量
clear all;
clc;clear;
% plate properties
Lx=4.5;%lenth
Ly=4.5;%width
h1=0.2;%high
E10=3e10;%youngs modulus
rho1=1750;%density
mu=0.3;%poissons ratio
D=E10*h1^3./(12*(1-mu^2));%flexaral rigidity
k1=5.5e7;%reaction coefficient
kp=k1./D;
sym f;
omega=2*pi*f;
a=omega.^2;
Lmda=kp-rho1*h1*a/D;
num1=5;
num2=5;
A1=zeros(num1+1);
A2=zeros(num1+1);
B1=zeros(num2+1);
B2=zeros(num2+1);
C=zeros(num1+1);
D=zeros(num1+1);
E=zeros(num2+1);
E1=zeros(num2+1);
F=zeros(num1+1);
P=zeros(num1+1);
Q=zeros(num1+1);
for n=0:num2
kn=n*pi./Lx;
for m=1:2:num1
km=m*pi./Ly;
C(m,n+1)=(km^2+mu*kn^2)./(Lx*Ly*((km^4+2*(km^2)*(kn^2)+kn^4)^2+Lmda));
D(m,n+1)=(kn^2+mu*km^2)./(Lx*Ly*((km^4+2*(km^2)*(kn^2)+kn^4)^2+Lmda));
A2(m,n+1)=(km^2+mu*kn^2)*D(m,n+1);
E(m,n+1)=(km^2+mu*kn^2)*C(m,n+1);
end
end
E1=E';
s1=sum(A2);
for m=1 num1+1);
A1(m,m)=s1(1,m);
end
for m=0:num1
km=m*pi./Lx;
for n=1:2:num2
kn=n*pi./Ly;
P(m+1,n)=(km^2+mu*kn^2)./(Lx*Ly*((km^4+2*(km^2)*(kn^2)+kn^4)^2+Lmda));
Q(m+1,n)=(kn^2+mu*km^2)./(Lx*Ly*((km^4+2*(km^2)*(kn^2)+kn^4)^2+Lmda));
B2(m+1,n)=(km^2+mu*kn^2)*P(m+1,n);
F(m+1,n)=(km^2+mu*kn^2)*Q(m+1,n);
end
end
s2=sum(B2,2);
for n=1num2+1);
B1(n,n)=s2(n,1);
end
mat=sym('[A1,E1;F,B1]');
solve(det(mat),'f')
|
|
