ЛР6 по ЧМ (1) (1)

Лабораторная работа №6.

Задание 1.

Код:

import numpy as np from random import randint def print_sys(A, f): for i in range(len(A)): for j in range(len(A[i])): if j == 0 and A[i][j] > 0: print(' ', end='') if A[i][j] >= 0 and j != 0: print('+{:0.3f}'.format(A[i][j].round(3))+chr(j + 97), end='') else: print('{:0.3f}'.format(A[i][j].round(3))+chr(j + 97), end='') print('=', end='') if f[i] >= 0: print(' ', end='') print(f'{f[i].round(3)}') def print_sol(solution): for i in range(len(solution)): print(f'{chr(i + 97)} = {solution[i]}') def kramer(A:np.ndarray, f:np.ndarray)->np.ndarray: if not np.linalg.det(A): raise ValueError('Unable to solve due to matrix singularity') solution=np.array(f,np.float64) delta=np.linalg.det(A) for col in range(len(f)): Ai=A.copy() Ai[:,col]=f solution[col]=np.linalg.det(Ai)/delta return solution if __name__ == '__main__': A = np.array([[randint(-9, 9) for _ in range(3)], [randint(-9, 9) for _ in range(3)], [randint(-9, 9) for _ in range(3)]]) f = np.array([randint(-9, 9) for _ in range(3)]) f = f.transpose() print_sys(A, f) print('-' * 24) print_sol(kramer(A, f))

Вывод:

import numpy as np from lab6_1 import print_sys, print_sol from random import randint def raise_max_row(matrix, current_diag_index): ind = current_diag_index + np.argmax( np.abs(matrix[current_diag_index:, current_diag_index])) # seatching for row with max element if ind != current_diag_index: matrix[current_diag_index, :], matrix[ind, :] = np.copy(matrix[ind, :]), np.copy( matrix[current_diag_index, :]) # swap rows if needed def s_gauss(matrix): n = matrix.shape[0] # transforming matrix to triangular form for column in range(n - 1): raise_max_row(matrix, column) for row in range(column + 1, n): # modify row frac = matrix[row, column] / matrix[column, column] matrix[row, :] -= matrix[column, :] * frac # check modified system for nonsingularity if np.any(np.diag(matrix) == 0): return None solution = [0.0 for _ in range(n)] for column in range(n - 1, -1, -1): solution[column] = matrix[column, -1] / matrix[column, column] matrix[:, -1] -= matrix[:, column] * solution[column] matrix[:, column] -= matrix[:, column] return solution def gauss(A, f): Af = np.c_[A, f] Af = np.array(Af, np.float64) solution = s_gauss(Af) if solution is None: print("System has infinite solutions") return solution if __name__ == '__main__': A = np.array([[randint(-9, 9) for _ in range(3)], [randint(-9, 9) for _ in range(3)], [randint(-9, 9) for _ in range(3)]]) print('Исходная матрица\n', A) f = np.array([randint(-9, 9) for _ in range(3)]) print(f' {f}') f = f.transpose() print_sys(A, f) print_sol(gauss(A, f))

Задание 3.

from lab6_1 import kramer from lab6_2 import gauss import numpy as np if __name__ == '__main__': A = np.random.rand(20, 20) X = np.random.rand(20, 1).T[0] f = A.dot(X) while np.linalg.cond(A) < 1000: A = np.random.rand(20, 20) print(A) print(f' {f}') print(f'cond(A)={np.linalg.cond(A)}') kramer_solution = kramer(A, f) gauss_solution = gauss(A, f) simpleAF_solution=np.linalg.inv(A).dot(f) print('Premade: \n', X) print('Gauss: \n', gauss_solution) print('Gauss difference: \n', tuple(abs(gauss_solution[i] - X[i]) for i in range(len(gauss_solution)))) print('Kramer: \n', kramer_solution) print('Kramer difference: \n', tuple(abs(kramer_solution[i] - X[i]) for i in range(len(kramer_solution)))) print('A^(-1)*f: \n', simpleAF_solution) print('A^(-1)*f difference\n ', tuple(abs(simpleAF_solution[i] - X[i]) for i in range(len(simpleAF_solution)))) print('Max:\n') print(f'Отличие по Крамеру: {max(abs(kramer_solution[i] - X[i]) for i in range(len(kramer_solution)))}') print(f'Отличие по гаусу: {max(abs(gauss_solution[i] - X[i]) for i in range(len(gauss_solution)))}') print(f'A^(-1)*f отличие: {max(abs(gauss_solution[i] - X[i]) for i in range(len(gauss_solution)))}')

