You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
184 lines
5.6 KiB
Plaintext
184 lines
5.6 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 15,
|
|
"id": "c8de76b1-c455-47ff-91c1-11abed30fcff",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import itertools"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 16,
|
|
"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",
|
|
"\n",
|
|
"# Construct the required matrix\n",
|
|
"\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": 17,
|
|
"id": "1fb8aa25-9a4f-408a-8705-3dd9d91c7089",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Solve Lx = f, where L is lower triangular\n",
|
|
"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",
|
|
" # Solve Ux = f, where U is upper triangular\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, where 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": 21,
|
|
"id": "6e399987-73ed-452c-927b-cc9908786520",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def do_first_task(n: int) -> None:\n",
|
|
" a = get_matrix_a(n)\n",
|
|
" f = get_f(n)\n",
|
|
" # HERE\n",
|
|
" l = crudepy.linalg.cholesky(a)\n",
|
|
" x = lu_solve(l, l.T, f)\n",
|
|
"\n",
|
|
" print(f\"a = \\n{a}\")\n",
|
|
" print(f\"f = \\n{f}\")\n",
|
|
" print(\"Solution: -----\")\n",
|
|
" print(f\"L = \\n{l}\")\n",
|
|
" print(f\"x = \\n{x}\")\n",
|
|
"\n",
|
|
" # Find the exact solution with numpy and compare results\n",
|
|
" exact_x = np.linalg.solve(a, f)\n",
|
|
" delta = exact_x - x\n",
|
|
"\n",
|
|
" print(f\"|delta|_1 = {crudepy.linalg.vec_norm_max(delta):#e}\")\n",
|
|
" print(f\"|delta|_2 = {crudepy.linalg.vec_norm_sum(delta):#e}\")\n",
|
|
" print(f\"|delta|_3 = {crudepy.linalg.vec_norm_euclid(delta):#e}\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 19,
|
|
"id": "47569073-4acd-42c0-a2dd-6c8c2deee9b4",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"n = 6"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 20,
|
|
"id": "300cd1ef-cbee-48fb-89ef-dd3779da1b24",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"ename": "NameError",
|
|
"evalue": "name 'crudepy' is not defined",
|
|
"output_type": "error",
|
|
"traceback": [
|
|
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
|
|
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
|
|
"Cell \u001b[0;32mIn[20], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mdo_first_task\u001b[49m\u001b[43m(\u001b[49m\u001b[43mn\u001b[49m\u001b[43m)\u001b[49m\n",
|
|
"Cell \u001b[0;32mIn[18], line 4\u001b[0m, in \u001b[0;36mdo_first_task\u001b[0;34m(n)\u001b[0m\n\u001b[1;32m 2\u001b[0m a \u001b[38;5;241m=\u001b[39m get_matrix_a(n)\n\u001b[1;32m 3\u001b[0m f \u001b[38;5;241m=\u001b[39m get_f(n)\n\u001b[0;32m----> 4\u001b[0m l \u001b[38;5;241m=\u001b[39m \u001b[43mcrudepy\u001b[49m\u001b[38;5;241m.\u001b[39mlinalg\u001b[38;5;241m.\u001b[39mcholesky(a)\n\u001b[1;32m 5\u001b[0m x \u001b[38;5;241m=\u001b[39m lu_solve(l, l\u001b[38;5;241m.\u001b[39mT, f)\n\u001b[1;32m 7\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124ma = \u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;132;01m{\u001b[39;00ma\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n",
|
|
"\u001b[0;31mNameError\u001b[0m: name 'crudepy' is not defined"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"do_first_task(n)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "9f33c8b8-c89f-43a8-9dfb-da5313ad0045",
|
|
"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
|
|
}
|