Matplotlib绘制B-H磁滞回线时线性拟合斜率方程异常问题求助

问题描述 #

用Matplotlib绘制磁滞回线时，拟合的斜率方程显示结果趋近于0，但图表本身显示正常。完整代码如下：

import matplotlib.pyplot as plt
import numpy as np

# Section 1: Calculating values of H for toroid A

# defining the values to find H for toroid A
N_1 = 3420
L_1 = 0.076
R_1 = 8

# defining function to find the value of H with given V_1
def H_toroid_A(V_1):
    return (V_1 * N_1) / (L_1 * R_1)

# Section 2: Calculating the delta list of uncertainties for H toroid A

# defining variables and uncertainties
V_1_values = [0.198, 0.266, 0.340, 0.402, 0.505]
delta_V_1_values = [0.050, 0.050, 0.050, 0.050, 0.1]  # uncertainties for V_1
delta_N_1 = 0  # uncertainty for N_1
delta_L_1 = 0.001  # uncertainty for L_1
delta_R_1 = 0.4  # uncertainty for R_1

# function to calculate uncertainties in H
def calculate_uncertainty_H(V_values, delta_V_values, delta_N, delta_L, delta_R):
    uncertainties = []
    for i in range(len(V_values)):
        uncertainty_H = np.sqrt(
            ((V_values[i] / (L_1 * R_1)) * delta_N)**2 +
            ((N_1 / (L_1 * R_1)) * delta_V_values[i])**2 +
            ((-V_values[i] * N_1 / (L_1**2 * R_1)) * delta_L)**2 +
            ((-V_values[i] * N_1 / (L_1 * R_1**2)) * delta_R)**2
        )
        uncertainties.append(uncertainty_H)
    return uncertainties

uncertainties_H = calculate_uncertainty_H(V_1_values, delta_V_values, delta_N_1, delta_L_1, delta_R_1)

# Section 3: Calculating values of B for toroid A

# defining the values to find B for toroid A
R_2 = 270000
C = 220e-9
N_2 = 820
A_2 = 9.025e-5

# defining function to find the value of B with given V_c
def B_toroid_A(V_c):
    return (R_2 * C * V_c) / (N_2 * A_2)

# Section 4: Calculating list of uncertainties for B toroid A

# defining variables and uncertainties
V_c_values = [0.0266, 0.0372, 0.0505, 0.0605, 0.076]
delta_V_c_values = [0.005, 0.005, 0.01, 0.01, 0.01]  # uncertainties for V_c
delta_N_2 = 0  # uncertainty for N_2
delta_C = 1.1e-8  # uncertainty for C
delta_A_2 = 5e-8  # uncertainty for A_2
delta_R_2 = 13500  # uncertainty for R_2

# function to calculate uncertainties in B
def calculate_uncertainty_B(V_values, delta_V_values, delta_N, delta_C, delta_A, delta_R):
    uncertainties = []
    for i in range(len(V_values)):
        uncertainty_B = np.sqrt(
            ((C * V_values[i]) / (N_2 * A_2))**2 * delta_R_2**2 +
            ((R_2 * V_values[i]) / (N_2 * A_2))**2 * delta_C**2 +
            ((R_2 * C) / (N_2 * A_2))**2 * delta_V_values[i]**2 +
            ((R_2 * C * V_values[i]) / (N_2 * A_2**2))**2 * delta_A_2**2 +
            ((R_2 * C * V_values[i]) / (N_2**2 * A_2))**2 * delta_N_2**2
        )
        uncertainties.append(uncertainty_B)
    return uncertainties

uncertainties_B = calculate_uncertainty_B(V_c_values, delta_V_c_values, delta_N_2, delta_C, delta_A_2, delta_R_2)

# Section Alpha: calculating list of values of H
H_values = [H_toroid_A(V_1) for V_1 in V_1_values]

# Section Beta: calculating list of values of B
B_values = [B_toroid_A(V_c) for V_c in V_c_values]

# Plotting
plt.figure(figsize=(8, 6))
plt.errorbar(H_values, B_values, xerr=uncertainties_H, yerr=uncertainties_B, fmt='o', label='Data Points')
plt.xlabel('H (Intensité du champ magnétique)')
plt.ylabel('B (Densité du flux magnétique)')
plt.title('Courbe Hysteresis')

# Performing linear regression with uncertainties
coefficients, cov_matrix = np.polyfit(H_values, B_values, 1, w=1/np.array(uncertainties_B), cov=True)
slope = coefficients[0]
intercept = coefficients[1]

# Plotting the linear fit
plt.plot(H_values, slope * np.array(H_values) + intercept, color='red', label=f'Linear Fit: B = {slope:.2f} * H + {intercept:.2f}')
plt.legend()

plt.grid(True)
plt.show()

错误原因分析 #

1. 数值格式化精度不足
   计算H和B的实际数量级：

   • H的取值范围约为1100~2800 A/m
   • B的取值范围约为0.02~0.06 T
     拟合斜率为ΔB/ΔH，数值约为2e-5量级，用.2f格式化时会被四舍五入为0.00，这是显示问题，并非拟合结果错误。

2. 拟合本身无逻辑错误
   代码中的np.polyfit调用逻辑正确，权重参数w=1/uncertainties_B符合加权线性回归的要求，拟合结果本身是准确的，只是显示时精度不够导致误解。

修复方案 #

• 调整格式化精度：将拟合方程的格式化字符串从.2f改为更高精度的格式，比如.6f（斜率）和.4f（截距），修改后的代码行如下：

  plt.plot(H_values, slope * np.array(H_values) + intercept, color='red', label=f'线性拟合: B = {slope:.6f} * H + {intercept:.4f}')
  

• 验证原始拟合值：在代码中添加print("拟合斜率：", slope)，直接查看斜率的原始数值，确认拟合结果的正确性。

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