отчет. Отчет По Лабораторной Работе 8 ПавлюкИ. УдоратинА.. Отчет по лабораторной работе 8
 Скачать 0.55 Mb.
|
  Отчет по лабораторной работе №8 Подготовили студенты группы 112-МКо Павлюк Илья Олегович и Удоратин Андрей Иванович 2 8.1 Решить при своем n СЛАУ с матрицей Гильберта. Правая часть – вектор из единиц.      Листинг программы: myA=[ [1.0, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10], [1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11], [1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12], [1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13], [1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14], [1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15], [1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16], [1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17], [1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18], [1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19] ] myB = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] def FancyPrint(A, B, selected): for row in range(len(B)): print("(", end='') for col in range(len(A[row])): print("\t{1:10.2f}{0}".format(" " if (selected is None or selected != (row, col)) else "*", A[row][col]), end='') print("\t) * (\tX{0}) = (\t{1:10.2f})".format(row + 1, B[row])) def SwapRows(A, B, row1, row2): A[row1], A[row2] = A[row2], A[row1] B[row1], B[row2] = B[row2], B[row1] def DivideRow(A, B, row, divider): A[row] = [a / divider for a in A[row]] B[row] /= divider def CombineRows(A, B, row, source_row, weight): A[row] = [(a + k * weight) for a, k in zip(A[row], A[source_row])] B[row] += B[source_row] * weight def Gauss(A, B): column = 0 while (column < len(B)): print("Ищем максимальный по модулю элемент в {0}-м столбце:".format(column + 1)) current_row = None for r in range(column, len(A)): if current_row is None or abs(A[r][column]) > abs(A[current_row][column]): current_row = r if current_row is None: print("решений нет") return None FancyPrint(A, B, (current_row, column)) if current_row != column: print("Переставляем строку с найденным элементом повыше:") SwapRows(A, B, current_row, column) FancyPrint(A, B, (column, column)) print("Нормализуем строку с найденным элементом:") DivideRow(A, B, column, A[column][column]) FancyPrint(A, B, (column, column)) print("Обрабатываем нижележащие строки:") for r in range(column + 1, len(A)): CombineRows(A, B, r, column, -A[r][column]) FancyPrint(A, B, (column, column)) column += 1 print("Матрица приведена к треугольному виду ") X = [0 for b in B] for i in range(len(B) - 1, -1, -1): X[i] = B[i] - sum(x * a for x, a in zip(X[(i + 1):], A[i][(i + 1):])) print("Получили ответ:") print("\n".join("X{0} =\t{1:10.2f}".format(i + 1, x) for i, x in enumerate(X))) return X print("Исходная система:") FancyPrint(myA, myB, None) print("Получили ответ:") Gauss(myA, myB)    Листинг программы: import math import copy a = [ [1.0, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10], [1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11], [1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12], [1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13], [1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14], [1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15], [1/7, 1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16], [1/8, 1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17], [1/9, 1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18], [1/10, 1/11, 1/12, 1/13, 1/14, 1/15, 1/16, 1/17, 1/18, 1/19] ] b = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] def isCorrectArray(a): for row in range(0, len(a)): if( len(a[row]) != len(b) ): print('Не соответствует размерность') return False for row in range(0, len(a)): if( a[row][row] == 0 ): print('Нулевые элементы на главной диагонали') return False return True def isNeedToComplete(x_old, x_new): eps = 0.0001 sum_up = 0 sum_low = 0 for k in range(0, len(x_old)): sum_up += ( x_new[k] - x_old[k] ) ** 2 sum_low += ( x_new[k] ) ** 2 return math.sqrt( sum_up / sum_low ) < eps def solution(a, b): if( not isCorrectArray(a) ): print('Ошибка в исходных данных') else: count = len(b) x = [1 for k in range(0, count) ] numberOfIter = 0 MAX_ITER = 100 while( numberOfIter < MAX_ITER ): x_prev = copy.deepcopy(x) for k in range(0, count): S = 0 for j in range(0, count): if( j != k ): S = S + a[k][j] * x[j] x[k] = b[k]/a[k][k] - S / a[k][k] if isNeedToComplete(x_prev, x) : break numberOfIter += 1 print('Количество итераций : ', numberOfIter) return x print( 'Получили ответ: ', solution(a, b) ) 8.2. Найти минимум предложенной функции градиентным методом. Например, наискорейшего спуска или дроблением шага.   Листинг программы: import numpy as np from sympy import * import math import matplotlib.pyplot as plt import mpl_toolkits.axisartist as axisartist x1, x2, t = symbols('x1, x2, t') def func(): return pow(x1, 2) + 4 * pow(x2, 2) - 6 * x1 - 32 * x2 def grad(data): f = func() grad_vec = [diff (f, x1), diff (f, x2)] grad = [] for item in grad_vec: grad.append(item.subs(x1, data[0]).subs(x2, data[1])) return grad def grad_len(grad): vec_len = math.sqrt(pow(grad[0], 2) + pow(grad[1], 2)) return vec_len def zhudian(f): t_diff = diff(f) t_min = solve(t_diff) return t_min def main(X0, theta): f = func() grad_vec = grad(X0) grad_length = grad_len (grad_vec) k = 0 data_x = [0] data_y = [0] while grad_length> theta: k += 1 p = -np.array(grad_vec) X = np.array(X0) + t*p t_func = f.subs(x1, X[0]).subs(x2, X[1]) t_min = zhudian(t_func) X0 = np.array(X0) + t_min*p grad_vec = grad(X0) grad_length = grad_len(grad_vec) print('grad', grad_length) data_x.append(X0[0]) data_y.append(X0[1]) print(k) fig = plt.figure() ax = axisartist.Subplot(fig, 111) fig.add_axes(ax) ax.axis["bottom"].set_axisline_style("-|>", size=1.5) ax.axis["left"].set_axisline_style("->", size=1.5) ax.axis["top"].set_visible(False) ax.axis["right"].set_visible(False) plt.title(r'$Градиентный \ метод - метод \ наискорейшего \ спуска$') plt.plot(data_x, data_y, label=r'$f(x_1,x_2)=x_1^2+4 \cdot x_2^2-6 \cdot x_1 - 32 \cdot x_2$') plt.legend() plt.scatter(1, 1, marker=(5, 1), c=5, s=1000) plt.grid() plt.xlabel(r'$x_1$', fontsize=20) plt.ylabel(r'$x_2$', fontsize=20) plt.show() if __name__ == '__main__': main([0, 0], 0.00001) |