% 步骤1：导入数据
t = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0];
y = [5.8955, 3.5639, 2.5173, 1.9790, 1.8990, 1.3938, 1.1359, 1.0096, 1.0343, 0.8435, 0.6856, 0.6100, 0.5392, 0.3946, 0.3903, 0.5474, 0.3459, 0.1730, 0.2211, 0.1704, 0.2636];
% 步骤2：定义双指数衰减模型
doubleExp = @(p, t) p(1)*exp(-p(2)*t) + p(3)*exp(-p(4)*t) + p(5);
% 步骤3：初始参数估计（关键步骤）
% C0: 取最后5个点的平均值作为本底估计
C0 = mean(y(end-4:end));
% 快衰减分量：取前5个点拟合单指数
early_t = t(1:5);
early_y = y(1:5) - C0;
p_early = polyfit(early_t, log(early_y), 1);
lambda1_0 = -p_early(1); % 快衰减常数
A1_0 = exp(p_early(2)); % 快衰减振幅
% 慢衰减分量：取中间段数据
mid_idx = (t >= 0.5) & (t <= 1.5);
mid_t = t(mid_idx);
mid_y = y(mid_idx) - A1_0*exp(-lambda1_0*mid_t) - C0;
p_mid = polyfit(mid_t, log(mid_y), 1);
lambda2_0 = -p_mid(1); % 慢衰减常数
A2_0 = exp(p_mid(2)); % 慢衰减振幅
% 初始参数向量 [A1, λ1, A2, λ2, C]
p0 = [A1_0, lambda1_0, A2_0, lambda2_0, C0];
% 步骤4：执行非线性最小二乘拟合
options = optimoptions('lsqcurvefit', 'Display', 'iter', 'MaxFunctionEvaluations', 2000);
lb = [0, 0, 0, 0, 0]; % 参数下界（所有参数非负）
ub = [10, 10, 10, 1, 1]; % 参数上界
[params, resnorm, residual, exitflag] = lsqcurvefit(doubleExp, p0, t, y, lb, ub, options);
% 步骤5：提取拟合参数
A1_fit = params(1);
lambda1_fit = params(2);
A2_fit = params(3);
lambda2_fit = params(4);
C_fit = params(5);
% 步骤6：计算半衰期
T_half1 = log(2)/lambda1_fit; % 快衰减半衰期
T_half2 = log(2)/lambda2_fit; % 慢衰减半衰期
% 步骤7：可视化结果
t_fit = linspace(min(t), max(t), 500);
y_fit = doubleExp(params, t_fit);
figure;
scatter(t, y, 40, 'filled', 'MarkerFaceColor', [0 0.447 0.741]);
hold on;
plot(t_fit, y_fit, 'LineWidth', 2, 'Color', [0.85 0.325 0.098]);
xlabel('时间');
ylabel('浓度');
title('放射性衰减双指数拟合');
legend('观测数据', '拟合曲线', 'Location', 'best');
grid on;
% 标注拟合参数
text(0.5, 4, sprintf('y = %.2f e^{-%.3f t} + %.2f e^{-%.3f t} + %.2f', ...
A1_fit, lambda1_fit, A2_fit, lambda2_fit, C_fit), ...
'FontSize', 10, 'BackgroundColor', [1 1 1 0.7]);
% 步骤8：评估拟合质量
SSE = resnorm; % 残差平方和
SST = sum((y - mean(y)).^2);
R_squared = 1 - SSE/SST;
fprintf('拟合参数:\n');
fprintf('A1 = %.4f, λ1 = %.4f (半衰期 = %.4f)\n', A1_fit, lambda1_fit, T_half1);
fprintf('A2 = %.4f, λ2 = %.4f (半衰期 = %.4f)\n', A2_fit, lambda2_fit, T_half2);
fprintf('C = %.4f\n', C_fit);
fprintf('拟合优度: R² = %.4f\n', R_squared);
生成的数据为什么只有19个 Iteration Func-count Resnorm step optimality
0 6 1.2679 4.6
1 12 0.992064 0.202205 0.955
2 18 0.972854 0.93646 7.42
3 24 0.630698 0.234115 0.901
4 30 0.559772 0.109448 1.5
5 36 0.502012 0.118323 1.04
6 42 0.498087 0.0088884 0.553
7 48 0.385485 0.46823 0.455
8 54 0.370825 0.056662 0.419
9 60 0.369192 0.00776557 0.326
10 66 0.344053 0.108125 0.257
11 72 0.316799 0.140779 0.176
12 78 0.294731 0.172189 0.099
13 84 0.283437 0.162168 0.0448
14 90 0.280038 0.110795 0.0162
15 96 0.279442 0.0549579 0.00453
16 102 0.279385 0.0185708 0.00106
17 108 0.279381 0.00479553 0.000236
18 114 0.279381 0.00110347 5.21e-05
最新发布
扫一扫
7193
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
