Matlab绘制频散曲线程序代码

1、functiondisper%绘制平板频散曲线%ticclc;clear;cl=5790;%材料纵波波速(钢板)cs=3200;%材料横波波速(钢板)dfd=0.01*1e3;fd0=(0.01:dfd/1e3:20)*1e3;%频厚积(MHz*mm)d_Q235=6;cps_min=2700;cpa_min=100;cp_max=10000;mode=3;%绘制的模式数precision=1e-8;cpa=zeros(length(fd0),mode);cps=zeros(length(fd0),mode);fori=1:length(fd0)fd=fd0(i);cp12n=ss(cps_m

2、in,cp_max,fd,cl,cs,mode);forj=1:ncp1=cp12(j,1);cp2=cp12(j,2);cps(i,j)=serfen(cp1,cp2,fd,cl,cs,precision);endcp12n=aa(cpa_min,cp_max,fd,cl,cs,mode);forj=1:ncp1=cp12(j,1);cp2=cp12(j,2);cpa(i,j)=aerfen(cp1,cp2,fd,cl,cs,precision);endendh=zeros(mode,2);%相速度figure(1)forj=1:2ifj=1cp=cps;color='b'e

3、lsecp=cpa;color='r'endfori=1:modecpp=cp(:,i);ind=find(cpp=0);ifisempty(ind)h(i,j)=plot(fd0(ind(end)+1:end)/d_Q235,cpp(ind(end)+1:end),color);elseh(i,j)=plot(fd0/d_Q235,cpp,color);endholdonendifj=2xlabel('f/(KHz)')ylabel('C_p/(kmsA-1)')title('6mm钢板相速度频散曲线')set(gca,'

4、;xtick',(0:0.6:20)*1e3/d_Q235,'xticklabel',(0:0.6:20)*1e3/d_Q235)xlim(0,1000);%set(gca,'ylim',0cp_max,'ytick',(0:cp_max/1e3)*1e3,.'yticklabel',0:cp_max/1e3)gridonhSGroup=hggroup;%要在子对象构建之后构建，构建后立即使用，否则将失效hAGroup=hggroup;set(h(:,1),'parent',hSGroup)set(h(:,

5、2),'parent',hAGroup)set(get(get(hSGroup,'Annotation'),'LegendInformation'),.'IconDisplayStyle','on');set(get(get(hAGroup,'Annotation'),'LegendInformation'),.'IconDisplayStyle','on');legend('对称卞式',反对称模式)endend%群速度figure(2

6、)forj=1:2ifj=1cp=cps;color='b'elsecp=cpa;color='r'endfori=1:modecpp=cp(:,i);ind=find(cpp=0);ifisempty(ind)fd=fd0(ind(end)+1:end)'cpp=cpp(ind(end)+1:end);elsefd=fd0'enddcdf=diff(cpp)/dfd;cg=cpp(1:end-1).A2./(cpp(1:end-1)-fd(1:end-1).*dcdf);h(i,j)=plot(fd(1:end-1)/d_Q235,cg,col

7、or);holdonendifj=2xlabel('f/(KHz)')ylabel('C_g/(kmsA-1)')title('6mm钢板群速度频散曲线')set(gca,'xtick',(0:0.6:20)*1e3/d_Q235,'xticklabel',(0:0.6:20)*1e3/d_Q235)xlim(0,1000);%set(gca,'ylim',05.5*1e3,'ytick',(0:0.5:5.5)*1e3,'yticklabel',0:0.5:5.5)

