Реализация базовой полиномиальной интерполяции [closed]

Я самостоятельно изучаю некоторые числовые алгоритмы из книги Numerical methods in Scientific computing пользователя Dahlquist / Bjorck. В задаче 4.1.2 предлагается написать программу, которая находит вектор коэффициентов c для полиномиального интерполянта. Я думаю, что моя математика верна, но сгенерированный таким образом многочлен соответствует точкам $ {(1955,100), (1960,117.7), ldots } $, что ошибочно, но я не могу определить проблему. Я бы хотел, чтобы кто-нибудь просмотрел мой код / ​​алгоритм.

Оригинальный вопрос из учебника для полноты картины.

Вопрос. (а) Напишите программу $ c = polyapp (x, y, n) $ который находит вектор коэффициентов $ c $ для полинома от $ p in mathcal {P} _n $, в основе сдвинутой мощности, так что $ y_i приблизительно p (x_i) $ за $ i = 1: m, m ge n $в смысле наименьших квадратов, или изучите программу, которая делает это почти так.

(b) Приведенные ниже данные показывают динамику валового внутреннего продукта (ВВП) Швеции, взятую из таблицы, составленной группой, связанной с конфедерацией шведских работодателей. (Данные выражены в ценах 1985 г. и масштабированы таким образом, что значение для 1950 г. равно 100.

x = np.array([1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990])
y = np.array([100.0, 117.7,139.3,179.3,219.3,249.1,267.5,291.5,326.4])

Для верхних пар данных вычислить и построить $ p (x) $, $ x в [1950,2000] $. Отметьте указанные точки данных. Сделай это для $ m = 9 $ и скажем $ п = 9 $ а затем для $ n = 8: -2: 2 $. Сохраните результаты, чтобы потом можно было сравнить их.

Моя попытка.

Сдвинутая основа власти для $ mathcal {P} _n $ дан кем-то,

$$ begin {align *} p_1 (x) & = 1 \ p_2 (x) & = (x — d) \ p_3 (x) & = (x — d) ^ 2 \ vdots \ p_n (x) & = (x — d) ^ n end {align *} $$

У нас есть $ m> п $ точки данных. Проблема интерполяции: $ p (x_i) = f (x_i) $ за $ i = 1: m $. В матричной форме $ M_n ( mathbf {p}) mathbf {c} = f $.

Нам интересно найти вектор коэффициентов $ mathbf {c} $ такая, что сумма квадратов остатков

$$ S ( mathbf {c}) = sum_i left[sum_{k=1}^{n}c_k p_k(x_i) — f(x_i) right]^ 2 $$

сводится к минимуму. Частная производная по $ c_k $ является,

$$ begin {align *} frac { partial S ( mathbf {c})} { partial c_k} & = 2 sum_ {i = 1} ^ {m} p_k (x_i) left ( sum_ {j = 1} ^ {n} c_j p_j (x_i) -f (x_i) right) \ & = 2 sum_ {i = 1} ^ {m} p_k (x_i) left ( sum_ {j = 1} ^ {n} c_j p_j (x_i) right) -2 n sum_ {i = 1} ^ {m} p_k (x_i) f (x_i) end {align *} $$

Функция $ S $ минимизируется в критической точке, где $ mathbf {c} _0 $ где вектор градиента и, следовательно, все частные производные равны нулю.

import numpy.linalg as linalg

def polyapp(x,y,n):
    m = len(x) # No of data points
    a = min(x) # Determine the interval [a,b]
    b = max(x)
    d = (a + b)/2 
    
    # We assume a shifted power basis for a polynomial p in P_n.
    def p(x,k):
        return (x - d)**k
    
    # We find the coefficient vector c for a polynomial p in P_n, such that 
    # y_i = p(x_i) in a least squares sense. The coefficient vector is given by the solution of a square linear
    # algebraic system.
    
    A = np.zeros((n,n))
    
    # The term multiplying c_j in the k'th equation is sum_{i=1}^m p_k(x_i) p_j(x_i)
    
    for k in range(0,n):
        for j in range(0,n):
            for i in range(0,m):
                A[k,j] = A[k,j] + p(x[i],k) * p(x[i],j)
                
    
    # Construct the right hand side vector 
    b = np.zeros(m)
    
    for k in range(0,n):
        for i in range(0,m):
            b[k] = b[k] + y[i]*p(x[i],k)
            
    c = linalg.solve(A,b)
    return A,c

def polyval(x,c,d,n):
    y = 0
    for j in range(0,n):
        y = y + c[j]*((x - d)**j)
    return y
```