Я самостоятельно изучаю некоторые числовые алгоритмы из книги 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
```