import numpy as np from random import randint noo = 0 def raise_max_row(matrix, current_diag_index): ind = current_diag_index + np.argmax( np.abs(matrix[current_diag_index:, current_diag_index])) # searching for row with max element if ind != current_diag_index: matrix[current_diag_index, :], matrix[ind, :] = np.copy(matrix[ind, :]), np.copy( matrix[current_diag_index, :]) # swap rows if needed def solve_gauss(matrix): global noo n = matrix.shape[0] # transforming matrix to triangular form for column in range(n - 1): raise_max_row(matrix, column) for row in range(column + 1, n): # modify row frac = matrix[row, column] / matrix[column, column] noo += 1 matrix[row, :] -= matrix[column, :] * frac noo += n # check modified system for nonsingularity if np.any(np.diag(matrix) == 0): return None solution = [0.0] * n for column in range(n - 1, -1, -1): solution[column] = matrix[column, -1] / matrix[column, column] noo += 1 matrix[:, -1] -= matrix[:, column] * solution[column] noo += 2 * n matrix[:, column] = np.zeros((1, n)) return solution def gauss(A, f): Af = np.c_[A, f] Af = np.array(Af, np.float64) solution = solve_gauss(Af) if solution is None: print("System has infinite solutions") return solution if __name__ == '__main__': noos = [] for i in range(3, 10): noo = 0 matr = [] for j in range(i): matr.append([randint(-9, 9) for _ in range(i)]) A = np.array(matr) f = np.array([randint(-9, 9) for _ in range(i)]) f = f.transpose() gauss(A, f) noos.append(noo) print('-' * 24) import matplotlib.pyplot as plt from matplotlib.widgets import Slider def update(): global a, c ax.cla() x = tuple((i for i in range(3, 10))) ax.plot(x, noos, 'ro', label='Real') x_i = list(np.linspace(3, 10, 1000)) ax.plot(x_i, tuple(i ** 3 * a + c for i in x_i), 'b', label="x^3") ax.plot(x_i, tuple(i ** 2 * a + c for i in x_i), 'b', label="x^2") ax.axvline(0, color='black') ax.axhline(0, color='black') ax.grid() ax.legend() fig, ax = plt.subplots() fig.subplots_adjust(bottom=0.3) x = tuple((i for i in range(3, 10))) fx = np.linspace(3, 10, 1000) ax.plot(fx, tuple(i ** 3 for i in fx), 'b', label="x^3") ax.plot(fx, tuple(i ** 2 for i in fx), 'b', label="x^2") ax.plot(x, noos, 'ro', label='Real') ax.axhline(color='k') ax.axvline(color='k') ax.legend() ax.grid() a = 2 / 3 nodes_slider_ax = fig.add_axes([0.25, 0.05, 0.65, 0.03]) nodes_slider = Slider(nodes_slider_ax, 'A', 0, 10, valinit=2 / 3) def temp(a_): global a a = a_ update() nodes_slider.on_changed(temp) nodes_slider_ax2 = fig.add_axes([0.25, 0.01, 0.65, 0.03]) nodes_slider2 = Slider(nodes_slider_ax2, 'B', 0, 100, valinit=20) c = 20 update() def temp2(c_): global c c = c_ update() nodes_slider2.on_changed(temp2) plt.show() # print_system(A, f) # print('-' * 24) # print_solution(gauss(A,f))