MATLAB求解瑞利波频散方程

% vr为未知变量，角频率w可自行确定
syms vr
f_cell = {};

for w = 1:1000     % 角频率w的范围
    N = 3; % 层数
    vp1 = [1200 1500 2500];  % 纵波波速
    vs1 = [1000 1200 1800];  % 横波波速
    ppg1 = [1800 2000 2000];  % 密度
    h1 = [6 3 inf];  % 厚度
    u = vs1.^2.*ppg1;  % 剪切模量
    I = [0 0 1 0; 0 0 0 1];

    % 计算存储各层传递矩阵
    k = w/vr;
    T_cell = cell(N-1, 1);
    for i = 1:N-1
        vp = vp1(i);
        vs = vs1(i);
        ppg = ppg1(i);
        h = h1(i);
        rp = sqrt(vr^2/vp^2 - 1);
        rs = sqrt(vr^2/vs^2 - 1);
        t = 1 - (vr^2) / (2 * vs^2);
        g = 1 - t;
        a = cos(rp * k * h);
        b = cos(rs * k * h);
        e = sin(rp * k * h) / rp;
        d = sin(rs * k * h) / rs;
        r = rp^2;
        s = rs^2;
        u1 = u(i);
        u2 = u(i+1);
        l = u2 / u1;
        T = (1 / g) * [a - b * t, e * t+d * s * (-e + d * s) * l * (b - a) * l; e * r + d * t, (b - a) * t, a - b * l, (-e * r + d * l) * l; e * r + d * t^2, (b - a) * t, (a - b * t) * l, (-e * r + d * t) * l; (a - b * t)e * t^2 + d * s, (-e + d * s) * l, (b - a) * l];
        T_cell{i} = T;
    end

    % 第N层J
    vpn = vp1(N);
    vsn = vs1(N);
    rpn = sqrt(vr^2/vpn^2 - 1);
    rsn = sqrt(vr^2/vsn^2 - 1);
    tn = 1 - (vr^2) / (2 * vsn^2);
    J = [1, 1i*rpn, 1i*rpn, tn; 1i*rsn, 1, tn, 1i*rsn];

    % 传递矩阵
    T1 = T_cell{1};
    T2 = T_cell{2};

    A = I * (T1 * T2) * J;
    f = det(A);
    f_cell{w} = f;
end

% 使用数值计算方法求解每个w对应的vr值
vr_values = zeros(1, 1000);
for w = 1:1000
    f = f_cell{w};
    vr_guess = 1000; % 初始猜测值
    vr_values(w) = newton_raphson(f, vr_guess);
end

vr_values % 输出每个w对应的vr值


% 牛顿拉弗森迭代方法
function x = newton_raphson(f, x0)
    f_d = diff(f); % 求导数
    x = x0; % 初始猜测值
    epsilon = 1e-6; % 迭代终止条件
    max_iter = 100; % 最大迭代次数
    iter = 0;

    while abs(subs(f, x)) > epsilon && iter < max_iter
        x = x - subs(f, x) / subs(f_d, x);
        iter = iter + 1;
    end
end