Matlab绘制频散曲线程序代码

1、function disper %绘制平板频散曲线 %tic clc;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);for i=1:length(fd0) fd=fd0(i)

2、;cp12 n=ss(cps_min,cp_max,fd,cl,cs,mode); for j=1:ncp1=cp12(j,1);cp2=cp12(j,2); cps(i,j)=serfen(cp1,cp2,fd,cl,cs,precision);endcp12 n=aa(cpa_min,cp_max,fd,cl,cs,mode); for j=1:ncp1=cp12(j,1);cp2=cp12(j,2); cpa(i,j)=aerfen(cp1,cp2,fd,cl,cs,precision);end endh=zeros(mode,2);%相速度figure(1)for j=1:2 if j

3、=1cp=cps;color=b;elsecp=cpa;color=r;endfor i=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);else h(i,j)=plot(fd0/d_Q235,cpp,color);endhold onendif j=2xlabel(f/(KHz)ylabel(C_p/(km sA-1)title(6mm 钢板 相速度频散曲线 ) set(gca,xtick,(0:0.6:20

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

5、rmation),.IconDisplayStyle,on); set(get(get(hAGroup,Annotation),LegendInformation),.IconDisplayStyle,on);legend(对称模式,反对称模式)endend%群速度figure(2)for j=1:2if j=1cp=cps;color=b;elsecp=cpa;color=r;endfor i=1:modecpp=cp(:,i);ind=find(cpp=0);if isempty(ind)fd=fd0(ind(end)+1:end); cpp=cpp(ind(end)+1:end);els

6、efd=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,color);hold onendif j=2xlabel(f/(KHz)ylabel(C_g/(km sA-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,0 5.

7、5*1e3,ytick,(0:0.5:5.5)*1e3,yticklabel,0:0.5:5.5) grid onhSGroup = 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),.I

8、conDisplayStyle,on);legend(对称模式,反对称模式)endend%tocend%对称模式function cp0 n=ss(cp_min,cp_max,fd,cl,cs,mode) cp2=cp_min;deter=33;cp0=zeros(mode,2);n=0;while cp2cp_max&n0&cp2cp_maxcp1=cp2;cp2=cp1+deter;y1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs);endif y1*y20n=n+1;cp0(n,:)=cp1 cp2;elseif y1=0&y2=0n=n+1;cp

9、0(n,:)=cp1 cp1;elseif y2=0&y1=0 n=n+1;cp0(n,:)=cp2 cp2; cp2=cp2+1;elseif y1=0&y2=0 n=n+1; cp0(n,:)=cp1 cp1; n=n+1;cp0(n,:)=cp2 cp2; cp2=cp2+1;endendendendfunction fs=smode(cp,fd,cl,cs)p=abs(sqrt(cp/cs)A2-1);q=abs(sqrt(cp/clF2-1);if cpcs&cpprecisiony1=smode(cp1,fd,cl,cs);y2=smode(cp2,fd,cl,cs); cp0=(

10、cp1+cp2)/2;y0=smode(cp0,fd,cl,cs);if y1*y00cp2=cp0;elseif y2*y00cp1=cp0;elseif y0=0breakelseif y1=0cp2=cp1;breakelseif y2=0cp1=cp2;breakendendcp=(cp2+cp1)/2;end %反对称模式function cp0 n=aa(cp_min,cp_max,fd,cl,cs,mode) cp2=cp_min;deter=33;cp0=zeros(mode,2);n=0;while cp2cp_max&n0&cp2cp_maxcp1=cp2;cp2=cp1+

11、deter;y1=amode(cp1,fd,cl,cs);y2=amode(cp2,fd,cl,cs);endif y1*y20n=n+1;cp0(n,:)=cp1 cp2;elseif y1=0&y2=0n=n+1;cp0(n,:)=cp1 cp1;elseif y2=0&y1=0n=n+1;cp0(n,:)=cp2 cp2;cp2=cp2+1;elseif y1=0&y2=0n=n+1;cp0(n,:)=cp1 cp1;n=n+1;cp0(n,:)=cp2 cp2;cp2=cp2+1;endendendendfunction fs=amode(cp,fd,cl,cs)p=abs(sqrt(cp/cs)A2-1);q=abs(sqrt(cp/clF2-1);if cpcs&cpprecision y1=amode(cp1,fd,cl,cs); y2=amode(cp