pythondata-sciencepolynomialsapproximationpolynomial-approximations

Approximating the sine function with polynomials using gaussian elimination


I'm trying to approximate the sine function over an interval with gaussian elimination using python. Using this code.

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()

But whenever I try making the increment smaller (to make the approximation function more accurately supposedly) python keeps giving me this error.

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')

Is this just because I have too many points so x_points[i] ** j becomes too big? or have I made a mistake somewhere? and even if I do make it work by making the increment larger some of the points don't match up with the sin function. 0.1 increment plot. Test_func is the sine function and f is the approximation function.

Does anyone know why this happens?

Here is one more screenshot over an increment of 0.07 within the same interval as the one in the code. 0.07 Increment plot. If there are other things that might help with this please let me know.


Solution

  • In such cases, it helps to use a debugger to find out what's going wrong. Or at the very least, use try-except blocks to catch the exception and print your variables to see where your code explodes.

    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
    

    In this case, your code breaks when i = 310 and j = 628. This is because x_points[i] ** j = 3.10 ** 628, which is too big to store in a double.

    (Note that I specifically catch the OverflowError exception because the error message tells me that's the error that was raised. It is not a good idea to have an unfiltered except clause that catches all exceptions, because that can hide actual problems in your code and prevent other functionality such as exiting with Ctrl+C)