diff --git a/exp-355.ipynb b/exp-355.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..1abe5e63d0a55871a8a7a8d7c7bb8b09b5bf1f82 --- /dev/null +++ b/exp-355.ipynb @@ -0,0 +1,179 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### IPython notebook for Example 3.5.5 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)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "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)\n", + "epsilon_a = np.linalg.norm(y_tilde-y)\n", + "print('epsilon_a=||y_tilde-y|| = {:.3f}'.format(epsilon_a))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "discrepancy = np.zeros(100)\n", + "alphas = np.zeros(100)\n", + "xi_alphas = np.zeros((100, n))\n", + "for i in range(1, 101):\n", + " alpha = 0.005*i\n", + " A_Tikh = np.concatenate((A, np.diag(alpha*np.ones(n))), axis=0)\n", + " y_tilde_Tikh = np.concatenate((y_tilde, np.zeros(n)), axis=0)\n", + " xi_alpha = np.linalg.lstsq(A_Tikh, y_tilde_Tikh, rcond=0)[0]\n", + " xi_alphas[i-1, :] = xi_alpha\n", + " alphas[i-1] = alpha\n", + " discrepancy[i-1] = np.linalg.norm(np.matmul(A, xi_alpha)-y_tilde)\n", + "\n", + "set_matplotlib_formats('svg')\n", + "plt.title(\"Tikhonov, discrepancy\")\n", + "plt.xlabel(r'$\\alpha$')\n", + "plt.ylabel(r'$||A\\tilde \\xi_{\\alpha}- \\tilde y||_2$')\n", + "plt.plot(alphas, discrepancy, 'k+')\n", + "plt.hlines(epsilon_a, xmin=alphas[0], xmax=alphas[-1])\n", + "plt.show()\n", + "plt.close()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "alpha_index = np.argmin(np.abs(discrepancy-epsilon_a))\n", + "print('alpha_d = {}'.format(alphas[alpha_index]))\n", + "\n", + "j = np.arange(1, n+1)\n", + "plt.title(\"Tikhonov, alpha_d = {}\".format(alphas[alpha_index]))\n", + "plt.plot(j, xi_alphas[alpha_index, :], 'kx')\n", + "plt.show()\n", + "plt.close()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "# CGNE\n", + "xi_tilde = np.zeros(n)\n", + "xi_ks = np.zeros((101, 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 = np.zeros(101)\n", + "discrepancy[0] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)\n", + "for k in range(1, 101):\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[k] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)\n", + " xi_ks[k, :] = xi_tilde\n", + "\n", + "plt.title(\"CGNE, discrepancy\")\n", + "plt.xlabel(r'$k$')\n", + "plt.ylabel(r'$||A\\tilde \\xi_k- \\xi||_2$')\n", + "plt.yscale('log')\n", + "plt.plot(discrepancy, 'k+')\n", + "plt.hlines(epsilon_a, xmin=0, xmax=discrepancy.shape[0])\n", + "plt.show()\n", + "plt.close()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "k_opt = np.argmin(np.abs(discrepancy-epsilon_a))\n", + "plt.title(\"CGNE, k_opt = {}\".format(k_opt))\n", + "plt.plot(j, xi_ks[k_opt, :], 'kx')\n", + "plt.show()\n", + "plt.close()\n" + ], + "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": 1 +} \ No newline at end of file