%% 分析监管收益D对演化的影响
function analyze_regulation_benefit_parameter()
    
    fprintf('\n开始监管收益D的敏感性分析...\n');
    
    % 第一组参数设置
    params1 = struct();
    params1.k1 = 1;
    params1.C2 = 50;
    params1.C3 = 20;
    params1.p1 = 200;
    params1.C = 40;
    params1.m = 50;
    params1.a = 0.5;
    params1.n = 80;
    params1.k2 = 1;
    
    % 第二组参数设置
    params2 = struct();
    params2.k1 = 1;
    params2.C2 = 30;
    params2.C3 = 20;
    params2.p1 = 200;
    params2.C = 40;
    params2.m = 50;
    params2.a = 0.5;
    params2.n = 80;
    params2.k2 = 1;
    
    % D值设置
    D_values = [50, 100, 150];
    
    % 创建图形窗口
    figure('Position', [100, 100, 1200, 600]);
    
    % 第一组分析
    subplot(1, 2, 1);
    analyze_D_group(params1, D_values, '第1组: 满足条件3');
    
    % 第二组分析
    subplot(1, 2, 2);
    analyze_D_group(params2, D_values, '第2组: 满足条件2和3');
    
    sgtitle('监管收益D的敏感性分析（图10）', 'FontSize', 14, 'FontWeight', 'bold');
    
    % 生成时间序列分析图
    generate_D_time_series(params1, params2, D_values);
end

function analyze_D_group(base_params, D_values, title_text)
    hold on;
    
    % 颜色和样式设置
    colors = {'r', 'g', 'b'};
    markers = {'+', 'x', '*'};
    line_styles = {'-', '--', ':'};
    
    % 时间跨度
    tspan = [0, 50];
    
    for D_idx = 1:length(D_values)
        base_params.D = D_values(D_idx);
        
        % 运行多次仿真
        num_sims = 10;
        
        for sim = 1:num_sims
            % 随机初始条件
            initial_state = [rand(); rand(); rand()];
            
            % 求解演化方程
            [t, trajectory] = ode45(@(t, state) evolutionary_dynamics_D(t, state, base_params), ...
                                    tspan, initial_state);
            
            % 绘制轨迹
            if sim == 1
                h = plot3(trajectory(:, 1), trajectory(:, 2), trajectory(:, 3), ...
                         'Color', colors{D_idx}, 'LineWidth', 2, ...
                         'LineStyle', line_styles{D_idx});
            else
                plot3(trajectory(:, 1), trajectory(:, 2), trajectory(:, 3), ...
                     'Color', [hex2rgb(colors{D_idx}), 0.3], 'LineWidth', 1, ...
                     'LineStyle', line_styles{D_idx});
            end
            
            % 标记终点
            if sim == 1
                scatter3(trajectory(end, 1), trajectory(end, 2), trajectory(end, 3), ...
                        100, colors{D_idx}, markers{D_idx}, 'LineWidth', 2);
            end
        end
    end
    
    % 设置图形属性
    xlabel('x (监管部门)');
    ylabel('y (共享方A)');
    zlabel('z (共享方B)');
    title(title_text);
    grid on;
    view(45, 30);
    xlim([0, 1]);
    ylim([0, 1]);
    zlim([0, 1]);
    
    % 添加图例
    legend(arrayfun(@(i) sprintf('D=%d', D_values(i)), 1:length(D_values), 'UniformOutput', false), ...
           'Location', 'best');
    
    hold off;
end

function generate_D_time_series(params1, params2, D_values)
    figure('Position', [100, 100, 1400, 600]);
    
    % 时间跨度
    tspan = [0, 50];
    initial_state = [0.5; 0.5; 0.5];
    
    % 第一组时间序列
    subplot(2, 2, 1);
    plot_D_time_series(params1, D_values, tspan, initial_state, 'x');
    title('第1组: 监管部门策略演化');
    
    subplot(2, 2, 2);
    plot_D_time_series(params1, D_values, tspan, initial_state, 'y');
    title('第1组: 共享方A策略演化');
    
    % 第二组时间序列
    subplot(2, 2, 3);
    plot_D_time_series(params2, D_values, tspan, initial_state, 'x');
    title('第2组: 监管部门策略演化');
    
    subplot(2, 2, 4);
    plot_D_time_series(params2, D_values, tspan, initial_state, 'z');
    title('第2组: 共享方B策略演化');
    
    sgtitle('监管收益D对策略演化时间序列的影响', 'FontSize', 14, 'FontWeight', 'bold');
end

function plot_D_time_series(params, D_values, tspan, initial_state, variable)
    hold on;
    
    colors = {'r', 'g', 'b'};
    line_widths = [2, 2, 2];
    
    for D_idx = 1:length(D_values)
        params.D = D_values(D_idx);
        
        % 求解演化方程
        [t, trajectory] = ode45(@(t, state) evolutionary_dynamics_D(t, state, params), ...
                                tspan, initial_state);
        
        % 选择要绘制的变量
        if variable == 'x'
            data = trajectory(:, 1);
        elseif variable == 'y'
            data = trajectory(:, 2);
        else
            data = trajectory(:, 3);
        end
        
        % 绘制时间序列
        plot(t, data, 'Color', colors{D_idx}, 'LineWidth', line_widths(D_idx));
    end
    
    xlabel('时间 t');
    ylabel(sprintf('策略概率 %s', variable));
    legend(arrayfun(@(i) sprintf('D=%d', D_values(i)), 1:length(D_values), 'UniformOutput', false), ...
           'Location', 'best');
    grid on;
    xlim([0, 50]);
    ylim([0, 1]);
    
    hold off;
end

function dstate_dt = evolutionary_dynamics_D(~, state, params)
    x = state(1);
    y = state(2);
    z = state(3);
    
    % 提取参数
    k1 = params.k1;
    k2 = params.k2;
    C2 = params.C2;
    C3 = params.C3;
    D = params.D;
    p1 = params.p1;
    C = params.C;
    m = params.m;
    a = params.a;
    n = params.n;
    
    % 复制动态方程
    dx_dt = x * (1 - x) * ((C2 * D * y * z) + D - D * y * z - C3 + C2 * z * y);
    dy_dt = y * (1 - y) * (k1 * p1 * z + C * z + m * x + a * n - n - C);
    dz_dt = z * (1 - z) * (a * n + C * y + m * x + k2 * p1 * y - n - C);
    
    dstate_dt = [dx_dt; dy_dt; dz_dt];
end

function rgb = hex2rgb(hex_color)
    % 辅助函数：将颜色字符转换为RGB值
    switch hex_color
        case 'r'
            rgb = [1, 0, 0];
        case 'g'
            rgb = [0, 1, 0];
        case 'b'
            rgb = [0, 0, 1];
        otherwise
            rgb = [0.5, 0.5, 0.5];
    end
end