如何用Matlab的interp1实现准确外推，获得无断点的光滑B(H)曲线？

我需要优化B(H)曲线，得到一条中间无断点的光滑曲线。当前曲线的断点是由外推至10⁴ A/m（此处斜率需等于μ₀）的操作导致的，尝试过用interp1在0至Hi(end)区间内插、在Hi(end)到近10⁴ H/m区间用spline方法外推，但始终无法得到光滑的外推效果，寻求解决方案。

原始实现代码：

mu0=4*pi*1e-7; % perméabilité du vide 
nu0=1/mu0; % réluctivité du vide 

H_data_ini = [ 0 0.3979 0.7958 1.1937 1.5916 1.9895 2.3874 2.7853 3.1832 3.9790 4.7748 5.9685 7.9580 9.9475 11.9370 15.9160 19.8950 23.8740 31.8320 39.7900 47.7480 59.6850 79.5800 99.4750 119.3700 159.1600 198.9500 238.7400 318.3200 397.9000 477.4800 636.6400 795.8000]

B_data_ini = [ 0    0.0030    0.0090    0.0170    0.0300    0.0520    0.0800    0.1110    0.1360    0.1910    0.2260    0.2660    0.3100    0.3420    0.3670    0.4030    0.4300    0.4520    0.4860    0.5130    0.5340    0.5610    0.5930 0.6190    0.6390    0.6700    0.6920    0.7100    0.7380    0.7550    0.7700    0.7940    0.8060 ]

Bi=0:0.0001:B_data_ini(end);
Hi = interp1(B_data_ini,H_data_ini,Bi);

Hii=B_data_ini(end)+(500:500:1e4);
Bii = interp1(Hi,Bi,Hii,'spline','extrap'); 

% Phytherm260
M_s = Bii(end) - mu0 * Hii(end); 
H_end=Hii(end)+(500:50:2e4);
B_end=mu0.*H_end+M_s;



%% Courbe B(H) expérimental

% Phytherm260 
H_data=[Hi Hii H_end]; 
B_data=[Bi Bii B_end]; 

figure; 
plot(H_data,B_data,'k','LineWidth',2, 'DisplayName', 'T = 20^°C');

解决方案 #

核心思路：基于磁化强度M(H)的连续建模 #

B(H)的物理本质是B = μ₀(H + M)，其中M为磁化强度。当H足够大时M趋近于饱和值Ms，此时B(H)的斜率自然等于μ₀。直接对B(H)分段外推易产生断点，转而对M(H)进行连续拟合+外推，能从根源保证曲线光滑过渡。

具体实现方法 #

1. 物理拟合方案（推荐）

利用铁磁材料的饱和磁化特性，用朗之万或反正切函数拟合M(H)，再推导B(H)：

mu0=4*pi*1e-7;

% 1. 从原始数据计算磁化强度M
M_data_ini = B_data_ini ./ mu0 - H_data_ini;

% 2. 定义拟合模型：M(H) = Ms * tanh(H/H0)（天然带饱和特性）
model = @(params, H) params(1) * tanh(H / params(2));

% 3. 初始参数猜测（根据实验数据调整）
Ms_guess = max(M_data_ini);
H0_guess = 1000;
initial_params = [Ms_guess, H0_guess];

% 4. 拟合数据
fit_params = lsqcurvefit(model, initial_params, H_data_ini, M_data_ini);
Ms_fit = fit_params(1);
H0_fit = fit_params(2);

% 5. 生成全程光滑的B(H)曲线
H_full = linspace(0, 2e4, 10000); % 足够多的采样点保证光滑
M_full = model(fit_params, H_full);
B_full = mu0 * (H_full + M_full);

% 6. 绘图验证
figure;
plot(H_full, B_full, 'k', 'LineWidth', 2, 'DisplayName', '光滑拟合曲线');
hold on;
plot(H_data_ini, B_data_ini, 'ro', 'MarkerSize', 4, 'DisplayName', '原始实验点');
legend;
xlabel('H (A/m)');
ylabel('B (T)');

2. 分段光滑过渡方案（无需拟合）

如果不想用拟合，可通过保证外推段与内插段的斜率连续来消除断点：

mu0=4*pi*1e-7;

% 内插段保持原有逻辑
Bi=0:0.0001:B_data_ini(end);
Hi = interp1(B_data_ini,H_data_ini,Bi);

% 计算内插段终点的斜率（dB/dH）
idx_end = find(Hi >= Hi(end)-100, 1, 'first');
slope_end = (Bi(end) - Bi(idx_end)) / (Hi(end) - Hi(idx_end));

% 外推到1e4 A/m的区间，用斜率线性过渡的方式生成B值
Hii_new = linspace(Hi(end), 1e4, 500);
Bii_new = zeros(size(Hii_new));
Bii_new(1) = Bi(end);

for i=2:length(Hii_new)
    dH = Hii_new(i) - Hii_new(i-1);
    % 斜率从内插段终点的slope_end平滑过渡到mu0
    current_slope = slope_end + (mu0 - slope_end)*(Hii_new(i) - Hi(end))/(1e4 - Hi(end));
    Bii_new(i) = Bii_new(i-1) + current_slope * dH;
end

% 饱和段保持原有逻辑
M_s = Bii_new(end) - mu0 * Hii_new(end);
H_end=linspace(1e4, 2e4, 500);
B_end=mu0.*H_end+M_s;

% 合并数据并绘图
H_data=[Hi Hii_new H_end]; 
B_data=[Bi Bii_new B_end]; 

figure; 
plot(H_data,B_data,'k','LineWidth',2, 'DisplayName', 'T = 20^°C');

方案优势 #

• 物理拟合方案：从铁磁材料的磁化规律出发，曲线全程光滑且符合物理特性，无断点风险
• 分段过渡方案：无需拟合，仅通过斜率连续过渡实现光滑外推，适合快速调整

内容的提问来源于stack exchange，提问作者hakim