计算并绘制频散曲线的matlab程序

function dispersion
%绘制平板频散曲线
%tic
clc;clear;
cl=5850;    %材料纵波波速（铝板）
cs=3127;    %材料横波波速（铝板）
dfd=0.01*1e3;
fd0=（0.01:dfd/1e3:20）*1e3;  %频厚积（MHz*mm）
d_lv=2.5;

cps_min=0;
cpa_min=0;
cp_max=10000;
mode=2; %绘制的模式数
precision=1e-8;

cpa=zeros（length（fd0）mode）;    
cps=zeros（length（fd0）mode）;
for i=1:length（fd0）
    fd=fd0（i）;
    [cp12 n]=ss（cps_mincp_maxfdclcsmode）;
    for j=1:n
        cp1=cp12（j1）;
        cp2=cp12（j2）;
        cps（ij）=serfen（cp1cp2fdclcsprecision）;
    end
    
    [cp12 n]=aa（cpa_mincp_maxfdclcsmode）;
    for j=1:n
        cp1=cp12（j1）;
        cp2=cp12（j2）;
        cpa（ij）=aerfen（cp1cp2fdclcsprecision）;
    end
end

h=zeros（mode2）;
%绘制相速度
figure（1）
for j=1:2
    if j==1
        cp=cps; %对称
        color=‘b‘;
    else
        cp=cpa; %反对称
        color=‘r‘;
    end
    for i=1:mode
        cpp=cp（:i）;
        ind=find（cpp==0）;
        if ~isempty（ind）
            h（ij）=plot（（fd0（ind（end）+1:end））/d_lvcpp（ind（end）+1:end）color）;
        else
            h（ij）=plot（fd0/d_lvcppcolor）;
        end
        hold on
    end
    if j==2
        xlabel（‘f/（KHz）‘）
        ylabel（‘C_{p}/（km?s^{-1}）‘）
        title（‘5mm铝板 相速度频散曲线‘）
        set（gca‘xtick‘（0:0.6:20）*1e3/d_lv‘xticklabel‘（0:0.6:20）*1e3/d_lv）
        xlim（[0 1000]）;%
        set（gca‘ylim‘[0 cp_max]‘ytick‘（0:cp_max/1e3）*1e3...
            ‘yticklabel‘0:cp_max/1e3）
        grid on
        hSGroup = 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（‘对称模式‘‘反对称模式‘）
    end
end

%绘制群速度
figure（2）
for j=1:2
    if j==1
        cp=cps;
        color=‘b‘;
    else
        cp=cpa;
        color=‘r‘;
    end
    for i=1:mode
        cpp=cp（:i）;
        ind=find（cpp==0）;
        if ~isempty（ind）
            fd=fd0（ind（end）+1:end）‘;
            cpp=cpp（ind（end）+1:end）;
        else
            fd=fd0‘;
        end
        dcdf=diff（cpp）/dfd;
        cg=cpp（1:end-1）.^2./（cpp（1:end-1）-fd（1:end-1）.*dcdf）;
        h（ij）=plot（fd（1:end-1）/d_lvcgcolor）;
        hold on
    end
    if j==2
        xlabel（‘f/（KHz）‘）
        ylabel（‘C_{g}/（km?s^{-1}）‘）
        title（‘5mm铝板 群速度频散曲线‘）
        set（gca‘xtick‘（0:0.6:20）*1e3/d_lv‘xticklabel‘（0:0.6:20）*1e3/d_lv）
        xlim（[0 1000]）;%
        set（gca‘ylim‘[0 5.5]*1e3‘ytick‘（0:0.5:5.5）*1e3‘yticklabel‘0:0.5:5.5）
        grid on
        hSGroup = hggroup;%要在子对象构建之后构建，构建后立即使用，否则将失效
        hAGroup = hggroup;
        set（h（:1）‘parent‘hSGroup）
        set（h（:2）‘parent‘hAGroup）
        set（get（get（hSGroup‘Annotation‘）‘LegendInformation‘）...
            ‘IconDisplaystyle‘‘on‘