From c0dfbe58ecb689885f6a9ab718870842ce5f7c74 Mon Sep 17 00:00:00 2001 From: Benjamin Berkels <berkels@aices.rwth-aachen.de> Date: Wed, 21 Apr 2021 12:03:53 +0200 Subject: [PATCH] added IPython notebook for Example 1.1.1 --- exp-111.ipynb | 135 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 exp-111.ipynb diff --git a/exp-111.ipynb b/exp-111.ipynb new file mode 100644 index 0000000..25b7d2d --- /dev/null +++ b/exp-111.ipynb @@ -0,0 +1,135 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### IPython notebook for Example 1.1.1 from the lecture" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import numpy as np\n", + "from numpy.linalg import lstsq, matrix_rank, norm\n", + "from IPython.display import display, Math\n", + "from sympy import Matrix, latex\n", + "\n", + "def my_latex(expr):\n", + " return latex(expr, mat_delim='(').replace('.0 &', '&').replace('.0\\\\', '\\\\')\n", + "\n", + "delta_t = 0.15\n", + "n = 6\n", + "t = np.arange(1, n + 1) * delta_t\n", + "xi_hat = np.array([1.2, 0.6, 1.6, 0.9])\n", + "xi = xi_hat[:-1]\n", + "\n", + "A_hat = np.stack((t, np.exp(t), t ** 3, np.sin(t)), axis=1)\n", + "A = A_hat[:, :-1]\n", + "display(Math(r'\\hat A = h \\begin{pmatrix} t_1 & e^{t_1} & t_1^3 & \\sin(t_1) \\\\ t_2 & e^{t_2} & t_2^3 & \\sin(t_2) \\\\ \\vdots &\\vdots & \\vdots & \\vdots \\end{pmatrix},\\quad ' + \\\n", + " r'A = h \\begin{pmatrix} t_1 & e^{t_1} & t_1^3 \\\\ t_2 & e^{t_2} & t_2^3 \\\\ \\vdots &\\vdots & \\vdots \\end{pmatrix},\\quad ' + \\\n", + " rf't = {my_latex(Matrix(t))}'))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "display(Math(r\"\\operatorname{rank}(\\hat A) = \" + f\"{matrix_rank(A_hat)}\"))\n", + "display(Math(r\"\\operatorname{rank}(A) = \" + f\"{matrix_rank(A)}\"))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "display(Math(r\"\\hat y=\\hat A\\hat\\xi,\\quad y=A\\xi, \\quad \" + rf\"\\hat\\xi = {my_latex(Matrix(xi_hat))},\\quad \\xi = {my_latex(Matrix(xi))}\"))\n", + "y_hat = np.matmul(A_hat, xi_hat)\n", + "y = np.matmul(A, xi)\n", + "\n", + "delta_y1 = 1e-2 * np.array([1, 0, -1, -1, -0.5, 1])\n", + "delta_y2 = 1e-2 * np.array([-1, 1, 1, -0.5, -2, 1])\n", + "\n", + "display(Math(rf\"\\delta y_1 = {my_latex(Matrix(delta_y1))},\\quad \\delta y_2 = {my_latex(Matrix(delta_y2))}\"))\n", + "display(Math(f\"\\|\\delta y_1\\|_2 = {norm(delta_y1):.3f}\" + r\",\\quad\" \\\n", + " f\"\\|\\delta y_2\\|_2 = {norm(delta_y2):.3f}\"))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "tilde_y1_hat = y_hat + delta_y1\n", + "tilde_y2_hat = y_hat + delta_y2\n", + "tilde_y1 = y + delta_y1\n", + "tilde_y2 = y + delta_y2\n", + "\n", + "tilde_xi1 = lstsq(A, tilde_y1, rcond=None)[0]\n", + "tilde_xi2 = lstsq(A, tilde_y2, rcond=None)[0]\n", + "\n", + "display(Math(r\"\\tilde \\xi_i=\\operatorname{arg\\,min}\\|A\\xi-(y+\\delta y_i)\\|\"))\n", + "\n", + "display(Math(rf\"\\|\\tilde \\xi_1-\\xi\\|_2 = {norm(tilde_xi1 - xi):.3f}\" + r\",\\quad\" \\\n", + " rf\"\\|\\tilde \\xi_2-\\xi\\|_2 = {norm(tilde_xi2 - xi):.3f}\"))\n", + "\n", + "tilde_xi1_hat = lstsq(A_hat, tilde_y1_hat, rcond=None)[0]\n", + "tilde_xi2_hat = lstsq(A_hat, tilde_y2_hat, rcond=None)[0]\n", + "\n", + "display(Math(r\"\\tilde{\\hat{\\xi}}_i=\\operatorname{arg\\,min}\\|\\hat A\\xi-(\\hat y+\\delta y_i)\\|\"))\n", + "\n", + "display(Math(r\"\\|\\tilde{\\hat{\\xi}}_1-\\hat{\\xi}\\|_2 = \" + f\"{norm(tilde_xi1_hat - xi_hat):.3f}\" + r\",\\quad\" \\\n", + " r\"\\|\\tilde{\\hat{\\xi}}_2-\\hat{\\xi}\\|_2 = \" + f\"{norm(tilde_xi2_hat - xi_hat):.1f}\"))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "display(Math(r\"\\frac{\\|\\tilde \\xi_1-\\xi\\|_2}{\\|\\delta y_1\\|_2} = \" + f\"{norm(tilde_xi1 - xi) / norm(delta_y1):.1f}\" + r\",\\quad\" \\\n", + " r\"\\frac{\\|\\tilde{\\hat{\\xi}}_1-\\hat{\\xi}\\|_2}{\\|\\delta y_1\\|_2} = \" + f\"{norm(tilde_xi1_hat - xi_hat) / norm(delta_y1):.1f}\"))\n", + "display(Math(r\"\\frac{\\|\\tilde \\xi_2-\\xi\\|_2}{\\|\\delta y_2\\|_2} = \" + f\"{norm(tilde_xi2 - xi) / norm(delta_y2):.1f}\" + r\",\\quad\" \\\n", + " r\"\\frac{\\|\\tilde{\\hat{\\xi}}_2-\\hat{\\xi}\\|_2}{\\|\\delta y_2\\|_2} = \" + f\"{norm(tilde_xi2_hat - xi_hat) / norm(delta_y2):.1e}\"))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "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.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file -- GitLab