From a0765e669c9ee08e3c4fbd991ac67479dbbb23e2 Mon Sep 17 00:00:00 2001
From: Benjamin Berkels <berkels@aices.rwth-aachen.de>
Date: Tue, 19 May 2020 14:28:30 +0200
Subject: [PATCH] Upload New File
---
exp-357.ipynb | 197 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 197 insertions(+)
create mode 100644 exp-357.ipynb
diff --git a/exp-357.ipynb b/exp-357.ipynb
new file mode 100644
index 0000000..3c61870
--- /dev/null
+++ b/exp-357.ipynb
@@ -0,0 +1,197 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### IPython notebook for Example 3.5.7 from the lecture"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from IPython.display import set_matplotlib_formats, display, Math\n",
+ "\n",
+ "\n",
+ "def k(s, t, d):\n",
+ " return d / np.power(d**2 + (s-t)**2, 3/2)\n",
+ "\n",
+ "\n",
+ "d = 0.1\n",
+ "n = 80\n",
+ "h = 1/n\n",
+ "\n",
+ "display(Math(r'\\text{Create }A\\text{ from Example 1.7.1, }\\xi_j=\\begin{cases}2&1\\leq j\\leq40\\\\1&41\\leq j\\leq80\\end{cases},\\ y=A\\xi'))\n",
+ "A = np.zeros((n, n))\n",
+ "for i in range(n):\n",
+ " for j in range(n):\n",
+ " A[i, j] = h * k((i+0.5)*h, (j+0.5)*h, d)\n",
+ "\n",
+ "xi = np.ones(n)\n",
+ "xi[:n//2] = 2\n",
+ "y = np.matmul(A, xi)\n",
+ "\n",
+ "display(Math(r'\\tilde y= y+\\delta y'))\n",
+ "np.random.seed(0)\n",
+ "y_tilde = y + 0.02*(np.random.rand(n)-0.5)"
+ ],
+ "outputs": [],
+ "execution_count": null
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "source": [
+ "discrepancy_tilde = np.zeros(100)\n",
+ "xi_tilde_norm = np.zeros(100)\n",
+ "discrepancy = np.zeros(100)\n",
+ "xi_norm = np.zeros(100)\n",
+ "for i in range(1, 101):\n",
+ " alpha = 0.005*i\n",
+ " A_Tikh = np.concatenate((A, alpha*np.eye(n)), axis=0)\n",
+ " y_tilde_Tikh = np.concatenate((y_tilde, np.zeros(n)), axis=0)\n",
+ " xi_tilde_alpha = np.linalg.lstsq(A_Tikh, y_tilde_Tikh, rcond=0)[0]\n",
+ " discrepancy_tilde[i-1] = np.linalg.norm(np.matmul(A, xi_tilde_alpha)-y_tilde)\n",
+ " xi_tilde_norm[i-1] = np.linalg.norm(xi_tilde_alpha)\n",
+ " y_Tikh = np.concatenate((y, np.zeros(n)), axis=0)\n",
+ " xi_alpha = np.linalg.lstsq(A_Tikh, y_Tikh, rcond=0)[0]\n",
+ " discrepancy[i-1] = np.linalg.norm(np.matmul(A, xi_alpha)-y)\n",
+ " xi_norm[i-1] = np.linalg.norm(xi_alpha)\n",
+ "\n",
+ "set_matplotlib_formats('svg')\n",
+ "plt.title(\"Tikhonov L-curve\")\n",
+ "plt.xlabel(r'$||A\\tilde \\xi_{\\alpha}- \\tilde y||_2^2$')\n",
+ "plt.ylabel(r'$||\\xi_{\\alpha}||_2^2$')\n",
+ "plt.plot(discrepancy_tilde**2, xi_tilde_norm**2, 'k+')\n",
+ "plt.show()\n",
+ "plt.close()"
+ ],
+ "outputs": [],
+ "execution_count": null
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "source": [
+ "plt.title(\"Tikhonov L-curves, zoom\")\n",
+ "plt.plot(discrepancy_tilde[:50]**2, xi_tilde_norm[:50]**2, 'k+', label=r'$(\\|A \\tilde \\xi_\\alpha - \\tilde y\\|_2^2, \\|\\tilde \\xi_\\alpha\\|_2^2)$')\n",
+ "plt.plot(discrepancy[:50]**2, xi_norm[:50]**2, 'k.', label=r'$(\\|A \\xi_\\alpha - y\\|_2^2, \\| \\xi_\\alpha\\|_2^2)$')\n",
+ "plt.legend()\n",
+ "plt.show()\n",
+ "plt.close()"
+ ],
+ "outputs": [],
+ "execution_count": null
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "source": [
+ "plt.title(\"Tikhonov L-curves, zoom\")\n",
+ "plt.plot(discrepancy_tilde[:50], xi_tilde_norm[:50], 'k+', label=r'$(\\|A \\tilde \\xi_\\alpha - \\tilde y\\|_2, \\|\\tilde \\xi_\\alpha\\|_2)$')\n",
+ "plt.plot(discrepancy[:50], xi_norm[:50], 'k.', label=r'$(\\|A \\xi_\\alpha - y\\|_2, \\| \\xi_\\alpha\\|_2)$')\n",
+ "plt.legend()\n",
+ "plt.show()\n",
+ "plt.close()"
+ ],
+ "outputs": [],
+ "execution_count": null
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "source": [
+ "def CGNE(A, y_tilde):\n",
+ " xi_tilde = np.zeros(n)\n",
+ " # xi_tilde is zero, so the initial residuum is the input data.\n",
+ " r = y_tilde.copy()\n",
+ " d = np.matmul(A.T, r)\n",
+ " p = d.copy()\n",
+ " p_norm = np.linalg.norm(p)\n",
+ " discrepancy_tilde = np.zeros(201)\n",
+ " discrepancy_tilde[0] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)\n",
+ " xi_tilde_norm = np.zeros(201)\n",
+ " for k in range(1, 201):\n",
+ " q = np.matmul(A, d)\n",
+ " beta = (np.linalg.norm(p)/np.linalg.norm(q))**2\n",
+ " xi_tilde += beta * d\n",
+ " r += -beta*q\n",
+ " p = np.matmul(A.T, r)\n",
+ " p_norm_new = np.linalg.norm(p)\n",
+ " gamma = (p_norm_new/p_norm)**2\n",
+ " p_norm = p_norm_new\n",
+ " d = p + gamma*d\n",
+ " discrepancy_tilde[k] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)\n",
+ " xi_tilde_norm[k] = np.linalg.norm(xi_tilde)\n",
+ " return discrepancy_tilde, xi_tilde_norm\n",
+ "\n",
+ "\n",
+ "discrepancy_tilde, xi_tilde_norm = CGNE(A, y_tilde)\n",
+ "\n",
+ "plt.title(\"CGNE L-curve\")\n",
+ "plt.xlabel(r'$||A\\tilde \\xi_{k}- \\tilde y||_2$')\n",
+ "plt.ylabel(r'$||\\xi_{k}||_2$')\n",
+ "plt.plot(discrepancy_tilde[2:], xi_tilde_norm[2:], 'k+')\n",
+ "plt.show()\n",
+ "plt.close()\n"
+ ],
+ "outputs": [],
+ "execution_count": null
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "source": [
+ "plt.title(\"CGNE L-curve log-log scale\")\n",
+ "plt.xlabel(r'$||A\\tilde \\xi_{k}- \\tilde y||_2$')\n",
+ "plt.ylabel(r'$||\\xi_{k}||_2$')\n",
+ "plt.loglog(discrepancy_tilde[2:], xi_tilde_norm[2:], 'k+')\n",
+ "plt.show()\n",
+ "plt.close()"
+ ],
+ "outputs": [],
+ "execution_count": null
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "source": [
+ "discrepancy, xi_norm = CGNE(A, y)\n",
+ "\n",
+ "plt.title(\"CGNE L-curve log-log scale, zoom\")\n",
+ "plt.loglog(discrepancy_tilde[10:80], xi_tilde_norm[10:80], 'k+', label=r'$(\\|A \\tilde \\xi_\\alpha - \\tilde y\\|_2, \\|\\tilde \\xi_\\alpha\\|_2)$')\n",
+ "plt.loglog(discrepancy[10:80], xi_norm[10:80], 'k.', label=r'$(\\|A \\xi_\\alpha - y\\|_2, \\| \\xi_\\alpha\\|_2)$')\n",
+ "plt.legend()\n",
+ "plt.show()\n",
+ "plt.close()\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
--
GitLab