【技巧篇】从实验数据到数学表达式：应力-应变曲线拟合的最佳实践

我们经常会需要做这样一件事：做完试验得到试验曲线了，需要通过拟合这些试验曲线得到一个数学表达式，那么以后别人在做相似的工作时，就不用重新做实验了，直接用我们给出的数学表达式就好啦。

而我们常常是做完试验，得到一大堆数据了。同时也通过机理分析或者看文献大概知道这类试验的曲线表达式大概是什么样子了。

但是表达式中的几个未知系数怎么也确定不下来，这篇就是介绍怎么通过一段matlab小程序解决这个问题。

试验数据

试验得到的应力-应变曲线如下：

 

数学表达式

拟合程序

别人的表达式里面的系数可能是拿公式算出来的，可能是拟合出来的，反正肯定是不适合我们自己的情况的，那这时候就要个性化地拟合出适合于我们自己试验结果的系数了。

close all
clear 
clc  
initial_data=[]%%我们的试验数据
A=initial_data
%% part3 截掉不用的数据 
A(:,2)=-A(:,2)
plot(A(:,1),A(:,2),'s');
%% 换成sigma-epsilon曲线

L0=200
A(:,1)=A(:,1)/L0
plot(A(:,1),A(:,2),'s')
epsilon=A(:,1)
sigma_test=A(:,2)
%% 换成sigma-X曲线
A1=A
col_LVDT=1
col_P=2
[Y_col,Ind_row]=max(A);
P_max1=Y_col(col_P);            
P_max1_hang=Ind_row(col_P);      
LVDT_max=A1(P_max1_hang,col_LVDT)
A1(:,1)=A1(:,1)/LVDT_max
plot(A1(:,1),A1(:,2),'s')
x=A1(:,1)
new=[epsilon,x,sigma_test];
% 绘图
figure;
plot(new(:,1), new(:,3), 's', 'Color', [1 0.75 0.8], 'LineWidth', 2);
% 设置坐标轴标签，字体为新罗马，字号为14
xlabel('strain', 'FontName', 'Times New Roman', 'FontSize', 14);
ylabel('stress-MPa', 'FontName', 'Times New Roman', 'FontSize', 14);
% 设置刻度线字体大小为 14
ax = gca;  % 获取当前坐标轴句柄
ax.FontSize = 14;  % 设置刻度线字体大小为 25
ax.TickLength = [0.02, 0.02];  % 设置刻度线长度

%% 筛点
downsampleFactor = 10; % 下采样因子，每隔10个点取一个值
data_downsampled = new(1:downsampleFactor:end, :);
plot(data_downsampled (:,1),data_downsampled (:,3),'s') 
plot(data_downsampled (:,2),data_downsampled (:,3),'s') 

%% model函数拟合
epsilon = new(:, 1);
x = new(:, 2);
sigma_test = new(:, 3);
% 分别获取两部分数据
data_part1 = new(x <= 1, :);
data_part2 = new(x > 1, :);
% 对于 x <= 1 部分数据的拟合
epsilon1 = data_part1(:, 1);
x1 = data_part1(:, 2);
sigma_test1 = data_part1(:, 3);
plot(x1,sigma_test1,'s') %%筛点后的p-delta图
% 定义模型函数
model1 = @(p, x1, epsilon1) (p(1) * p(2) * epsilon1) ./ (p(2) - 1 + x1.^p(2));
% 初始猜测值 [aa, n]
aa1_intial=42.95
n=1.05
initial_guess1 = [aa1_intial, n];
% 进行非线性拟合
params1 = lsqcurvefit(@(p, x1) model1(p, x1, epsilon1), initial_guess1, x1, sigma_test1);
% 对于 x > 1 部分数据的拟合
epsilon2 = data_part2(:, 1);
x2 = data_part2(:, 2);
sigma_test2 = data_part2(:, 3);
plot(x2,sigma_test2,'s') %%筛点后的p-delta图
% 定义模型函数
model2 = @(p, x2, epsilon2) (p(1) * epsilon2) ./ (p(2) * (x2 - 1).^2 + x2);
% 初始猜测值 [aa, bb]
initial_guess2 = [42.95, 1.57];
% 进行非线性拟合
params2 = lsqcurvefit(@(p, x2) model2(p, x2, epsilon2), initial_guess2, x2, sigma_test2);
% 显示结果
aa1 = params1(1);
n = params1(2);
aa2 = params2(1);
bb = params2(2);
fprintf('拟合结果:\naa1 = %.4f\nn = %.4f\naa2 = %.4f\nbb = %.4f\n', aa1, n, aa2, bb);
%% 公式值与计算值对比
%上升段
sigma = zeros(size(epsilon1));
for i = 1:length(epsilon1)
    sigma(i) = (aa1 * n * epsilon1(i)) / (n - 1 + x1(i)^n);
end
plot(epsilon1,  sigma,'s')
hold on 
plot(epsilon1,  sigma_test1,'s')
A_UP=[epsilon1,sigma];
%下降段
% 初始化 sigma_down 数组
sigma_down = zeros(size(epsilon2));
% 计算 sigma_down
for i = 1:length(epsilon2)
    sigma_down(i) = (aa2 * epsilon2(i)) / (bb * (x2(i) - 1)^2 + x2(i));
end
% sigma_down 是一个单列数据
plot(epsilon2,  sigma_down,'s')
hold on 
plot(epsilon2,  sigma_test2,'s')

%% 拟合的结果
% 合并 epsilon 列
epsilon_combined = [epsilon1; epsilon2];
% 合并 sigma 列
sigma_combined = [sigma; sigma_down];
% 合并 x 列
x_combined = [x1; x2];
% 形成最终的 end_RRAY
end_RRAY = [epsilon_combined, x_combined, sigma_combined];
plot(epsilon_combined, sigma_combined,'s')
hold on 
plot(epsilon,  sigma_test,'s')
plot(x_combined, sigma_combined,'s')
hold on 
plot(x,  sigma_test,'s')
% 设置坐标轴标签，字体为新罗马，字号为14
xlabel('strain', 'FontName', 'Times New Roman', 'FontSize', 14);
ylabel('stress-MPa', 'FontName', 'Times New Roman', 'FontSize', 14);
% 设置刻度线字体大小为 14
ax = gca;  % 获取当前坐标轴句柄
ax.FontSize = 14;  % 设置刻度线字体大小为 25
ax.TickLength = [0.02, 0.02];  % 设置刻度线长度

拟合结果

蓝线为拟合曲线，红线为试验曲线