{ "cells": [ { "cell_type": "code", "execution_count": 20, "id": "c8de76b1-c455-47ff-91c1-11abed30fcff", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import itertools" ] }, { "cell_type": "code", "execution_count": 21, "id": "b6479073-0cff-478f-868d-26fedd310af5", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "def get_a_ij(i: int, j: int, n: int):\n", " return 10.0 / (i + j + 1)\n", "\n", "def get_matrix_a(n: int) -> np.array:\n", " a = np.zeros((n, n))\n", " from_0_to_n = [i for i in range(n)]\n", " for i, j in itertools.product(from_0_to_n, from_0_to_n):\n", " a[i, j] = get_a_ij(i, j, n)\n", " return a\n", "\n", "\n", "def get_f(n: int) -> np.array:\n", " return np.sum(get_matrix_a(n), axis=1)" ] }, { "cell_type": "code", "execution_count": 22, "id": "1fb8aa25-9a4f-408a-8705-3dd9d91c7089", "metadata": {}, "outputs": [], "source": [ "def lower_solve(l: np.array, f: np.array):\n", " assert l.ndim == 2 and l.shape[0] == l.shape[1]\n", " dim = l.shape[0]\n", "\n", " x = np.zeros(dim)\n", " for i in range(dim):\n", " tmp = f[i]\n", " for j in range(i):\n", " tmp -= l[i, j] * x[j]\n", " x[i] = tmp / l[i, i]\n", "\n", " return x\n", "\n", "\n", "def upper_solve(u: np.array, f: np.array):\n", " assert u.ndim == 2 and u.shape[0] == u.shape[1]\n", " dim = u.shape[0]\n", "\n", " x = np.zeros(dim)\n", " for i in range(dim - 1, -1, -1):\n", " tmp = f[i]\n", " for j in range(i + 1, dim):\n", " tmp -= u[i, j] * x[j]\n", " x[i] = tmp / u[i, i]\n", "\n", " return x\n", "\n", "\n", "def lu_solve(l, u, f):\n", " # Ax = f, где A = LU\n", " # L(Ux) = f, Ux = y <=> y = Ux\n", " y = lower_solve(l, f)\n", " x = upper_solve(u, y)\n", " return x" ] }, { "cell_type": "code", "execution_count": 23, "id": "81937d82-febb-46b0-954d-e2d631f124b8", "metadata": {}, "outputs": [], "source": [ "def cholesky(A):\n", " n = A.shape[0]\n", " L = np.zeros_like(A)\n", " for i in range(n):\n", " for j in range(i+1):\n", " s = 0\n", " for k in range(j):\n", " s += L[i][k] * L[j][k]\n", "\n", " if (i == j):\n", " L[i][j] = (A[i][i] - s) ** 0.5\n", " else:\n", " L[i][j] = (1.0 / L[j][j] * (A[i][j] - s))\n", " return L" ] }, { "cell_type": "code", "execution_count": 24, "id": "06916637-b82c-4e0b-9792-6558bced33b8", "metadata": {}, "outputs": [], "source": [ "# Векторные и подчиненные им матричные нормы\n", "\n", "def vec_norm_max(vector: np.array):\n", " # |v|_1 = max(abs(v_i))\n", " return np.abs(vector).max()\n", "\n", "\n", "def vec_norm_sum(vector: np.array):\n", " # |v|_2 = sum(abs(v_i))\n", " return np.abs(vector).sum()\n", "\n", "def vec_norm_euclid(vector: np.array):\n", " return np.sqrt((vector ** 2.0).sum())\n", "\n", "\n", "def row_sum_norm(matrix: np.array):\n", " # |A|_1 = max( sum(a_ij) for j ) over 1 <= i <= n\n", " return np.sum(np.abs(matrix), axis=1).max()\n", "\n", "\n", "def column_sum_norm(matrix: np.array):\n", " # |A|_2 = max( sum(a_ij) for i ) over 1 <= j <= n\n", " return np.sum(np.abs(matrix), axis=0).max()\n", "\n", "\n", "def spectral_norm(matrix: np.array):\n", " matrix_h = matrix.conj().T\n", " eigenvalues = np.linalg.eigvals(np.matmul(matrix_h, matrix))\n", " return np.sqrt(eigenvalues.max())\n", "\n", "\n", "def mu(matrix_a: int, norm):\n", " inv_matrix_a = np.linalg.inv(matrix_a)\n", " norm_a = norm(matrix_a)\n", " norm_inv_a = norm(inv_matrix_a)\n", " return norm_a * norm_inv_a\n" ] }, { "cell_type": "code", "execution_count": 25, "id": "6e399987-73ed-452c-927b-cc9908786520", "metadata": {}, "outputs": [], "source": [ "def task1(n: int) -> None:\n", " a = get_matrix_a(n)\n", " f = get_f(n)\n", " l = cholesky(a)\n", " x = lu_solve(l, l.T, f)\n", "\n", " print(\"Исходные данные:\")\n", " print(f\"a = \\n{a}\")\n", " print(f\"f = \\n{f}\")\n", " print(\"Решение:\")\n", " print(f\"x = \\n{x}\")\n", "\n", " exact_x = np.linalg.solve(a, f)\n", " delta = exact_x - x\n", "\n", " print(f\"|delta|_1 = {vec_norm_max(delta):#e}\")\n", " print(f\"|delta|_2 = {vec_norm_sum(delta):#e}\")\n", " print(f\"|delta|_3 = {vec_norm_euclid(delta):#e}\")" ] }, { "cell_type": "code", "execution_count": 26, "id": "47569073-4acd-42c0-a2dd-6c8c2deee9b4", "metadata": {}, "outputs": [], "source": [ "n = 6" ] }, { "cell_type": "code", "execution_count": 27, "id": "300cd1ef-cbee-48fb-89ef-dd3779da1b24", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Исходные данные:\n", "a = \n", "[[10. 5. 3.33333333 2.5 2. 1.66666667]\n", " [ 5. 3.33333333 2.5 2. 1.66666667 1.42857143]\n", " [ 3.33333333 2.5 2. 1.66666667 1.42857143 1.25 ]\n", " [ 2.5 2. 1.66666667 1.42857143 1.25 1.11111111]\n", " [ 2. 1.66666667 1.42857143 1.25 1.11111111 1. ]\n", " [ 1.66666667 1.42857143 1.25 1.11111111 1. 0.90909091]]\n", "f = \n", "[24.5 15.92857143 12.17857143 9.95634921 8.45634921 7.36544012]\n", "Решение:\n", "x = \n", "[1. 1. 1. 1. 1. 1.]\n", "|delta|_1 = 2.136471e-10\n", "|delta|_2 = 5.829489e-10\n", "|delta|_3 = 3.120671e-10\n" ] } ], "source": [ "task1(n)" ] }, { "cell_type": "code", "execution_count": null, "id": "0dab679e-f64b-482d-864d-2943f0e793a4", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5" } }, "nbformat": 4, "nbformat_minor": 5 }