我怎样才能找到系数？

Question

V (cm3) T (◦C)
1.032 10.0
1.094 29.5
1.156 50.0
1.215 69.5
1.273 90.0

编写一个Python程序，通过假设求出系数a和b 该液体体积随关系式 V = 1 + aT + bT 2 变化。 我必须使用高斯消去法。

我不是CS专业的，所以对我来说真的很难。

有一个示例代码，但我未能将其转化为这个方程。

这是我来的最远的地方

import numpy as np
import matplotlib.pyplot as plt
import os,sys
def GEshow(A,b,ptol, verbose):

    peps = max(ptol, 5 * eps)                       
    print("Tolerance value in pivoting is %e." % (peps))
    row,column=np.shape(A)
    if row  != column:
        print('A matrix needs to be square! Exiting.')
        os._exit(os.EX_OK)
    if verbose:
        print('Gaussian Elimination with Augmented system:')
        print(np.concatenate((A, b), axis=1))              

    for k in range(column):                                  
        print("--- Stage %d" % k)
        pivot = A[k,k]
        if abs(pivot)<ptol:
            print("Zero Pivot Encountered. Exiting.")
            os._exit(os.EX_OK)
        for i in range(k+1,column):        
            multiply=A[i,k]/pivot
            A[i,k:column] = A[i,k:column] - (multiply)*A[k,k:column] 
            b[i] = b[i] -  (multiply) * b[k]                         
        if verbose:
            print('After elimination in column %d with pivot = %f' %(k,pivot))
            print(np.concatenate((A, b), axis=1))          

    x=np.copy(b) 
    x[column-1] = x[column-1]/A[column-1,column-1]  
    if verbose:
        print("Backward substitution")
        print("x%d=%f" % (column-1,x[column-1]))     
    for k in range(column-2,-1,-1):                  
        x[k] = (x[k] - np.dot(A[k,k+1:column],x[k+1:column]))/A[k,k]  
        if verbose:
            print("x%d=(%f-%f)/%f)=%f" % (k,b[k],np.dot(A[k,k+1:column], x[k+1:column]),A[k,k],x[k]))
    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])

A = np.zeros((len(T), 5))
A[:, 0] = 1
A[:, 1] = T
A[:, 2] = T**2
A[:, 3] = 0 
A[:, 4] = 0 

A = np.append(A, [[1, 1, 1]], axis=0)

eps = sys.float_info.epsilon
mycoefficients = GEshow(A, np.append(V, [0]), 5 * eps, 1)

mycoefficients = mycoefficients[:-1]

print(f"The coefficients are {mycoefficients}")

T_new = 30.0
V_new = mycoefficients[0] + mycoefficients[1]*T_new + mycoefficients[2]*T_new**2
print(f"At {T_new} degrees Celsius, the volume is {V_new} cm3")

plt.figure(figsize=(6, 6))

plt.scatter(T, V, color='blue', label='Original data')

T_line = np.linspace(T.min(), T.max(), 500)
V_line = mycoefficients[0] + mycoefficients[1]*T_line + mycoefficients[2]*T_line**2
plt.plot(T_line, V_line, color='red', label='Fitted polynomial')

plt.xlabel('Temperature (°C)')
plt.ylabel('Volume (cm3)')
plt.legend()
plt.grid(True)
plt.show()

Answer 1

在这里，我使用简单代数来计算直线，并使用二次回归来计算曲线。请注意直线如何精确拟合，而曲线则不然：

import numpy as np
import matplotlib.pyplot as plt

data = np.array([
    [1.032, 10.0],
    [1.094, 29.5],
    [1.156, 50.0],
    [1.215, 69.5],
    [1.273, 90.0]
])

def find_line(x,y):
    m = (y[-1]-y[0]) / (x[-1]-x[0])
    b = y[0] - m*x[0]
    return m,b

def find_quadratic(x,y):
    x2 = x*x
    xbar = x.mean()
    x2bar = (x*x).mean()
    ybar = y.mean()

    sxx = (x-xbar).sum()**2
    sxy = ((x-xbar)*(y-ybar)).sum()
    sxx2 = ((x-xbar)*(x2-x2bar)).sum()
    sx2x2 = ((x2-x2bar)**2).sum()
    sx2y = ((x2-x2bar)*(y-ybar)).sum()
    print(sxx,sxy,sxx2,sx2x2,sx2y)

    b = (sxy*sx2x2 - sx2y*sxx2) / (sxx*sx2x2 - sxx2*sxx2)
    c = (sx2y*sxx - sxy*sxx2) / (sxx*sx2x2 - sxx2*sxx2)
    a = ybar - b * xbar - c * x2bar
    return a,b,c

x = data[:,1]
y = data[:,0]

m,b = find_line(x,y)
mycoefficients = find_quadratic(x,y)

plt.figure(figsize=(6, 6))

plt.scatter(x, y, color='blue', label='Original data')

T_line = np.linspace(x.min(), x.max(), 500)
V_curve = mycoefficients[0] + mycoefficients[1]*T_line + mycoefficients[2]*T_line**2
V_line = m*T_line + b
plt.plot(T_line, V_curve, color='red', label='Fitted polynomial')
plt.plot(T_line, V_line, color='green', label='Fitted line')

plt.xlabel('Temperature (°C)')
plt.ylabel('Volume (cm3)')
plt.legend()
plt.grid(True)
plt.show()

输出：