1

我正在尝试使用 python 用高斯消除来近似正弦函数。使用此代码。

from copy import deepcopy
def convert_to_row_eschelon(A_, B_):
    A = deepcopy(A_)
    B = deepcopy(B_)
    dim = len(A)
    for cc in range(dim):
        # pivot_row = A[cc]
        for r in range(cc + 1, dim):
            leading_term = A[r][cc]
            for c in range(cc, dim):
                # print(A[r][c], A[r][cc])
                A[r][c] = A[r][c] - A[cc][c] * leading_term / A[cc][cc]
            B[r] = B[r] - B[cc] * leading_term / A[cc][cc]
    return A, B

def back_sub(matrix_pair):
    A = matrix_pair[0]
    B = matrix_pair[1]
    res = [None] * len(B)
    
    for i in range(len(B) - 1, -1, -1):
        def f(j):
            return A[i][j] * res[j]
        res[i] = (B[i] - sum([f(k) for k in range(i + 1, len(B))])) / A[i][i]
    return res


def gaussian_elimination(A, B):
    return back_sub(convert_to_row_eschelon(A, B))
A = [
      [1, 2, 3],
      [4, 5, 7],
      [23, 12, 12]
]
B = [4, 6, 7]
fig = 10
# print(convert_to_row_eschelon(A, B))
def make_polynomial(x_points, y_points):
    # A[x_point index used][degree]
    degree = len(x_points)
    A = []
    for i in range(degree):
        A.append([])
        for j in range(degree):
            A[i].append(x_points[i] ** j) # This is line 45 
    coeff = gaussian_elimination(A, y_points)
    def f(x):
        coeff_f = coeff
        res = 0
        for i in range(len(coeff_f)):
            res += x ** i * coeff_f[i]
        return res
    return f

def generate_x(start, finish, increment):
    x_points = []
    curr = start
    while curr < finish:
        x_points.append(curr)
        curr += increment
    return x_points
from math import sin, pi
start = 0 # These are the intervals
finish = 2 * pi
increment = 0.01
def test_func(x):
    return  sin(x)
# Creating the polynomial
x_val_f = generate_x(start, finish, increment)
x_val_test = generate_x(start, finish, 0.01)

f = make_polynomial(x_val_f, [test_func(i) for i in x_val_f])
print(f(3))
y_val_f = [f(i) for i in x_val_f]
y_val_test = [test_func(i) for i in x_val_test]
error = sum([abs(y_val_f[i] - y_val_test[i]) for i in range(len(y_val_f))]) / len(y_val_f)
print('average error : {}'.format(error))
# plotting f
import matplotlib.pyplot as plt


plt.plot(x_val_test, y_val_test, label = "test_func")
plt.scatter(x_val_f, y_val_f, label = "f(x)", s = 10)
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.ylim(-1,1) 
plt.title('Graph')
plt.legend()
plt.show()

但是，每当我尝试使增量更小（据说使近似函数更​​准确）时，python 都会不断给我这个错误。

File "c:/Users/username/Desktop/Curve fitting.py", line 45, in make_polynomial
A[i].append(x_points[i] ** j)
        OverflowError: (34, 'Result too large')

难道这只是因为我有太多的积分所以x_points[i] ** j变得太大了吗？还是我在某个地方犯了错误？即使我确实通过使增量更大来使其工作，但某些点与 sin 函数不匹配。 0.1 增量图。Test_func 是正弦函数，f 是近似函数。

有谁知道为什么会这样？

这是在与代码中的相同间隔内增加 0.07 的另一个屏幕截图。0.07 增量图。如果还有其他可能对此有所帮助的事情，请告诉我。

4

1 回答 1

2

在这种情况下，使用调试器找出问题所在是有帮助的。或者至少，使用try-except块来捕获异常并打印您的变量以查看您的代码在哪里爆炸。

try:
    for i in range(degree):
        A.append([])
        for j in range(degree):
            A[i].append(x_points[i] ** j) # This is line 45 
except:
    print(f"i = {i}; j = {j}")

### output: 
i = 310; j = 628

i = 310在这种情况下，您的代码会在和时中断j = 628。这是因为x_points[i] ** j = 3.10 ** 628, 它太大而无法存储在双精度中。

于 2020-07-07T15:12:51.393 回答