计算方法 | 曲线拟合的最小二乘法

就是找个破多项式，让每个点和真实数据的差值的平方最小，原理就是这么简单。

1.写正则方程组（法方程组）

2.写结果向量

3.解方程组

丢代码：

import numpy as np
import math
# data = eval(input("请输入数据："))
data = [(0.6,7.05),(1.3,12.2),(1.64,14.4),(1.8,15.2),(2.1,17.4),(2.3,19.6),(2.44,20.2)]
n = 2
m = len(data)
fa_mat = [[0 for j in range(n)] for i in range(n)]
for i in range(len(fa_mat)):
    for j in range(len(fa_mat[0])):
        _sum = 0
        for k in range(m):
            _sum += data[k][0] ** (i+j)
        fa_mat[i][j] = _sum
print("法方程组如下：")
print(np.array(fa_mat))
# 计算结果向量
b = [[0] for i in range(n)]
for i in range(n):
    _sum = 0
    for k in range(m):
        _sum += ((data[k][0] ** i) * data[k][1])
        # print("_sum += (({} ** {}) * {})".format(data[k][0],i,data[k][1]))
    b[i][0] = _sum
print("结果向量如下：")
print(np.array(b))
print("下面是解方程组：\n\n")
# 解线性方程组 使用三角分解法 ======================================
mat = fa_mat[:]
b = b[:]
height = len(mat)
width = len(mat[0])

# 1 直接分解矩阵
# 初始化LU
U = [[0 for i in range(width)] for j in range(height)]
L = [[0 if i!=j else 1 for i in range(width) ] for j in range(height)]

print(L,U)

# 从U第一行开始，对于L来说是从第一列开始
for i in range(height):
    # 依次计算LU的行、列
    for j in range(i, width):
        _sum = 0
        for k in range(i):
            _sum += L[i][k] * U[k][j]
            print("U{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(i,j,i,k,k,j,L[i][k],U[k][j],L[i][k] * U[k][j]))
        U[i][j] = mat[i][j] - _sum
    # L的i.j是反着的
    for j in range(i, width):
        _sum = 0
        if j > i:
            for k in range(i):
                _sum += L[j][k] * U[k][i]
                print("L{}{}_sum += L[{}][{}] * U[{}][{}] = {} * {} = {}".format(j,i,j,k,k,i,L[j][k],U[k][i],L[j][k] * U[k][i]))
            L[j][i] = (mat[j][i] - _sum) / U[i][i]
    print("L = \n{}".format(np.array(L)))
    print("U = \n{}".format(np.array(U)))
    
print("分解完毕，开始计算Ly=b")
# 分解完毕 解Ly =b
y = [[0] for i in range(height)]
# print(y)
for i in range(height):
    _sum = 0
    for k in range(i):
        _sum += L[i][k] * y[k][0]
    y[i][0] = (b[i][0] - _sum)          # 这里不用除以系数因为对角线为1
    print(y)
print("计算完毕，开始计算Ux=y")
x = [[0] for i in range(height)]
for i in range(height-1,-1,-1):
    _sum = 0
    for k in range(height-1,i,-1):
        _sum += U[i][k] * x[k][0]
    x[i][0] = (y[i][0] - _sum) / U[i][i]
    print(x)
print("最终结果x: \n{}".format(np.array(x)))

over.