Метод наименьших квадратов аппроксимации функций

Прежде чем запускать основной код, нужно сформировать код для метода Гаусса и проинициализировать его.

Вспомогательный код для решения СЛАУ методом Гаусса (подробнее в этой статье):

import numpy

def Metod_Gaussa(arr, brr):
    for k in range(arr.shape[0] - 1):
        #поиск строки с максимальным элементом
        max_elem = 0
        str = 0
        for i in range (k, arr.shape[0]):
            if abs(arr[i,k]) > abs(max_elem):
                max_elem = arr[i,k]
                str = i
        #меняем местами строки квадратной матрицы
        change = numpy.repeat(arr[k], 1)
        arr[k], arr[str] = arr[str], change
        #меняем местами элементы вектора-столбца
        change = numpy.repeat(brr[k], 1)
        brr[k], brr[str] = brr[str], change
        #делим полученную строку на max_elem
        arr[k] = arr[k] / max_elem
        brr[k] = brr[k] / max_elem
        #домножаем строку на коэффициенты и вычитаем ее из остальных строк
        for i in range (k + 1, arr.shape[0]):
            factor = arr[i,k]
            arr[i] = arr[i] - arr[k] * factor
            brr[i] = brr[i] - brr[k] * factor

    #находим неизвестные
    arg = [brr[brr.shape[0] - 1] / (arr[arr.shape[0] - 1, arr.shape[0] - 1])]
    for i in range(arr.shape[0] - 2, -1, -1):
        n = brr[i]
        for j in range(len(arg)):
            n = n - arg[j] * arr[i, arr.shape[0] - 1 - j]
        arg.append(n)

    #переворачиваем значения в списке
    otv = []
    for i in reversed(arg): otv.append(i)
    return otv

Этот код следует добавить в отдельный файл в ту же папку PyCharm, что и основной код:

Основной код:

#библиотека для параметризации функции и построения графиков
from sympy import *

#для построения набора точек на графике
import matplotlib.pyplot as plt

#для решения СЛАУ
from Project1 import Gaussa
import numpy

x = symbols('x')
K = 3 #степень многочлена

#опорные точки
xi = [1, 1.25, 1.25, 1.5, 1.75, 2]
yi = [5.47863093565943e-837, 0.000244140625000000, 1, 0.0156250000000000, 0.177978515625000, 1]

def approksim(xi, yi, K):
    #поиск констант
    A = numpy.zeros([K + 1, K + 1])
    for i in range(K + 1):
        for j in range(K + 1):
            for t in xi: A[i, j] += t ** (2 * K - i - j)
    a = numpy.zeros([K + 1])
    for i in range(K + 1):
        for j in range(len(xi)): a[i] += yi[j] * xi[j] ** (K - i)
    const = []
    for j in Gaussa.Metod_Gaussa(A, a): const.append(float(j))
    #составление уравнения
    f = 0
    for i in range(len(const)):
        f += const[i] * x ** (len(const) - 1 - i)
    return f

#значения линейной функции, полученной МНК, в тех же точках
f1 = []
for i in xi: f1.append(approksim(xi, yi, 1).evalf(subs={'x': i}))
#отклонения (невязки в точках)
e1 = []
for i in range(len(xi)): e1.append(abs(f1[i] - yi[i]))

#значения квадратичной функции, полученной МНК, в тех же точках
f2 = []
for i in xi: f2.append(approksim(xi, yi, 2).evalf(subs={'x': i}))
#отклонения (невязки в точках)
e2 = []
for i in range(len(xi)): e2.append(abs(f2[i] - yi[i]))

#значения кубической функции, полученной МНК, в тех же точках
f3 = []
for i in xi: f3.append(approksim(xi, yi, 3).evalf(subs={'x': i}))
#отклонения (невязки в точках)
e3 = []
for i in range(len(xi)): e3.append(abs(f3[i] - yi[i]))

#все функции на одном графике
xfi = numpy.linspace(xi[0], xi[len(xi) - 1], 100)
f1i = [approksim(xi, yi, 1).subs(x, a) for a in xfi]
f2i = [approksim(xi, yi, 2).subs(x, a) for a in xfi]
f3i = [approksim(xi, yi, 3).subs(x, a) for a in xfi]
plt.plot(xfi, f1i, 'b')
plt.plot(xfi, f2i, 'r')
plt.plot(xfi, f3i, 'g')
plt.plot(xi, yi, 'k*')
plt.grid()
plt.show()

Результат (линейная, квадратичная и кубическая аппроксимации на одном графике):