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 res = linalg.solve(AA, bb) print(res)