diff --git a/numeric_lab3/newton-jakobi.py b/numeric_lab3/newton-jakobi.py new file mode 100644 index 0000000..bfecf53 --- /dev/null +++ b/numeric_lab3/newton-jakobi.py @@ -0,0 +1,187 @@ +def f1(x1, x2): + return x2**2 + 0.8 * x2**2 + 0.1 - x1 + + +def f2(x1, x2): + return 2*x1*x2 + 0.1 - x2 + + +def newton_jakobi(): + eps = 0.001 + argument_vector = [1, 1] + prev_vector = [0, 0] + identity_matrix = initialize_matrix() + f = [0, 0] + while True: + f[0] = f1(argument_vector[0], argument_vector[1]) + f[1] = f2(argument_vector[0], argument_vector[1]) + jacobian = build_jacobian(argument_vector, f) + reversed_jacobian = reverse_matrix(jacobian, identity_matrix) + for i in range(len(f)): + x_old = argument_vector[i] + prev_vector[i] = x_old + increment = 0.0 + for j in range(2): + increment += reversed_jacobian[i][j] * f[i] + argument_vector[i] = prev_vector[i] - increment + if abs((argument_vector[0] - prev_vector[0])/argument_vector[0]) < eps: + break + elif abs((argument_vector[1] - prev_vector[1])/argument_vector[1]) < eps: + break + + print(argument_vector) + print(f1(argument_vector[0], argument_vector[1])) + print(f2(argument_vector[0], argument_vector[1])) + + +def build_jacobian(argument_vector, function_vector): + h = 1 + args = [0, 0] + f = [0, 0] + jacobian = [[0, 0], [0, 0]] + for i in range(2): + for j in range(2): + for k in range(2): + args[k] = argument_vector[k] + args[j] = argument_vector[j] + h + f[0] = f1(args[0], args[1]) + f[1] = f2(args[0], args[1]) + jacobian[i][j] = (f[i] - function_vector[i])/h + return jacobian + + +def initialize_matrix(): + matrix = [[0 for _ in range(4)] for _ in range(4)] + for i in range(len(matrix)): + for j in range(len(matrix[i])): + if i == j: + matrix[i][j] = 1 + else: + matrix[i][j] = 0 + return matrix + + +def reverse_matrix(matrix, identity_matrix): + size = len(matrix) + reversed_matrix = [[0 for _ in range(size)] for _ in range(size)] + for i in range(size): + solution = roots(i, matrix, identity_matrix) + for j in range(len(solution)): + reversed_matrix[j][i] = solution[j] + return reversed_matrix + + +def roots(counter, matrix, identity_matrix): + size = len(matrix) + coefficients = [[0 for _ in range(size)] for _ in range(size)] + free_members = [0 for _ in range(size)] + argument_positions = [0 for _ in range(size)] + result_coefficients = [[0 for _ in range(size)] for _ in range(size)] + result_free_members = [0 for _ in range(size)] + free_members, coefficients, argument_positions = initialize_system( + free_members, identity_matrix, counter, argument_positions, coefficients, matrix) + result_free_members, free_members, coefficients, result_coefficients = direct_way( + result_free_members, free_members, coefficients, result_coefficients, argument_positions, size) + result = reverse_way(result_free_members, result_coefficients) + result, argument_positions = order_vector(result, argument_positions) + return result + + +def order_vector(result, argument_positions): + for i in range(len(result)): + if argument_positions != i: + arg = argument_positions[i] + value = result[i] + result[i] = result[arg] + result[arg] = value + argument_positions[i] = argument_positions[arg] + argument_positions[arg] = arg + return result, argument_positions + + +def direct_way(result_free_members, free_members, coefficients, result_coefficients, argument_positions, size): + for i in range(size): + coefficients, free_members, argument_positions, result_coefficients = optimize_matrix( + i, coefficients, free_members, argument_positions, result_coefficients, size) + result_free_members[i] = free_members[i] / coefficients[i][i] + for j in range(i + 1, size): + free_members[j] = free_members[j] - coefficients[j][i] * result_free_members[i] + for k in range(i + 1, size): + result_coefficients[i][k] = coefficients[i][k] / coefficients[i][i] + coefficients[j][k] = coefficients[j][k] - coefficients[j][i] * result_coefficients[i][k] + return result_free_members, free_members, coefficients, result_coefficients + + +def reverse_way(result_free_members, result_coefficients): + size = len(result_free_members) + solution = [0 for _ in range(size)] + for i in range(size -1, -1, -1): + sum = 0.0 + for j in range(i + 1, size): + sum += result_coefficients[i][j] * solution[j] + solution[i] = result_free_members[i] - sum + return solution + + +def initialize_system(free_members, identity_matrix, counter, argument_positions, coefficients, matrix): + size = len(matrix) + for i in range(size): + free_members[i] = identity_matrix[i][counter] + argument_positions[i] = i + for j in range(size): + coefficients[i][j] = matrix[i][j] + return free_members, coefficients, argument_positions + + +def optimize_matrix(r, coefficients, free_members, argument_positions, result_coefficients, size): + max_coefficient = coefficients[r][r] + max_row = r + max_col = r + for i in range(size): + for j in range(size): + if max_coefficient < abs(coefficients[i][j]): + max_coefficient = abs(coefficients[i][j]) + max_row = i + max_col = j + free_members = swap_array_values(free_members, r, max_row) + for l in range(size): + coefficients = swap_matrix_values_row(coefficients, r, max_row, l) + argument_positions = swap_argument_positions(argument_positions, r, max_col) + for m in range(size): + if m < r: + result_coefficients = swap_matrix_values_columns(result_coefficients, m, r, max_col) + else: + coefficients = swap_matrix_values_columns(coefficients, m, r, max_col) + return coefficients, free_members, argument_positions, result_coefficients + + +def swap_matrix_values_columns(matrix, r, fc, sc): + temp = matrix[r][fc] + matrix[r][fc] = matrix[r][sc] + matrix[r][sc] = temp + return matrix + + +def swap_argument_positions(argument_positions, r, max_col): + temp = argument_positions[r] + argument_positions[r] = argument_positions[max_col] + argument_positions[max_col] = temp + return argument_positions + + +def swap_array_values(free_members, fc, sc): + temp = free_members[fc] + free_members[fc] = free_members[sc] + free_members[sc] = temp + return free_members + + +def swap_matrix_values_row(matrix, fc, sc, col): + temp = matrix[fc][col] + matrix[fc][col] = matrix[sc][col] + matrix[sc][col] = temp + return matrix + + +if __name__ == '__main__': + newton_jakobi()