Upload New File

352e4078 · Benjamin Berkels · 178acef1 · 352e4078
Commit 352e4078 authored 5 years ago by Benjamin Berkels
--- a/exp-232.ipynb
+++ b/exp-232.ipynb
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### IPython notebook for Example 2.3.2 from the lecture"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "import numpy as np\n",
+        "from IPython.display import display, Math\n",
+        "\n",
+        "display(Math(r'''A=\\begin{pmatrix} 1 & 1 \\\\ 0 & 0 \\\\ 0 & 1 \\end{pmatrix}'''))\n",
+        "A = np.array([[1., 1.], [0., 0.], [0., 1.]])\n",
+        "display(Math(r'''y= \\begin{pmatrix} 0.01 \\\\ 1 \\\\ 0 \\end{pmatrix}'''))\n",
+        "y = np.array([0.01, 1., 0.])\n",
+        "display(Math(r'''\\delta y= \\begin{pmatrix} 0 \\\\ 0 \\\\ 0.01 \\end{pmatrix}'''))\n",
+        "delta_y = np.array([0., 0., 0.01])"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "print('rank(A) =', np.linalg.matrix_rank(A))\n",
+        "print('cond(A) = {:.2f}'.format(np.linalg.cond(A)))"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "display(Math(r'U^TAV=\\Sigma'))\n",
+        "U, sigma, VT = np.linalg.svd(A)\n",
+        "display(Math(r'P_Ay=(u_1^Ty)u_1+(u_2^Ty)u_2'))\n",
+        "PAy = np.dot(U[:, 0], y)*U[:, 0] + np.dot(U[:, 1], y)*U[:, 1]\n",
+        "display(Math(r'\\text{We have }P_Ay=\\begin{pmatrix}0.01 \\\\ 0 \\\\ 0 \\end{pmatrix}\\text{, since }\\text{range}(A)= \\text{span} \\{e_1,e_3\\}'))"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "print('||P_Ay||/||y|| = {:.2f}'.format(np.linalg.norm(PAy)/np.linalg.norm(y)))\n",
+        "data_error = np.linalg.norm(delta_y)/np.linalg.norm(y)\n",
+        "print('||delta y||/||y|| = {:.2f}'.format(data_error))"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "A_dagger = np.linalg.pinv(A)\n",
+        "\n",
+        "display(Math(r'''\\xi = A^\\dagger y'''))\n",
+        "xi = np.matmul(A_dagger, y)\n",
+        "display(Math(r'''\\tilde\\xi = A^\\dagger (y+\\tilde y)'''))\n",
+        "xi_tilde = np.matmul(A_dagger, y + delta_y)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "solution_error = np.linalg.norm(xi_tilde-xi)/np.linalg.norm(xi)\n",
+        "print('||xi_tilde - xi||/||xi|| = {:.2f}'.format(solution_error))\n",
+        "print('Q(y,delta_y) = {:.2f}'.format(solution_error/data_error))\n"
+      ],
+      "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": 1
+}
\ No newline at end of file
+%% Cell type:markdown id: tags:
+
+### IPython notebook for Example 2.3.2 from the lecture
+
+%% Cell type:code id: tags:
+
+``` python
+import numpy as np
+from IPython.display import display, Math
+
+display(Math(r'''A=\begin{pmatrix} 1 & 1 \\ 0 & 0 \\ 0 & 1 \end{pmatrix}'''))
+A = np.array([[1., 1.], [0., 0.], [0., 1.]])
+display(Math(r'''y= \begin{pmatrix} 0.01 \\ 1 \\ 0 \end{pmatrix}'''))
+y = np.array([0.01, 1., 0.])
+display(Math(r'''\delta y= \begin{pmatrix} 0 \\ 0 \\ 0.01 \end{pmatrix}'''))
+delta_y = np.array([0., 0., 0.01])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+print('rank(A) =', np.linalg.matrix_rank(A))
+print('cond(A) = {:.2f}'.format(np.linalg.cond(A)))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+display(Math(r'U^TAV=\Sigma'))
+U, sigma, VT = np.linalg.svd(A)
+display(Math(r'P_Ay=(u_1^Ty)u_1+(u_2^Ty)u_2'))
+PAy = np.dot(U[:, 0], y)*U[:, 0] + np.dot(U[:, 1], y)*U[:, 1]
+display(Math(r'\text{We have }P_Ay=\begin{pmatrix}0.01 \\ 0 \\ 0 \end{pmatrix}\text{, since }\text{range}(A)= \text{span} \{e_1,e_3\}'))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+print('||P_Ay||/||y|| = {:.2f}'.format(np.linalg.norm(PAy)/np.linalg.norm(y)))
+data_error = np.linalg.norm(delta_y)/np.linalg.norm(y)
+print('||delta y||/||y|| = {:.2f}'.format(data_error))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+A_dagger = np.linalg.pinv(A)
+
+display(Math(r'''\xi = A^\dagger y'''))
+xi = np.matmul(A_dagger, y)
+display(Math(r'''\tilde\xi = A^\dagger (y+\tilde y)'''))
+xi_tilde = np.matmul(A_dagger, y + delta_y)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+solution_error = np.linalg.norm(xi_tilde-xi)/np.linalg.norm(xi)
+print('||xi_tilde - xi||/||xi|| = {:.2f}'.format(solution_error))
+print('Q(y,delta_y) = {:.2f}'.format(solution_error/data_error))
+```