马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我想画图,但是子函数中t是syms,画图用plot要求double,怎么转化啊?我把程序贴出来大家帮我弄弄,小妹多谢了!!!!!!!!
主函数
nr=9;
ns=5;
Lx=1.0;
Ly=1.5;
L=0.75;
c=340;
f=200;
k=2*pi*f/c;
omga=2*pi*f;
Ds=2790;
t=0.002;
v=0.3;
xp0=0.5;
yp0=1;
xp1=1.25;
yp1=1;
xb=0.5;
yb=1;
Da=1.225;
E=9.0e10;
Fload=1;
D=E*t^3/(12*(1-v^2));
sp=sqrt(Ds*t*omga^2/D);
kb=sqrt(sp);
na = 2*(nr+1) + 2*(ns+1);
t=0:0.1:0.75;
z=dis(x,y,nr,ns,Lx,Ly,kb,k,xp0,yp0,xp1,yp1,L,D,xb,Fload);
plot(t,real(z))
子函数
function z=dis(x,y,nr,ns,Lx,Ly,kb,k,xp0,yp0,xp1,yp1,L,D,xb,Fload)
syms t
pusi1=exp(-kb*t*j);
pusi2=exp(kb*(t-L));
pusi3=exp(-kb*t*j);
pusi4=exp(-kb*t);
stru=[pusi1 pusi2 pusi3 pusi4];
ws=[ -.68423e-7+.34252e-6*i
-.74339e-7-.71546e-7*i
.65116e-7-.34011e-6*i
-.39375e-7+.87761e-8*i];
A1=stru*ws;
pa=[7.0052-1.1545*i
-4.7501+2.1268*i
-1.5536+1.3057*i
2.0828+.47272*i
17.979+25.869*i
22.266-3.8239*i
9.6485-5.9311*i
-53.923+35.741*i
-249.20-1.7802*i
195.63-141.19*i
435.52+403.67*i
-259.41+261.85*i
1234.6-168.58*i
-1224.5+135.76*i
5980.3-20377.*i
8742.0-2439.5*i
22017.-12866.*i
-24877.+7328.0*i
.22660e7+.47022e7*i
37481.-10130.*i
1.7396-.47984*i
-2.6366+.12197*i
1.9372+.29824*i
.45464-.11423*i
2.4070-3.3590*i
1.2118-.42206*i
7.0606+2.7055*i
.24693-.89894e-1*i
6.5731+12.141*i
.11083+.68151e-1*i
-.35012e-1-3.1229*i
.33063e-1+.34949e-2*i];
if xp1 == xp0
cosp = 1;
sinp = 0;
L = yp1 - yp0;
else
L = sqrt( (xp1-xp0)^2 + (yp1-yp0)^2 );
cosp = (xp1-xp0)/L;
sinp = (yp1-yp0)/L;
end
% 建立波函数矩阵
for ir = 1:2*(nr+1);
r = fix(ir/2) - 1 + rem(ir,2);
%--------- r segment begin --------------------
% Kr-x,Kr-y for r
if (rem(ir,2) == 1)
krx=r*pi/Lx;
kry=sqrt(k^2 - krx^2);
fir=cos(krx*(xp0+t*cosp))*exp(-j*kry*(yp0+t*sinp-Ly));
else
krx=r*pi/Lx;
kry=-sqrt(k^2 - krx^2);
fir = cos(krx*(xp0+t*cosp)) * exp(-j*kry*(yp0+t*sinp));
end;
war(ir)= fir/(D*(((krx*cosp+kry*sinp)^4-kb^4)));
end
%
for js = 1:2*(ns+1);
s = fix(js/2) - 1 + rem(js,2);
%--------- s segment begin --------------------
% Ks-x,Ks-y for s
if (rem(js,2) == 1)
ksy= s*pi/Ly;
ksx= sqrt(k^2 - ksy^2);
fjs = cos(ksy*(yp0+t*sinp))*exp(-j*ksx*(xp0+t*cosp-Lx));
else
ksy= s*pi/Ly;
ksx=-sqrt(k^2 - ksy^2);
fjs = cos(ksy*(yp0+t*sinp))*exp(-j*ksx*(xp0+t*cosp));
end;
was(js)= fjs/(D*(((ksx*cosp+ksy*sinp)^4-kb^4)));
end
wa=[war was];
A2=wa*pa;
fdk3=Fload/(4.0*D*kb^3);
fdk2=Fload/(4.0*D*kb^2);
absf=sqrt((t-xb)*(t-xb));
%
wf=j*fdk3*(exp(-j*kb*absf)-j*exp(-kb*absf));
z=A1+A2+wf;
return |