Matlab绘制频散曲线程序代码

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_min,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);end[cp12n]=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';elsecp=cpa;color='r';endfori=1:modecpp=cp(:,i);ind=find(cpp==0);if~isempty(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}/(km·s^{-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',[0cp_max],'ytick',(0:cp_max/1e3)*1e3,...'yticklabel',0:cp_max/1e3)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'),'LegendInformation'),...'IconDisplayStyle','on');legend('对称模式','反对称模式')endend%群速度figure(2)forj=1:2ifj==1cp=cps;color='b';elsecp=cpa;color='r';endfori=1:modecpp=cp(:,i);ind=find(cpp==0);if~isempty(ind)fd=fd0(ind(end)+1:end)';cpp=cpp(ind(end)+1:end);elsefd=fd0';enddcdf=diff(cpp)/dfd;cg=cpp(1:end-1).^2./(cpp(1:end-1)-fd(1:end-1).*dcdf);h(i,j)=plot(fd(1:end-1)/d_Q235,cg,color);holdonendifj==2xlabel('f/(KHz)')ylabel('C_{g}/(km·s^{-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)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'),'LegendInformation'),...'IconDisplayStyle','on');legend('对称模式','反对称模式')endend%tocend%对称模式function[cp0n]=ss(cp_min,cp_max,fd,cl,cs,mode)cp2=cp_min;deter=33;cp0=zeros(mode,2);n=0;whilecp2cp_max&&nmodecp1=cp2;cp2=cp1+deter;y1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs);whiley1*y20&&cp2cp_maxcp1=cp2;cp2=cp1+deter;y1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs);endify1*y20n=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(n,:)=[cp2cp2];cp2=cp2+1;endendendendfunctionfs=smode(cp,fd,cl,cs)p=abs(sqrt((cp/cs)^2-1));q=abs(sqrt((cp/cl)^2-1));ifcp=cs%p和q都是复数fs=-4*p*q*sinh(pi*fd/cp*q)*cosh(pi*fd/cp*p)+(-p^2-1)^2*sinh(pi*fd/cp*p)*cosh(pi*fd/cp*q);elseifcpcs&&cp=cl%p是实数，q是复数fs=-4*p*q*sinh(pi*fd/cp*q)*cos(pi*fd/cp*p)+(p^2-1)^2*sin(pi*fd/cp*p)*cosh(pi*fd/cp*q);elsefs=4*p*q*sin(pi*fd/cp*q)*cos(pi*fd/cp*p)+(p^2-1)^2*sin(pi*fd/cp*p)*cos(pi*fd/cp*q);endendfunctioncp=serfen(cp1,cp2,fd,cl,cs,precision)whilecp2-cp1precisiony1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs);cp0=(cp1+cp2)/2;y0=smode(cp0,fd,cl,cs);ify1*y00cp2=cp0;elseify2*y00cp1=cp0;elseify0==0breakelseify1==0cp2=cp1;breakelseify2==0cp1=cp2;breakendendcp=(cp2+cp1)/2;end%反对称模式function[cp0n]=aa(cp_min,cp_max,fd,cl,cs,mode)cp2=cp_min;deter=33;cp0=zeros(mode,2);n=0;whilecp2cp_max&&nmodecp1=cp2;cp2=cp1+deter;y1=amode(cp1,fd,cl,cs);y2=amode(cp2,fd,cl,cs);whiley1*y20&&cp2cp_maxcp1=cp2;cp2=cp1+deter;y1=amode(cp1,fd,cl,cs);y2=amode(cp2,fd,cl,cs);endify1*y20n=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(n,:)=[cp2cp2];cp2=cp2+1;endendendendfunctionfs=amode(cp,fd,cl,cs)p=abs(sqrt((cp/cs)^2-1));q=abs(sqrt((cp/cl)^2-1));ifcp=cs%p和q都是复数fs=-4*p*q*sinh(pi*fd/cp*p)*cosh(pi*fd/cp*q)+(-p^2-1)^2*sinh(pi*fd/cp*q)*cosh(pi*fd/cp*p);elseifcpcs&&cp=cl%p是实数，q是复数fs=4*p*q*sin(pi*fd/cp*p)*cosh(pi*fd/cp*q)+(p^2-1)^2*sinh(pi*fd/cp*q)*cos(pi*fd/cp*p);elsefs=4*p*q*sin(pi*fd/cp*p)*cos(pi*fd/cp*q)+(p^2-1)^2*sin(pi*fd/cp*q)*cos(pi*fd/cp*p);endendfunctioncp=aerfen(cp1,cp2,fd,cl,cs,precision)whilecp2-cp1precisiony1=amode(cp1,fd,cl,cs);y2=amode(cp2,fd,cl,cs);cp0=(cp1+cp2)/2;y0=amode(cp0,fd,cl,cs);ify1*y00cp2=cp0;elseify2*y00cp1=cp0;elseify0==0breakelseify1==0cp2=cp1;breakelseify2==0cp1=cp2;breakendendcp=(cp2+cp1)/2;end