使用高斯消元法用多项式逼近正弦函数

Question

我正在尝试使用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 函数不匹配。 。 test_func 是正弦函数，f 是近似函数。

有谁知道为什么会发生这种情况？

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

Answer 1

在这种情况下，使用“调试器”有助于找出问题所在。或者至少，使用 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 OverflowError:
    print(f"i = {i}; j = {j}")

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

在这种情况下，您的代码会在 

i = 310

和

j = 628

时中断。这是因为

x_points[i] ** j = 3.10 ** 628

太大，无法存储在双精度中。（请注意，我专门捕获了

OverflowError

异常，因为错误消息告诉我这是引发的错误。使用未过滤的

except

子句捕获

all

异常并不是一个好主意，因为这可以隐藏代码中的实际问题并阻止其他功能，例如使用 Ctrl+C 退出）