In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
5.1.14 Функция $y=f(x)$ задана таблицей значений $y_0, y_1, \dots, y_n$ в точках $x_0, x_1, \dots, x_n$. Используя метод наименьших квадратов (МНК), найти многочлен $P_m(x)=a_0+a_1x+...+a_mx^m$ наилучшего среднеквадратичного приближения оптимальной степени $m=m^*$¶
In [2]:
x_1 = np.arange(-3, 2.5, 0.45)
y_1 = np.array([0.262, -1.032, -1.747, -1.981, -0.564, 0.774, 2.400, 2.131, 2.2, -0.393, -1.815, -0.788, 8.030])
In [15]:
np.sum(y_1*x_1**2)
Out[15]:
In [26]:
mnk(x_1, y_1, x_1.shape[0], 2)
Out[26]:
In [24]:
-3.9*0.792+38.025*0.723
Out[24]:
In [27]:
sigma(x_1, y_1, mnk(x_1, y_1, x_1.shape[0], 2), x_1.shape[0], 2)
Out[27]:
In [16]:
def mnk(x, y, n, m):
b = np.zeros(m+1)
A = np.zeros((m+1, m+1))
for j in range(m + 1):
b[j] = np.sum(y * (x ** j))
for k in range(m + 1):
A[j, k] = np.sum(x ** (k + j))
a = np.linalg.solve(A, b)
return a
def P(a, m, t):
p = 0
for j in range(m+1):
p += a[j] * (t ** j)
return p
def sigma(x, y, a, n, m):
d = 0
for k in range(n):
d += (P(a, m, x[k]) - y[k]) ** 2
d = np.sqrt(d / (n - m))
return d
In [29]:
sigmas = []
min_e = np.inf
m = 0
cur_e = sigma(x_1, y_1, mnk(x_1, y_1, x_1.shape[0], m), x_1.shape[0], m)
sigmas.append(cur_e)
while (cur_e < min_e) or (m < 10):
min_e = cur_e
m += 1
cur_e = sigma(x_1, y_1, mnk(x_1, y_1, x_1.shape[0], m), x_1.shape[0], m)
sigmas.append(cur_e)
print("σ = ", sigmas)
In [30]:
plt.plot(range(len(sigmas)), sigmas, '.-')
plt.ylabel('$\sigma_m$', size=18)
plt.xlabel('$m$', size=18)
plt.show()
In [32]:
m_s = np.argmin(sigmas)
print("m* = ", m_s)
In [56]:
plt.figure(figsize=(12, 10))
plt.plot(x_1, y_1, 'or', label="$y=f(x)$")
for m in range(m_s +1):
plt.plot(x_1, P(mnk(x_1, y_1, x_1.shape[0], m), m, x_1), label="$y=P_%d(x)$"%(m))
plt.legend(fontsize=14)
plt.grid(True)
plt.show()
In [46]:
b_t = np.zeros(5+1)
A_t = np.zeros((5+1, 5+1))
for j in range(5 + 1):
b_t[j] = np.sum(y_1 * (x_1 ** j))
for k in range(5 + 1):
A_t[j, k] = np.sum(x_1 ** (k + j))
print("Определитель матрицы: ", np.linalg.det(A_t))
print("Матрица не сингулярная.")
In [47]:
A_t
Out[47]:
In [48]:
Q, R = np.linalg.qr(A_t)
a_t = np.dot(np.dot(np.linalg.inv(R), Q.T), b_t)
print("a = ", a_t)
In [49]:
print("Погрешность при использовании QR-разложения матрицы по сравнению с обычным методом: ", np.max(np.abs(a_t - mnk(x_1, y_1, x_1.shape[0], 5))))
5.3.7 Зависимость между величинами $x$ и $y$ описывается функцией $y=f(x, a, b)$, где $a$ и $b$ – неизвестные параметры. Найти эти параметры, сведя исходную задачу к линейной задаче метода наименьших квадратов.¶
In [50]:
x_2 = np.arange(-4, 4.1, 0.8)
y_2 = np.array([-6.47, -3.2086, -2.3433, -2.2767, -1.4114, 1.85, 9.105, 21.951, 41.986, 70.806, 110.01])
t_2 = (x_2 + 2) ** 3
In [55]:
for i in t_2:
print(i)
In [12]:
temp_2 = mnk(t_2, y_2, y_2.shape[0], 1)
print("a = ", temp_2[0], ", b = ", temp_2[1])
5.7.7 Дана кусочно-гладкая функция $y=f(x)$. Сравнить качество приближения функции кусочно-линейной и глобальной интерполяциями¶
In [2]:
a_7 = -5
b_7 = 5
def func_7(x):
return x * (np.abs(x) - 4)
x_7 = np.linspace(a_7, b_7, num=51)
y_7 = func_7(x_7)
print("x = ", x_7)
print("\ny = ", y_7)
In [3]:
def lin_interp(x, y):
int_x = np.array([])
int_y = np.array([])
for i in range(x.shape[0] - 1):
b = (y[i] - y[i+1]) / (x[i] - x[i+1])
a = y[i] - b*x[i]
t_x = 0.2 * np.random.sample(3) + x[i]
int_x = np.append(int_x, t_x)
int_y = np.append(int_y, a + b * t_x)
return int_y, int_x
In [4]:
intp_y, intp_x = lin_interp(x_7, y_7)
In [5]:
def new_interp(x, y, t):
n = x.shape[0]
f = np.zeros((n, n))
f[:, 0] = y
for k in range(1, n):
for i in range(n-k):
f[i, k] = (f[i+1, k-1] - f[i, k-1]) / (x[i+k] - x[i])
s = y[0]
r = 1
for k in range(0, n - 1):
r = r * (t - x[k])
s += f[0, k + 1] * r
return s
In [16]:
plt.figure(figsize=(12, 10))
plt.plot(x_7[1:-1], y_7[1:-1], 'or', label="$y=f(x)$")
plt.plot(intp_x, intp_y, label="$Линейная \ интерполяция$")
intp_x_s = np.sort(intp_x)
plt.plot(intp_x_s[30:-30], new_interp(x_7, y_7, intp_x_s)[30:-30], label="$Интерполяционный \ многочлен \ в \ форме \ Ньютона$")
plt.legend(fontsize=14)
plt.grid(True)
plt.show()
In [32]:
d_lin_int = np.abs(intp_y - func_7(intp_x))
d_new_int = np.abs(new_interp(x_7, y_7, intp_x_s) - func_7(intp_x))
fig, ax = plt.subplots(2, 1, figsize = (12, 20))
ax[0].plot(intp_x, d_lin_int, '.r', label="$Линейная \ интерполяция$")
ax[0].plot(intp_x_s, d_new_int, '.', label="$Интерполяционный \ многочлен \ в \ форме \ Ньютона$")
ax[0].legend(fontsize=14, loc=2)
ax[0].set_xlabel("$x$", size=16)
ax[0].set_ylabel("$\Delta$", size=16)
ax[1].plot(intp_x, d_lin_int, '.r', label="$Линейная \ интерполяция$")
ax[1].plot(intp_x_s[30:-30], d_new_int[30:-30], '.', label="$Интерполяционный \ многочлен \ в \ форме \ Ньютона$")
ax[1].legend(fontsize=14, loc=2)
ax[1].set_xlabel("$x$", size=16)
ax[1].set_ylabel("$\Delta$", size=16)
ax[0].grid(True)
ax[1].grid(True)
fig.show()
Как видно по графикам, линейная интерполяция даёт лучшее приближение для заданной функции. Глобальная интерполяция в некоторых точках даёт большую погрешность, но в некоторых точках приближение имеет погрешность меньше, чем при линейной интерполяции. Такая погрешность в глобальной интерполяции объясняется неправильным выбором узлов интерполяции для данной функции.¶
Вывод: глобальную интерполяцию выгодно использовать при определенных видах функций, к тому же, нужно правильно выбирать узлы интерполяции. При хорошем подборе параметров глобальная интерполяция даёт лучшее приближение функции на точках, не лежащих в исходном отрезке интерполяции, чем линейная.¶
5.9.7 Дана функция $y=f(x)$. Приблизить $f(x)$ на отрезке $[a, b]$ методом глобальной интерполяции и указанным в индивидуальном варианте сплайном (с отсутствием узла). На одном чертеже построить графики приближающей функции и функции $f(x)$. Сравнить качество приближения при разном количестве узлов интерполяции.¶
In [17]:
def interp_spline_nak(x, y, t):
n = x.shape[0] - 1
h = np.zeros(n)
A = np.zeros((n + 1, n + 1))
B = np.zeros(n + 1)
s = np.zeros(n + 1)
S = np.zeros(n)
for i in range(n):
h[i] = x[i + 1] - x[i]
for i in range(1, n):
A[i, i - 1] = h[i - 1] ** (-1)
A[i, i] = 2 * (h[i - 1] ** (-1) + h[i] ** (-1))
A[i, i + 1] = h[i] ** (-1)
B[i] = 3 * (h[i - 1] ** (-2) * (y[i] - y[i - 1]) + h[i] ** (-2) * (y[i + 1] - y[i]))
A[0, 0] = h[0] ** (-2)
A[0, 1] = h[0] ** (-2) - h[1] ** (-2)
A[0, 2] = -h[1] ** (-2)
B[0] = - 2 * h[0] ** (-3) * (y[0] - y[1]) + 2 * h[1] ** (-3) * (y[1] - y[2])
A[n, n - 2] = h[n - 2] ** (-2)
A[n, n - 1] = h[n - 2] ** (-2) - h[n - 1] ** (-2)
A[n, n] = -h[n - 1] ** (-2)
B[n] = - 2 * h[n - 2] ** (-3) * (y[n - 2] - y[n - 1]) + 2 * h[n - 1] ** (-3) * (y[n - 1] - y[n])
s = np.linalg.solve(A, B)
for i in range(1, n + 1):
t_a = ((t - x[i]) ** 2) * (2 * (t - x[i - 1]) + h[i - 1]) * y[i - 1] * (h[i - 1] ** (-3))
t_b = ((t - x[i - 1]) ** 2) * (2 * (x[i] - t) + h[i - 1]) * y[i] * (h[i - 1] ** (-3))
t_c = ((t - x[i]) ** 2) * (t - x[i - 1]) * s[i - 1] * (h[i - 1] ** (-2))
t_d = ((t - x[i - 1]) ** 2) * (t - x[i]) * s[i] * (h[i - 1] ** (-2))
S[i - 1] = t_a + t_b + t_c + t_d
return S
In [18]:
def func_9(x):
return 8 * np.exp(x) * np.cos(x ** 2)
In [22]:
def my_plot_func(x, y, t, plot=True):
plt_y = []
for i in t:
temp_y = interp_spline_nak(x, y, i)
t_pos = np.argwhere((i < x_9) == True)
if t_pos.shape[0]:
plt_y.append(temp_y[t_pos[0, 0] - 1])
else:
plt_y.append(temp_y[-1])
if plot:
plt.figure(figsize=(12, 10))
plt.plot(x, y, '*r', MarkerSize=14, label="$Узлы \ интерполяции$")
plt.plot(t, plt_y, label="$S_3(x)$")
plt.plot(t, func_9(t), label="$f(x)$")
plt.legend(fontsize=14)
plt.grid(True)
plt.show()
else:
return plt_y
In [23]:
a_9 = 1
b_9 = 3.75
x_9 = np.linspace(a_9, b_9, 5)
y_9 = func_9(x_9)
t_a = np.arange(a_9, b_9 + 0.1, 0.01)
my_plot_func(x_9, y_9, t_a)
In [24]:
x_9 = np.linspace(a_9, b_9, 20)
y_9 = func_9(x_9)
t_a = np.arange(a_9, b_9 + 0.1, 0.01)
my_plot_func(x_9, y_9, t_a)
При 5 узлах интерполяции построенная функция плохо приближает настоящую, в то время как 20-ти узлов хватает для хорошего приближения.¶
In [26]:
knots = range(4, 26)
deltas = []
for j in knots:
x_9 = np.linspace(a_9, b_9, j)
y_9 = func_9(x_9)
t_a = np.arange(a_9, b_9 + 0.1, 0.01)
intp_y_9 = my_plot_func(x_9, y_9, t_a, plot=False)
deltas.append(np.linalg.norm(intp_y_9 - func_9(t_a)))
In [30]:
plt.figure(figsize = (12, 10))
plt.plot(knots, deltas, 'o-')
plt.xlabel("$Количество \ узлов$", size=15)
plt.ylabel("$Погрешность$", size=15)
plt.grid(True)
plt.show()