8、gridonhSGroup=hggroup;%要在子对象构建之后构建，构建后立即使用，否则将失效hAGroup=hggroup;set(h(:,1),'parent',hSGroup)set(h(:,2),'parent',hAGroup)set(get(get(hSGroup,'Annotation'),'LegendInformation'),.'IconDisplayStyle','on');set(get(get(hAGroup,'Annotation'),'Lege

9、ndInformation'),.'IconDisplayStyle','on');legend('对称卞式',反对称模式)endend%tocend%对称模式functioncp0n=ss(cp_min,cp_max,fd,cl,cs,mode)cp2=cp_min;deter=33;cp0=zeros(mode,2);n=0;whilecp2<cp_max&&n<modecp1=cp2;cp2=cp1+deter;y1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs);w

10、hiley1*y2>0&&cp2<cp_maxcp1=cp2;cp2=cp1+deter;y1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs);endify1*y2<0n=n+1;cp0(n,:)=cp1cp2;elseify1=0&&y2=0n=n+1;cp0(n,:)=cp1cp1;elseify2=0&&y1=0n=n+1;cp0(n,:)=cp2cp2;cp2=cp2+1;elseify1=0&&y2=0n=n+1;cp0(n,:)=cp1cp1;n=n+1;cp0(

11、n,:)=cp2cp2;cp2=cp2+1;endendendendfunctionfs=smode(cp,fd,cl,cs)p=abs(sqrt(cp/cs)A2-1);q=abs(sqrt(cp/cl)A2-1);ifcp<=cs%p和q都是复数fs=-4*p*q*sinh(pi*fd/cp*q)*cosh(pi*fd/cp*p)+(-pA2-1)A2*sinh(pi*fd/cp*p)*cosh(pi*fd/cp*q);elseifcp>cs&&cp<=cl%p是实数，q是复数fs=-4*p*q*sinh(pi*fd/cp*q)*cos(pi*fd/cp*

12、p)+(pA2-1)A2*sin(pi*fd/cp*p)*cosh(pi*fd/cp*q);elsefs=4*p*q*sin(pi*fd/cp*q)*cos(pi*fd/cp*p)+(pA2-1)A2*sin(pi*fd/cp*p)*cos(pi*fd/cp*q);endendfunctioncp=serfen(cp1,cp2,fd,cl,cs,precision)whilecp2-cp1>precisiony1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs);cp0=(cp1+cp2)/2;y0=smode(cp0,fd,cl,cs);ify1*y

13、0<0cp2=cp0;elseify2*y0<0cp1=cp0;elseify0=0breakelseify1=0cp2=cp1;breakelseify2=0cp1=cp2;breakendendcp=(cp2+cp1)/2;end%反对称模式functioncp0n=aa(cp_min,cp_max,fd,cl,cs,mode)cp2=cp_min;deter=33;cp0=zeros(mode,2);n=0;whilecp2<cp_max&&n<modecp1=cp2;cp2=cp1+deter;y1=amode(cp1,fd,cl,cs);y2=

14、amode(cp2,fd,cl,cs);whiley1*y2>0&&cp2<cp_maxcp1=cp2;cp2=cp1+deter;y1=amode(cp1,fd,cl,cs);y2=amode(cp2,fd,cl,cs);endify1*y2<0n=n+1;cp0(n,:)=cp1cp2;elseify1=0&&y2=0n=n+1;cp0(n,:)=cp1cp1;elseify2=0&&y1=0n=n+1;cp0(n,:)=cp2cp2;cp2=cp2+1;elseify1=0&&y2=0n=n+1;cp0(n

15、,:)=cp1cp1;n=n+1;cp0(n,:)=cp2cp2;cp2=cp2+1;endendendendfunctionfs=amode(cp,fd,cl,cs)p=abs(sqrt(cp/cs)A2-1);q=abs(sqrt(cp/cl)A2-1);ifcp<=cs%p和q都是复数fs=-4*p*q*sinh(pi*fd/cp*p)*cosh(pi*fd/cp*q)+(-pA2-1)A2*sinh(pi*fd/cp*q)*cosh(pi*fd/cp*p);elseifcp>cs&&cp<=cl%p是实数，q是复数fs=4*p*q*sin(pi*fd/cp*p)*cosh(pi*fd/cp*q)+(pA2-1)A2*sinh(pi*fd/cp*q)*cos(pi*fd/cp*p);fs=4*p*q*sin(pi*fd/cp*p)*cos(pi*fd/cp*q)+(pA2-142*sin(pi*fd/cp*q)*cos(pi*fd/cp*p);endendfunctioncp=aerfen(cp1,cp2,fd,cl,cs,pre