clear all
clc
close all
%%%%%%%%%%%%%%%%%%%%%%% 基础参数
a = 1;
h = 0.2756;
b = 0.5;
K = 5;
r = 1;
c = 0.5;
d1 = 0.5;
d2 = 0.5;
tau = 0;
L = 4*pi;
dx = 0.1;
n = round(L/dx);
NX = linspace(0, L, n);
Tend = 100;
dt = 0.001;
%%%%%%%%%%%%%%%%%%%%%%% 扩展参数敏感性分析设置
% 原始参数 + beta, m1, m2
parameters = {'r', 'K', 'h', 'b', 'c', 'd1', 'd2', 'beta', 'm1', 'm2'};
base_values = [1, 5, 0.2756, 0.5, 0.5, 0.5, 0.5, 0.15, 0.5, 0.3];
variations = 0.5:0.1:1.5; % 参数变化范围 (50% 到 150%)
% 存储结果
results = struct();
%%%%%%%%%%%%%%%%%%%%%%% 进行参数敏感性分析
for param_idx = 1:length(parameters)
param_name = parameters{param_idx};
fprintf('正在分析参数 %s 的敏感性...\n', param_name);
mean_u_final = [];
mean_v_final = [];
amp_u_final = [];
amp_v_final = [];
for var_idx = 1:length(variations)
% 设置当前参数值
current_params = base_values;
current_params(param_idx) = base_values(param_idx) * variations(var_idx);
% 提取参数
r_temp = current_params(1);
K_temp = current_params(2);
h_temp = current_params(3);
b_temp = current_params(4);
c_temp = current_params(5);
d1_temp = current_params(6);
d2_temp = current_params(7);
beta_base_temp = current_params(8);
m1_base_temp = current_params(9);
m2_base_temp = current_params(10);
% 运行模拟
[mean_u, mean_v, amp_u, amp_v] = run_simulation_extended(...
r_temp, K_temp, h_temp, b_temp, c_temp, d1_temp, d2_temp, a, ...
beta_base_temp, m1_base_temp, m2_base_temp, L, dx, Tend, dt);
mean_u_final(var_idx) = mean_u;
mean_v_final(var_idx) = mean_v;
amp_u_final(var_idx) = amp_u;
amp_v_final(var_idx) = amp_v;
end
% 存储结果
results.(param_name).variations = variations;
results.(param_name).mean_u = mean_u_final;
results.(param_name).mean_v = mean_v_final;
results.(param_name).amp_u = amp_u_final;
results.(param_name).amp_v = amp_v_final;
end
%%%%%%%%%%%%%%%%%%%%%%% 绘制完整参数敏感性分析图
figure('Position', [100, 100, 1600, 1200]);
% 绘制平均密度敏感性
subplot(2,2,1)
hold on
colors = lines(length(parameters));
for i = 1:length(parameters)
param_name = parameters{i};
plot(results.(param_name).variations, results.(param_name).mean_u, ...
'o-', 'LineWidth', 2, 'Color', colors(i,:), 'MarkerSize', 5, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('猎物平均密度')
title('参数对猎物平均密度的敏感性')
legend('Location', 'bestoutside')
grid on
subplot(2,2,2)
hold on
for i = 1:length(parameters)
param_name = parameters{i};
plot(results.(param_name).variations, results.(param_name).mean_v, ...
'o-', 'LineWidth', 2, 'Color', colors(i,:), 'MarkerSize', 5, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('捕食者平均密度')
title('参数对捕食者平均密度的敏感性')
legend('Location', 'bestoutside')
grid on
% 绘制振幅敏感性
subplot(2,2,3)
hold on
for i = 1:length(parameters)
param_name = parameters{i};
plot(results.(param_name).variations, results.(param_name).amp_u, ...
's-', 'LineWidth', 2, 'Color', colors(i,:), 'MarkerSize', 5, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('猎物空间振幅')
title('参数对猎物空间振幅的敏感性')
legend('Location', 'bestoutside')
grid on
subplot(2,2,4)
hold on
for i = 1:length(parameters)
param_name = parameters{i};
plot(results.(param_name).variations, results.(param_name).amp_v, ...
's-', 'LineWidth', 2, 'Color', colors(i,:), 'MarkerSize', 5, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('捕食者空间振幅')
title('参数对捕食者空间振幅的敏感性')
legend('Location', 'bestoutside')
grid on
sgtitle('完整参数敏感性分析（包含beta、m1、m2）', 'FontSize', 16, 'FontWeight', 'bold');
%%%%%%%%%%%%%%%%%%%%%%% 绘制敏感性系数图
figure('Position', [100, 100, 1400, 600]);
% 计算敏感性系数
sensitivity_u = zeros(length(parameters), 1);
sensitivity_v = zeros(length(parameters), 1);
for i = 1:length(parameters)
param_name = parameters{i};
mean_u = results.(param_name).mean_u;
mean_v = results.(param_name).mean_v;
variations = results.(param_name).variations;
% 在基准值处的敏感性
base_idx = find(variations == 1, 1);
if base_idx > 1 && base_idx < length(variations)
sensitivity_u(i) = (mean_u(base_idx+1) - mean_u(base_idx-1)) / (variations(base_idx+1) - variations(base_idx-1));
sensitivity_v(i) = (mean_v(base_idx+1) - mean_v(base_idx-1)) / (variations(base_idx+1) - variations(base_idx-1));
else
sensitivity_u(i) = (mean_u(end) - mean_u(1)) / (variations(end) - variations(1));
sensitivity_v(i) = (mean_v(end) - mean_v(1)) / (variations(end) - variations(1));
end
end
% 绘制敏感性系数条形图
subplot(1,2,1)
bar(sensitivity_u)
set(gca, 'XTickLabel', parameters, 'XTickLabelRotation', 45)
ylabel('敏感性系数')
title('猎物平均密度敏感性系数')
grid on
subplot(1,2,2)
bar(sensitivity_v)
set(gca, 'XTickLabel', parameters, 'XTickLabelRotation', 45)
ylabel('敏感性系数')
title('捕食者平均密度敏感性系数')
grid on
sgtitle('完整参数敏感性系数比较', 'FontSize', 16, 'FontWeight', 'bold');
%%%%%%%%%%%%%%%%%%%%%%% 绘制beta、m1、m2的详细分析
figure('Position', [100, 100, 1200, 800]);
% 选择要详细分析的参数
special_params = {'beta', 'm1', 'm2'};
colors_special = [1 0 0; 0 0.5 0; 0 0 1]; % 红、绿、蓝
subplot(2,3,1)
hold on
for i = 1:length(special_params)
param_name = special_params{i};
plot(results.(param_name).variations, results.(param_name).mean_u, ...
'o-', 'LineWidth', 2.5, 'Color', colors_special(i,:), 'MarkerSize', 6, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('猎物平均密度')
title('beta/m1/m2对猎物密度的影响')
legend('Location', 'best')
grid on
subplot(2,3,2)
hold on
for i = 1:length(special_params)
param_name = special_params{i};
plot(results.(param_name).variations, results.(param_name).mean_v, ...
'o-', 'LineWidth', 2.5, 'Color', colors_special(i,:), 'MarkerSize', 6, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('捕食者平均密度')
title('beta/m1/m2对捕食者密度的影响')
legend('Location', 'best')
grid on
subplot(2,3,3)
hold on
for i = 1:length(special_params)
param_name = special_params{i};
plot(results.(param_name).variations, results.(param_name).amp_u, ...
's-', 'LineWidth', 2.5, 'Color', colors_special(i,:), 'MarkerSize', 6, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('猎物空间振幅')
title('beta/m1/m2对猎物振幅的影响')
legend('Location', 'best')
grid on
subplot(2,3,4)
hold on
for i = 1:length(special_params)
param_name = special_params{i};
plot(results.(param_name).variations, results.(param_name).amp_v, ...
's-', 'LineWidth', 2.5, 'Color', colors_special(i,:), 'MarkerSize', 6, ...
'DisplayName', param_name);
end
xlabel('参数变化倍数')
ylabel('捕食者空间振幅')
title('beta/m1/m2对捕食者振幅的影响')
legend('Location', 'best')
grid on
% 绘制相图 (beta vs m1)
subplot(2,3,5)
beta_var = results.beta.variations;
m1_var = results.m1.variations;
[BB, MM] = meshgrid(beta_var, m1_var);
Z_phase = zeros(length(beta_var), length(m1_var));
% 计算相图数据
for i = 1:length(beta_var)
for j = 1:length(m1_var)
% 使用线性插值近似
beta_val = base_values(8) * beta_var(i);
m1_val = base_values(9) * m1_var(j);
% 简化的稳定性判断（基于线性稳定性分析）
% 这里使用一个简化的标准：捕食者密度是否能够维持
if beta_val > 0.1 && m1_val < 0.8
Z_phase(j,i) = 1; % 共存状态
elseif beta_val > 0.2
Z_phase(j,i) = 2; % 捕食者占优
else
Z_phase(j,i) = 3; % 猎物占优或灭绝
end
end
end
imagesc(beta_var, m1_var, Z_phase)
xlabel('beta变化倍数')
ylabel('m1变化倍数')
title('beta-m1参数相图')
colorbar
colormap(jet(3))
% 绘制三维敏感性图
subplot(2,3,6)
% 选择两个参数进行三维可视化
[X, Y] = meshgrid(results.beta.variations, results.m2.variations);
Z = zeros(size(X));
for i = 1:size(X,1)
for j = 1:size(X,2)
beta_mult = X(i,j);
m2_mult = Y(i,j);
% 查找对应的结果（简化处理）
beta_idx = find(results.beta.variations == beta_mult);
m2_idx = find(results.m2.variations == m2_mult);
if ~isempty(beta_idx) && ~isempty(m2_idx)
Z(i,j) = results.beta.mean_v(beta_idx) + results.m2.mean_v(m2_idx);
end
end
end
surf(X, Y, Z)
xlabel('beta变化倍数')
ylabel('m2变化倍数')
zlabel('捕食者密度')
title('beta-m2对捕食者密度的联合影响')
colorbar
sgtitle('beta、m1、m2参数详细分析', 'FontSize', 16, 'FontWeight', 'bold');
%%%%%%%%%%%%%%%%%%%%%%% 结果汇总
fprintf('\n=== 完整敏感性分析结果汇总 ===\n');
[~, max_sense_u] = max(abs(sensitivity_u));
[~, max_sense_v] = max(abs(sensitivity_v));
fprintf('对猎物密度最敏感的参数: %s (敏感性系数: %.4f)\n', parameters{max_sense_u}, sensitivity_u(max_sense_u));
fprintf('对捕食者密度最敏感的参数: %s (敏感性系数: %.4f)\n', parameters{max_sense_v}, sensitivity_v(max_sense_v));
% 特别关注beta, m1, m2的结果
fprintf('\n--- beta, m1, m2 敏感性分析 ---\n');
beta_idx = find(strcmp(parameters, 'beta'));
m1_idx = find(strcmp(parameters, 'm1'));
m2_idx = find(strcmp(parameters, 'm2'));
fprintf('beta 敏感性 - 猎物: %.4f, 捕食者: %.4f\n', sensitivity_u(beta_idx), sensitivity_v(beta_idx));
fprintf('m1 敏感性 - 猎物: %.4f, 捕食者: %.4f\n', sensitivity_u(m1_idx), sensitivity_v(m1_idx));
fprintf('m2 敏感性 - 猎物: %.4f, 捕食者: %.4f\n', sensitivity_u(m2_idx), sensitivity_v(m2_idx));
%%%%%%%%%%%%%%%%%%%%%%% 扩展的模拟函数
function [mean_u_final, mean_v_final, amp_u_final, amp_v_final] = run_simulation_extended(...
r, K, h, b, c, d1, d2, a, beta_base, m1_base, m2_base, L, dx, Tend, dt)
% 空间划分格数
n = round(L/dx);
NX = linspace(0, L, n);
% 扩散系数
r1 = d1*dt/dx^2;
r2 = d2*dt/dx^2;
% 边界条件索引
ind = 2:(n-1);
IndR = 3:n;
IndL = 1:(n-2);
% 时间步数
Nstep = round(Tend/dt);
% 初始条件 - 使用传入的beta_base, m1_base, m2_base
Eq = [0.6, 0.1, beta_base, m1_base, m2_base]; % 注意顺序调整
U1 = zeros(n, 1);
U2 = zeros(n, 1);
for i = 1:n
U1(i) = Eq(1) + 0.01 * sin(i*dx);
U2(i) = Eq(2) + 0.01 * cos(i*dx);
end
% 使用传入的基准值生成空间变化的参数
beta = beta_base + 0.05 * cos(NX);
m1 = m1_base + 0.01 * cos(NX);
m2 = m2_base + 0.005 * cos(NX);
% 主循环
for p = 1:Nstep
% 扩散项
DU1 = zeros(n,1);
DU1(ind) = U1(IndR) - 2*U1(ind) + U1(IndL);
DU1(1) = 2*(U1(2) - U1(1));
DU1(n) = 2*(U1(n-1) - U1(n));
DU2 = zeros(n,1);
DU2(ind) = U2(IndR) - 2*U2(ind) + U2(IndL);
DU2(1) = 2*(U2(2) - U2(1)) - 2*c*dx*U2(1);
DU2(n) = 2*(U2(n-1) - U2(n)) - 2*c*dx*U2(n);
% 反应项
f = r * U1(ind) .* (1 - U1(ind)./K) - (beta(ind)' .* U1(ind) .* U2(ind)) ./ (a + U1(ind)) - h * U1(ind);
g = m1(ind)' .* U2(ind) + (b * beta(ind)' .* U1(ind) .* U2(ind)) ./ (a + U1(ind)) - m2(ind)' .* U2(ind) .* U2(ind);
% 更新
U1_new = zeros(n,1);
U2_new = zeros(n,1);
U1_new(ind) = U1(ind) + r1 * DU1(ind) + dt * f;
U2_new(ind) = U2(ind) + r2 * DU2(ind) + dt * g;
% 边界条件
U1_new(1) = U1_new(2);
U1_new(n) = U1_new(n-1);
U2_new(1) = U2_new(2) / (1 + c*dx);
U2_new(n) = U2_new(n-1) / (1 + c*dx);
% 非负性
U1_new = max(U1_new, 0);
U2_new = max(U2_new, 0);
U1 = U1_new;
U2 = U2_new;
% 稳定性检查
if mod(p, 1000) == 0
if any(isnan(U1)) || any(isnan(U2)) || max(U1) > 100 || max(U2) > 100
warning('模拟不稳定，参数: r=%.2f, beta=%.3f, m1=%.3f', r, beta_base, m1_base);
break;
end
end
end
% 计算统计量
mean_u_final = mean(U1);
mean_v_final = mean(U2);
amp_u_final = max(U1) - min(U1);
amp_v_final = max(U2) - min(U2);
end检查该代码参数敏感性分析是否正确？
本文提供了解决Windows环境下使用Cocos2d-lib开发新项目时遇到的致命错误LNK1104的方法。通过将路径cocos2dx/Debug.win32加入到配置属性中，可以成功编译项目。
扫一扫
796
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
