求助：用高斯消元法求解液体体积公式V=1+aT+bT²系数的Python代码

实验数据 #

V (cm³) T (°C)
1.032 10.0
1.094 29.5
1.156 50.0
1.215 69.5
1.273 90.0

需求说明 #

已知液体体积与温度满足关系式 V = 1 + aT + bT²，需用Python编写程序，必须通过高斯消元法求解系数a和b。以下是目前编写的代码，需要修正：

原代码问题分析 #

原代码存在几个关键错误：

• 构造的矩阵不符合最小二乘的高斯消元要求：数据点数量（5个）远多于待求系数（2个），需将原方程转化为正规方程组再求解
• 错误构造高维度矩阵并添加无效行，导致矩阵维度不匹配，高斯消元无法正常运行
• 原公式中常数项固定为1，无需求解，但原代码仍将其作为待求系数处理

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
import os, sys

def GEshow(A, b, ptol, verbose):
    eps = sys.float_info.epsilon
    peps = max(ptol, 5 * eps)                       
    print(f"选主元的容差值为 {peps:e}")
    row, column = np.shape(A)
    if row != column:
        print('矩阵A必须是方阵！程序退出。')
        os._exit(os.EX_OK)
    if verbose:
        print('高斯消元增广矩阵：')
        print(np.concatenate((A, b.reshape(-1, 1)), axis=1))              

    for k in range(column):                                  
        print(f"--- 消元阶段 {k}")
        pivot = A[k, k]
        if abs(pivot) < ptol:
            print("遇到零主元，程序退出。")
            os._exit(os.EX_OK)
        for i in range(k+1, row):         
            multiply = A[i, k] / pivot
            A[i, k:] = A[i, k:] - multiply * A[k, k:] 
            b[i] = b[i] - multiply * b[k]                         
        if verbose:
            print(f"第{k}列消元完成，主元值为 {pivot:.6f}")
            print(np.concatenate((A, b.reshape(-1, 1)), axis=1))          

    x = np.copy(b) 
    x[column-1] = x[column-1] / A[column-1, column-1]  
    if verbose:
        print("回代求解阶段")
        print(f"x{column-1} = {x[column-1]:.6f}")      
    for k in range(column-2, -1, -1):                  
        x[k] = (x[k] - np.dot(A[k, k+1:], x[k+1:])) / A[k, k]  
        if verbose:
            dot_val = np.dot(A[k, k+1:], x[k+1:])
            print(f"x{k} = ({b[k]:.6f} - {dot_val:.6f}) / {A[k, k]:.6f} = {x[k]:.6f}")
    return x

# 实验数据
T = np.array([10.0, 29.5, 50.0, 69.5, 90.0])
V = np.array([1.032, 1.094, 1.156, 1.215, 1.273])

# 转化方程：V - 1 = aT + bT²，令y = V - 1
y = V - 1

# 构造正规方程组：将超定问题转化为可解的方阵问题
sum_T = np.sum(T)
sum_T2 = np.sum(T**2)
sum_T3 = np.sum(T**3)
sum_T4 = np.sum(T**4)
sum_yT = np.sum(y * T)
sum_yT2 = np.sum(y * T**2)

# 正规方程组的系数矩阵A和右侧向量B
A = np.array([[sum_T2, sum_T3],
              [sum_T3, sum_T4]])
B = np.array([sum_yT, sum_yT2])

# 用高斯消元法求解a和b
eps = sys.float_info.epsilon
a, b = GEshow(A, B, 5 * eps, 1)

print(f"\n求解得到的系数：a = {a:.6f}, b = {b:.6f}")

# 验证计算
T_new = 30.0
V_new = 1 + a * T_new + b * T_new**2
print(f"当温度为{T_new}℃时，体积预测值为 {V_new:.4f} cm³")

# 绘图展示拟合效果
plt.figure(figsize=(6, 6))
plt.scatter(T, V, color='blue', label='原始数据')
T_line = np.linspace(T.min(), T.max(), 500)
V_line = 1 + a * T_line + b * T_line**2
plt.plot(T_line, V_line, color='red', label='拟合曲线')
plt.xlabel('温度 (℃)')
plt.ylabel('体积 (cm³)')
plt.legend()
plt.grid(True)
plt.show()

修正说明 #

1. 方程转化：将原公式变形为V - 1 = aT + bT²，仅保留a、b两个待求参数，简化求解逻辑
2. 正规方程组构造：针对超定数据构造最小二乘对应的正规方程组，将问题转化为高斯消元可处理的方阵求解问题
3. 维度匹配修正：构造2×2的系数矩阵和长度为2的向量，匹配待求系数数量，确保高斯消元正常运行
4. 细节优化：修复原代码中向量维度错误，调整输出为中文提示，提升可读性

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