HOUSEHOLDER

import math import scipy.linalg as sla import numpy as np import matplotlib.pyplot as plt from math import copysign, hypot #hypot = euclidian norm from numpy import linalg %matplotlib inline N = 7 M = 4 def func(x): return x*math.sin(x) x = np.zeros(shape=(M), dtype=np.float128(0.)) y = np.zeros(shape=(M), dtype=np.float128(0.)) for i in range(0, M): x[i] = -1. + 2. / (M - 1) * i y[i] = func(x[i]) #Матрица степеней полиномов Лежандра G = np.zeros((max(2, N), max(2, N)), dtype=np.float128(0.)) G[0][0] = 1 # g_0 = 1*x^0 G[1][1] = 1 # g_1 = 1*x^1 for i in range(2, N): # np.roll - Roll array elements along a given axis. G[i] = (np.roll(G[i - 1], 1) * (2 * i - 1.) - (i - 1) * G[i - 2]) / i px = np.zeros((max(2, N), M), dtype=np.float128(0.)) for i in range(0, M): px[0][i] = 1. for i in range(1, N): for j in range(0, M): px[i][j] = px[i - 1][j] * x[j] #print(px) G_px = np.dot(G,px) #print(G_px) def householder_reflection(A, b): (num_rows, num_cols) = np.shape(A) R = np.copy(A) for cnt in range(num_rows - 1): x = R[cnt:, cnt] e = np.zeros_like(x) e[0] = copysign(np.linalg.norm(x), -A[cnt, cnt]) u = x + e v = u / np.linalg.norm(u) Q_cnt = np.identity(num_rows) Q_cnt[cnt:, cnt:] -= 2.0 * np.outer(v, v) R = np.dot(Q_cnt, R) b = np.dot(Q_cnt, b) return (R, b) A = G_px.dot(G_px.T) # A = ggT - Gram's matrix b = G_px.dot(y.T) # b = gyT - also mult the b vector (AA, bb) = householder_reflection(A, b) # QR decomp of original system ans = sla.solve_triangular(AA, bb) #вычисление функци Psi for i in range(0, N): G[i] *= ans[i] #вычисление Psi(x) для данного x def G_at_x(x): w = np.zeros(shape=(N), dtype=np.float128(0.)) w[0] = 1. for i in range(1, N): w[i] = w[i - 1] * x #строка степеней x return np.sum(np.dot(G,w)) #значение Psi в точке х eps = 1e-4 drawx = np.arange(x[0], x[M - 1], eps) #построение графика def get_graph(xx, f): drawy = np.array(xx) for i in range(0, drawy.shape[0]): drawy[i] = f(drawx[i]) return drawy plt.figure(num=None, figsize=(12,10), dpi=200, facecolor='w', edgecolor='k') #отрисовка Psi drawy = get_graph(drawx, G_at_x) curve1, = plt.plot(drawx, drawy, 'm-', linewidth=1.) #отрисовка исходной функции drawy = get_graph(drawx, func) curve2, = plt.plot(drawx, drawy, 'g-', linewidth=1.) #узлы curve3, = plt.plot(x, y, 'ko') plt.legend(("Psi", "Function"),loc='upper right') plt.setp(curve1, linewidth=3) plt.setp(curve2, linewidth=3) plt.title('Approximation Results') plt.show()