diff --git a/Lab-1/viktorme-report.ipynb b/Lab-1/viktorme-report.ipynb new file mode 100644 index 0000000..c650812 --- /dev/null +++ b/Lab-1/viktorme-report.ipynb @@ -0,0 +1,374 @@ +{ + "nbformat": 4, + "nbformat_minor": 0, + "metadata": { + "colab": { + "name": "lab1_matrixalgorithms.ipynb", + "provenance": [], + "collapsed_sections": [] + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + } + }, + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "6RgtXlfYO_i7", + "colab_type": "text" + }, + "source": [ + "# **Lab 1: Matrix algorithms**\n", + "**Viktor Meyer - DD2363 Methods in Scientific Computing**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "9x_J5FVuPzbm", + "colab_type": "text" + }, + "source": [ + "# **Abstract**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "6UFTSzW7P8kL", + "colab_type": "text" + }, + "source": [ + "This lab is an exercise in matrix operations. The mandatory part includes scalar product, matrix-vector product and matrix-matrix product. There is also an extra assignment euclidean norm or euclidian distance." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "28xLGz8JX3Hh", + "colab_type": "text" + }, + "source": [ + "# **Environment**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "D2PYNusD08Wa", + "colab_type": "text" + }, + "source": [ + "To have access to the neccessary modules you have to run this cell." + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "Xw7VlErAX7NS", + "colab_type": "code", + "colab": {} + }, + "source": [ + "# Load neccessary modules.\n", + "from google.colab import files\n", + "\n", + "import time\n", + "import numpy as np\n", + "\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib import tri\n", + "from matplotlib import axes\n", + "from mpl_toolkits.mplot3d import Axes3D" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "gnO3lhAigLev", + "colab_type": "text" + }, + "source": [ + "# **Introduction**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "l5zMzgPlRAF6", + "colab_type": "text" + }, + "source": [ + "The methods to be implemented were covered in the first lectures and definitions are available in the lecture notes: LN-DD2363-part1.pdf\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "WeFO9QMeUOAu", + "colab_type": "text" + }, + "source": [ + "# **Methods**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "I5Kf3qRyY5J5", + "colab_type": "text" + }, + "source": [ + "### Scalar Product ###" + ] + }, + { + "cell_type": "code", + "metadata": { + "colab_type": "code", + "outputId": "5463ca4e-73b9-40c0-bd06-958e76b35400", + "id": "PbnT724f5mfQ", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 68 + } + }, + "source": [ + "# LN-DD2363-part1.pdf, pp. 7\n", + "x = np.array([0, 3, -7])\n", + "y = np.array([2, 3, 1])\n", + "\n", + "def scalar(x, y):\n", + " result = 0\n", + " assert(x.size == y.size)\n", + " for i in range(x.size):\n", + " result += x[i]*y[i]\n", + " \n", + " return result\n", + "\n", + "print(\"x = \", x)\n", + "print(\"y = \", y)\n", + "print('scalar(x, y) =',scalar(x, y))" + ], + "execution_count": 65, + "outputs": [ + { + "output_type": "stream", + "text": [ + "x = [ 0 3 -7]\n", + "y = [2 3 1]\n", + "scalar(x, y) = 2\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "colab_type": "text", + "id": "k9b0yn9qF_Jq" + }, + "source": [ + "### Matrix-Vector product ###" + ] + }, + { + "cell_type": "code", + "metadata": { + "colab_type": "code", + "outputId": "2ad783ee-22d5-46bb-a184-44e6036d1ee1", + "id": "401b1-yLGFaH", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 85 + } + }, + "source": [ + "# LN-DD2363-part1.pdf, pp. 22\n", + "x = np.array([2, 1, 0])\n", + "A = np.matrix(\"1 -1 2; 0 -3 1\")\n", + "\n", + "def matrixvectorproduct(x, A):\n", + " b = np.zeros(A.shape[0])\n", + " assert(A.shape[1] == x.size)\n", + " for colIndex in range(A.shape[0]):\n", + " for rowIndex in range(A.shape[1]):\n", + " b[colIndex] += x[rowIndex]*A[colIndex, rowIndex]\n", + " return b\n", + "\n", + "print(\"x = \", x)\n", + "print(\"A = \", A)\n", + "print('matrixvectorproduct(x, A) =', matrixvectorproduct(x, A))" + ], + "execution_count": 66, + "outputs": [ + { + "output_type": "stream", + "text": [ + "x = [2 1 0]\n", + "A = [[ 1 -1 2]\n", + " [ 0 -3 1]]\n", + "matrixvectorproduct(x, A) = [ 1. -3.]\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "colab_type": "text", + "id": "lVE-mL2SDDZu" + }, + "source": [ + "### Matrix-Matrix product ###" + ] + }, + { + "cell_type": "code", + "metadata": { + "colab_type": "code", + "outputId": "8e285832-68da-43ad-bde6-4f37d93b265a", + "id": "BBZekOW0DGxN", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 136 + } + }, + "source": [ + "# LN-DD2363-part1.pdf, pp. 24\n", + "A = np.matrix(\"0 4 -2; -4 -3 0\")\n", + "B = np.matrix(\"0 1; 1 -1; 2 3\")\n", + "\n", + "def matrixmatrixproduct(A, B):\n", + " C = np.matrix(np.zeros((A.shape[0], B.shape[1])))\n", + " assert(A.shape[1] == B.shape[0])\n", + " for rowIndex in range(A.shape[0]):\n", + " for colIndex in range(B.shape[1]):\n", + " for offset in range(A.shape[1]):\n", + " C[rowIndex, colIndex] += A[rowIndex, offset]*B[offset, colIndex]\n", + " return C\n", + "\n", + "print(\"A = \", A)\n", + "print(\"B = \", B)\n", + "print('matrixmatrixproduct(A, B) =', matrixmatrixproduct(A, B))" + ], + "execution_count": 67, + "outputs": [ + { + "output_type": "stream", + "text": [ + "A = [[ 0 4 -2]\n", + " [-4 -3 0]]\n", + "B = [[ 0 1]\n", + " [ 1 -1]\n", + " [ 2 3]]\n", + "matrixmatrixproduct(A, B) = [[ 0. -10.]\n", + " [ -3. -1.]]\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "colab_type": "text", + "id": "l9zEl6bG6lg3" + }, + "source": [ + "### Euclidian distance ###" + ] + }, + { + "cell_type": "code", + "metadata": { + "colab_type": "code", + "outputId": "e8351a3e-d446-481d-9bd4-5ef5eb35bff9", + "id": "b52PSKkj6nbm", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 68 + } + }, + "source": [ + "# LN-DD2363-part1.pdf, pp. 24\n", + "x = np.array([-1, 2, 3])\n", + "y = np.array([4, 0, -3])\n", + "\n", + "def euclideandistance(x, y):\n", + " result = 0\n", + " assert(x.size == y.size)\n", + " for i in range(x.size):\n", + " result += (x[i]-y[i])*(x[i]-y[i])\n", + " result = np.sqrt(result)\n", + " return result\n", + "\n", + "print(\"x = \", x)\n", + "print(\"y = \", y)\n", + "print('euclideandistance(x, y) =', euclideandistance(x, y))" + ], + "execution_count": 68, + "outputs": [ + { + "output_type": "stream", + "text": [ + "x = [-1 2 3]\n", + "y = [ 4 0 -3]\n", + "euclideandistance(x, y) = 8.06225774829855\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "SsQLT38gVbn_", + "colab_type": "text" + }, + "source": [ + "# **Results**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "RLwlnOzuV-Cd", + "colab_type": "text" + }, + "source": [ + "The results appear to be correctly implemented and match what is found in e.g the examples at: \n", + "* http://tutorial.math.lamar.edu/Classes/CalcII/DotProduct.aspx\n", + "* https://mathinsight.org/matrix_vector_multiplication\n", + "* http://rosalind.info/glossary/euclidean-distance/" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_4GLBv0zWr7m", + "colab_type": "text" + }, + "source": [ + "# **Discussion**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "6bcsDSoRXHZe", + "colab_type": "text" + }, + "source": [ + "The results were as expected since the operations were rather simple and basic in terms of linear algebra." + ] + } + ] +} \ No newline at end of file diff --git a/Lab-2/lab2_matrixfactorization.ipynb b/Lab-2/lab2_matrixfactorization.ipynb new file mode 100644 index 0000000..d3cb389 --- /dev/null +++ b/Lab-2/lab2_matrixfactorization.ipynb @@ -0,0 +1,566 @@ +{ + "nbformat": 4, + "nbformat_minor": 0, + "metadata": { + "colab": { + "name": "lab2_matrixfactorization.ipynb", + "provenance": [], + "collapsed_sections": [] + }, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + }, + "accelerator": "GPU" + }, + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "6RgtXlfYO_i7", + "colab_type": "text" + }, + "source": [ + "# Lab 2: Matrix factorization #\n", + "**Viktor Meyer - DD2363 Methods in Scientific Computing**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "9x_J5FVuPzbm", + "colab_type": "text" + }, + "source": [ + "# **Abstract**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "6UFTSzW7P8kL", + "colab_type": "text" + }, + "source": [ + "This lab is an exercise in matrix factorization. The mandatory part includes sparse matrix vector product, QR factorization and direct solver (Ax=b). There is also an extra assignment least squares problem, QR eigenvalue algorithm or blocked matrix-matrix product." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "28xLGz8JX3Hh", + "colab_type": "text" + }, + "source": [ + "# **Environment**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "D2PYNusD08Wa", + "colab_type": "text" + }, + "source": [ + "To have access to the neccessary modules you have to run this cell." + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "Xw7VlErAX7NS", + "colab_type": "code", + "colab": {} + }, + "source": [ + "# Load neccessary modules.\n", + "from google.colab import files\n", + "\n", + "import time\n", + "import numpy as np\n", + "\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib import tri\n", + "from matplotlib import axes" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "gnO3lhAigLev", + "colab_type": "text" + }, + "source": [ + "# **Introduction**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "l5zMzgPlRAF6", + "colab_type": "text" + }, + "source": [ + "The methods to be implemented were covered in the during lectures and definitions are available in the lecture notes: LN-DD2363-part3.pdf. \n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "WeFO9QMeUOAu", + "colab_type": "text" + }, + "source": [ + "# **Methods**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "I5Kf3qRyY5J5", + "colab_type": "text" + }, + "source": [ + "### Sparse Matrix-Vector Product ###" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "eRBfBO2AtXdb", + "colab_type": "text" + }, + "source": [ + "Sparse matrices are matrices that contain a lot of zeros and are not dense in information. This is a problem since they essentially waste space. CRS is a datastructure that can be employed to make more efficient use of computer memory. \n", + "\n", + "The smvp (SparseMatrixVectorProduct) method below takes a CRS structure as input and performs matrix vector product calculation:" + ] + }, + { + "cell_type": "code", + "metadata": { + "colab_type": "code", + "id": "PbnT724f5mfQ", + "colab": {} + }, + "source": [ + "# LN-DD2363-part3.pdf, pp. 7\n", + "\n", + "def smvp(val, col_idx, row_ptr, x):\n", + " b = np.matrix(np.ndarray(shape=(len(x), 1)))\n", + " \n", + " for i in range(len(x)):\n", + " b[i] = 0\n", + " for j in range(int(row_ptr[i]), int(row_ptr[i+1])):\n", + " b[i, 0] += val[j] * x[int(col_idx[j])]; \n", + " \n", + " return b" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "fjImmRK-t2xF", + "colab_type": "text" + }, + "source": [ + "The crs method below takes a matrix A as input and converts it to CRS form:" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "scHcB7wBgzWh", + "colab_type": "code", + "colab": {} + }, + "source": [ + "def crs(A):\n", + " val = np.array([])\n", + " col_idx = np.array([])\n", + " row_ptr = np.array([])\n", + "\n", + " for dy in range(A.shape[1]):\n", + " addedRowPtr = False\n", + " for dx in range(A.shape[0]):\n", + " if(A[dy, dx] != 0):\n", + " col_idx = np.append(col_idx, dx)\n", + " if not (addedRowPtr):\n", + " row_ptr = np.append(row_ptr, len(val))\n", + " addedRowPtr = True\n", + " val = np.append(val, A[dy, dx])\n", + "\n", + " row_ptr = np.append(row_ptr, len(val))\n", + "\n", + " return np.array([val, col_idx, row_ptr]) " + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "LLNy5emTuGrz", + "colab_type": "text" + }, + "source": [ + "We test the SVMP method with a simple matrix to ensure that the SVMP is equal to ordinary matrix vector product:" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "75ulld-xuE3X", + "colab_type": "code", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 172 + }, + "outputId": "5fac676e-9775-4213-bcc9-5e6e0c0c1fd6" + }, + "source": [ + "def smvp_test():\n", + " A = np.matrix(\"2 2; 2 2\")\n", + " x = np.matrix(\"2; 2\")\n", + " crsform = crs(A)\n", + " print(\"A = \", A)\n", + " print(\"x = \", x)\n", + " print(\"b = \", smvp(crsform[0], crsform[1], crsform[2], x))\n", + " print(\"b(dense) = \", A*x)\n", + "\n", + "smvp_test()" + ], + "execution_count": 99, + "outputs": [ + { + "output_type": "stream", + "text": [ + "A = [[2 2]\n", + " [2 2]]\n", + "x = [[2]\n", + " [2]]\n", + "b = [[8.]\n", + " [8.]]\n", + "b(dense) = [[8]\n", + " [8]]\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "FikvZiv9uVeC", + "colab_type": "text" + }, + "source": [ + "We also test the SVMP method and compare it to ordinary vector product to see how it behaves for large matrices:" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "q-R-tnP-g3ny", + "colab_type": "code", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 279 + }, + "outputId": "56dd177a-4a5f-4722-cdd4-09ed3fbc7817" + }, + "source": [ + "def smvp_plot():\n", + " sizes = np.array([])\n", + " errors = np.array([])\n", + "\n", + " for n in range(128):\n", + " A = np.random.rand(n, n)\n", + " x = np.random.rand(n, 1)\n", + "\n", + " crsform = crs(A)\n", + "\n", + " sparseresult = smvp(crsform[0], crsform[1], crsform[2], x)\n", + " denseresult = A*x\n", + " \n", + " error = np.abs(np.sum(denseresult-sparseresult))\n", + " \n", + " sizes = np.append(sizes, n*n)\n", + " errors = np.append(errors, error)\n", + "\n", + " plt.figure()\n", + " plt.xlabel('Matrix size')\n", + " plt.ylabel('Error')\n", + " plt.plot(sizes, errors, 'o-')\n", + "\n", + "smvp_plot()" + ], + "execution_count": 100, + "outputs": [ + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEGCAYAAACpXNjrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5zVdb3v8ddnLsBwHVBAGCCwEG8o\n2CSWXbwkkHY25K7U095ZufOxT9bpSmH1yJ1dIN1na53dtnyUO+2mVobstJCQTm0LEQJEFBJBZUbu\nMFyHYS6f88fvu4Y1a9Zas2ZmXWfez8djPea3vr/f+n2/84OZz3zv5u6IiIhkU1mhCyAiIn2PgouI\niGSdgouIiGSdgouIiGSdgouIiGRdRaELUCxOP/10nzx5cqGLISJSUtauXbvP3Ucnpiu4BJMnT2bN\nmjWFLoaISEkxs1eSpatZTEREsk7BRUREsk7BRUREsk7BRUREsk7BRUREsk6jxURE+oEl6+q5c9kW\nXmtoZHx1FQvmTGP+zJqc5afgIiLSxy1ZV8+tj2yksbkVgPqGRm59ZCNAzgKMmsVERPq4O5dtaQ8s\nMY3Nrdy5bEvO8sxpcDGzl81so5mtN7M1IW2UmS03sxfD15Eh3czsO2a21cyeNbOL4u5zY7j+RTO7\nMS79jeH+W8NnLV0eIiL90WsNjd1Kz4Z81Fwud/cZ7l4b3i8EVrj7VGBFeA/wLmBqeN0M3ANRoABu\nA2YBFwO3xQWLe4CPxn1ubhd5iIj0O+Orq7qVng2FaBabB9wfju8H5selP+CRVUC1mY0D5gDL3f2A\nux8ElgNzw7nh7r7Ko+00H0i4V7I8RET6nc/NPgtLSKuqLGfBnGk5yzPXwcWBJ8xsrZndHNLGuvvO\ncLwLGBuOa4AdcZ+tC2np0uuSpKfLowMzu9nM1pjZmr1793b7mxMRKQXnjB+OA9VVlQAY8OVrzinp\n0WJvdfd6MxsDLDezzfEn3d3NzHNZgHR5uPu9wL0AtbW1OS2HiEihPLFpN2bwxGfeTsPxZmbf9Ud2\nHz6R0zxzWnNx9/rwdQ/wa6I+k92hSYvwdU+4vB6YGPfxCSEtXfqEJOmkyUNEpN954vldzJxYzZhh\ngzhr7DDmnDeWH/35ZY6caM5ZnjkLLmY2xMyGxY6B2cBzwFIgNuLrRuDRcLwU+GAYNXYJcCg0bS0D\nZpvZyNCRPxtYFs4dNrNLwiixDybcK1keIiL9St3B4zxXf5jZ553Rnvbxy6dy+EQLP1n1as7yzWWz\n2Fjg12F0cAXwM3f/nZk9AzxsZjcBrwDvD9c/DlwNbAWOAx8GcPcDZvY14Jlw3e3ufiAcfwz4EVAF\n/Da8ABanyENEpF9Z/vxuAObEBZfpE0Yw7Yxh3LlsM3f8bnNOZuznLLi4+zbgwiTp+4Erk6Q7cEuK\ne90H3JckfQ1wfqZ5iIj0NV0t6/LEpt1MHTOUKacP6fCZ7fuO0RZ6mnMxY18z9EVESlRsWZf6hkac\nU0Fiybqo+/ngsZOsfvkAs8/rOGD2zmVbONnS1iEt2zP2FVxEREpUV8u6rNi8h9Y2Z/a5Z3S4Jh8z\n9hVcRERKVFdB4olNuzhj+CCm14zocD4fM/YVXEREStToYQOTpo+oquAti1bwxPO7OXKimaUbXutw\nfsGcaVRVlndIy/aMfS25LyJSQuI78FPN/D7U2EJDYwsAx062duqsj33N5f4uCi4iIiUicV8WgHKD\n4VWVNBxvZnx1FQePn+T4yeT9MPHBY/7MmpJe/kVERLIkWQd+q8PgARWs+8psAKYsfCzpZ3O5vH4y\n6nMRESkRmYzyKsTy+skouIiIlIhMAkc+OuszoeAiIlIiFsyZ1uW+LPNn1rDo2unUVFdhQE11FYuu\nnZ7T/pVk1OciIlIi5p5/Bp9+CIYOrOBYU0vKUV657qzPhIKLiEiJeGHnYRz41/ddyNzzz+jy+kJS\ns5iISIl4tu4QABdMGNHFlYWn4CIiUiKerTvE6UMHMG7EoEIXpUtqFhMRKYCulspPZmN9A9NrRhD2\nySpqqrmIiORZV0vlJ3OsqYWte45ywYTq/BW0FxRcRETyrKul8pPZ9Nph2rw0+ltAwUVEJO96sp/K\ns3UNAJ2Wzy9WCi4iInmWaqa9WbQ22KWLn+zURLax/hBnDB/EmOHF35kPCi4iInm3YM40Kso6d8q3\nOSn7YDbWHSqZJjFQcBERybv5M2t4/eghVJQZBpQnGf0V3wdzqLGZbfuOKbiIiEhq7s7OQyd4X+1E\nti++hjZPvu1XrA9mU300eXJ6iYwUAwUXEZG8277vGIdPtDBjYlQT6Wq142djwaVEOvNBwUVEJO/W\n74hGfs2YOBJIvkx+udG+2vHGukNMHFXFqCED8lvQXtAMfRGRPFu/o4EhA8p5w5ihQOc97YcOrOBI\nU0v7+Q11DVxYQk1ioOAiIpJ3G3Y0MH3CCMrjRozFL5N/+EQz77hjJXcs28Ld182g7mAj/3DJ6wpV\n3B5Rs5iISB6daG7l+Z2H25vEkhk+qJJbLn8Df/zbXt5+x0oAfvCnbWmXhyk2Ci4iInn0/M7DNLc6\nMyamb+YaUVUJwNGmFgD2HT3Z5fpjxUTBRUQkj9a/GuvMTx9c7v79i53Sulp/rJgouIiI5NGGugbO\nGD6IM7rYk6Un648Vk5wHFzMrN7N1Zvab8H6KmT1tZlvN7CEzGxDSB4b3W8P5yXH3uDWkbzGzOXHp\nc0PaVjNbGJeeNA8RkUJbv6Ohy1oLdD33pdjlo+bySeCFuPffAu5y9zcAB4GbQvpNwMGQfle4DjM7\nF7geOA+YC/xHCFjlwHeBdwHnAjeEa9PlISJSMAeOneSV/ce5MIPgkmzuS1Vlefvcl2KX0+BiZhOA\na4AfhPcGXAH8MlxyPzA/HM8L7wnnrwzXzwMedPcmd98ObAUuDq+t7r7N3U8CDwLzushDRKRgNuzI\nrL8FoqHJi66dTk11FQbUVFex6NrpXe5WWSxyPc/lbuDzwLDw/jSgwd1bwvs6IPakaoAdAO7eYmaH\nwvU1wKq4e8Z/ZkdC+qwu8ujAzG4GbgaYNGlSD749EZHMrd/RQJllvuFX/NyXUpOz4GJm7wb2uPta\nM7ssV/n0hrvfC9wLUFtbm3zlOBGRblqyrr59tv346ioWzJnG/Jk1rN/RwFljhzFkYN+fv57L7/BS\n4O/M7GpgEDAc+DZQbWYVoWYxAYgN2q4HJgJ1ZlYBjAD2x6XHxH8mWfr+NHmIiPRKYuC4/OzRrNy8\nt8P7X62tb9/GuL6hkQW/2MBX/2sTB483M3hAOUvW1ZdsjSRTOetzcfdb3X2Cu08m6pB/0t0/AKwE\n3hsuuxF4NBwvDe8J5590dw/p14fRZFOAqcBq4BlgahgZNiDksTR8JlUeIiI9tmRdPbc+spH6hsb2\nTb1+surVTu9jgSWmuc05eLwZgOMnW0tqMmRPFWKeyxeAz5jZVqL+kR+G9B8Cp4X0zwALAdx9E/Aw\n8DzwO+AWd28NtZKPA8uIRqM9HK5Nl4eISI/duWxLp8DRE6U0GbKnzFNsUtPf1NbW+po1awpdDBEp\nYlMWPka2fmMasH3xNVm6W+GY2Vp3r01M1wx9EZEMZTqBsfOmxT2/V6lScBERydCCOdOoKEsfOqoq\ny/nAJZPa56dUV1VSWW6drimVyZA91ffHw4lIv5VqSHBPzZ9Zw13Lt/DaoRO0tHrS0WLJ8sh2OUqB\ngouI9EmxkV3xQ4JvfWQjQI9/se872sSrBxv55JVT+dQ7z8r4c6U8GbKn1CwmIn1SspFdvR2ltXLz\nHtzhneeM7W3x+jwFFxHpk3KxZP2KF/ZwxvBBnDd+eI/v0V8ouIhIn5TtJeubWlr504t7ueKcMUTr\n40o6Ci4i0ictmDOt05Dg3ozSWrXtAMdOtvLOc8b0vnD9gDr0RaRPumjSSBwYNqiCIydaGDKgnG+8\np+OS9d0ZxbXihd0MqizjLa8/PU/fQWlTzUVE+qQ//G0PAP/18bfynpk1lJUZc847o/18snXCUq35\n5e6seGEPb5s6mkEJG3hJcgouItInrdy8h8mnDWby6UO4/k0TOXKihcc27mw/353RZJt3HaG+oVFN\nYt2g4CIifc6J5lb+/NJ+LpsWBYOLp4zizNFDeHD1q+3XZDqabMm6eq77/l8AuGv5i31+NeNsUXAR\nkT7nL9v209TSxuVnR8HFzLj+TRNZ88pBXtx9BICRQwYk/Wz8aLJY09nhE9HGtrsOn+gXy+Vng4KL\niPQ5f9i8h0GVZcyaMqo97e8vmkCZwfz/eIopCx/jwLGTnUaTVZZbh9FkuZiI2V8ouIhIn+LurNyy\nl0tff3qHzvc/vbgPA441tbYvm19uMHJwJUYUWAZVlHHVuadm3+diImZ/oaHIItKnbN93jFcPHOej\nbz+zQ/qdy7bQmrAZS4vD4AEVrPvKbNa+cpC/v+fPXPLNFRxtamHM8IEp8+jry+Vng2ouItKnrNyy\nF4DLzhrdIb2rWsiOA8cpN+NIUwsO7D7chEOnJfb7w3L52aCai4j0CbEJkfUNjVSUGWtfOcjEUYPb\nz4+vrqI+SYCJ1UKimk3nfSaHDqxgyMCKfrVcfjYouIhIUejunifx14+oquTYyRaaQ7tXS5t3Wl5/\nwZxpHZbgh461kFQ1m0ONzay/bXZWvsf+RM1iIlJw3Zktn+z6hsbm9sASkziqa/7MGhZdO719h8ia\n6ioWXXtqOZhsL3TZ36nmIiIFl27Ib7LaS7Lrk0msjaTbtKurmo10j4KLiBRcd4f8ZjoUuDu1jljQ\n6W/bEeeKgouIFFxXne2ZXh+vJ7WO/rgdca6oz0VECm7BnGlUlmc+5HfejPGd0irLrH1CZGJ/iuSf\nai4iUnDzZ9bwizWv8ueXDuCAAbfPOy/l3itmMGxgOUMHVbLr0Ak1YRUhBRcRKQr7jjbztrNGc9Nb\np3Djfas5beiphSVjo8Nine3ucLLV+cLcsxVQipSaxUSk4A4eO8mW3UeYNWUUl5w5iqEDK1j+/O72\n88lGhzW1tGkBySKm4CIiBffMyweAaN+VgRXlvOOs0fz+hT20tUVzV7SAZOlRcBGRgnvm5QMMqCjj\nggkjALjq3LHsPdLE+roGQBMcS5GCi4jkxZJ19Vy6+EmmLHyMSxc/2WH2/ertB5gxsZqBFdES+ZdP\nG0N5mbU3jV1/8cRO99MEx+KWs+BiZoPMbLWZbTCzTWb21ZA+xcyeNrOtZvaQmQ0I6QPD+63h/OS4\ne90a0reY2Zy49LkhbauZLYxLT5qHiBRGuuVdjja18Nxrhzts7DVicCWzpoxi+fO7aWtzVm7ew5AB\n5YwbMUhDjUtELkeLNQFXuPtRM6sE/tvMfgt8BrjL3R80s+8BNwH3hK8H3f0NZnY98C3gOjM7F7ge\nOA8YD/zezM4KeXwXuAqoA54xs6Xu/nz4bLI8RKQA0i3vMmrIAFrbnIvjggvAGcMH8ueX9nPmFx8H\notrL4msvyFuZpXdyVnPxyNHwtjK8HLgC+GVIvx+YH47nhfeE81eamYX0B929yd23A1uBi8Nrq7tv\nc/eTwIPAvPCZVHmISAGk65B/5uUDlJcZF00a2Z6+ZF09j23c1eHaR9fVa+/6EpLTPhczKzez9cAe\nYDnwEtDg7i3hkjogVq+tAXYAhPOHgNPi0xM+kyr9tDR5JJbvZjNbY2Zr9u7d25tvVUTSSNch//T2\nA5w/fjhDBp5qSLlz2RaaWto6XNvYrKHHpSSnwcXdW919BjCBqKZxdi7z6y53v9fda929dvTo0V1/\nQER6ZMGcaQxIWN4F4FhTC6u3H+ClvUc71Eo09Lj05WW0mLs3ACuBNwPVZhb7E2UCEPsfVQ9MBAjn\nRwD749MTPpMqfX+aPESkAObPrGHmpGqMaGmXIQOiUWENjc0AHG1q7bB/i4Yel75cjhYbbWbV4biK\nqOP9BaIg895w2Y3Ao+F4aXhPOP+ku3tIvz6MJpsCTAVWA88AU8PIsAFEnf5Lw2dS5SEiBdDW5ry0\n9zjXXDCO7YuvoXpw5wGc8Zt7LZgzjarK8g7nNfS4tORytNg44H4zKycKYg+7+2/M7HngQTP7OrAO\n+GG4/ofAj81sK3CAKFjg7pvM7GHgeaAFuMXdWwHM7OPAMqAcuM/dN4V7fSFFHiJSABvrD7HvaBNX\nnjMG6LrZS3urlL6cBRd3fxaYmSR9G1H/S2L6CeB9Ke71DeAbSdIfBx7PNA8RKYwVL+ymzOAdZ0XB\nJZP9W7S3SmnTDH0RybkVm/dw0aSRjBoSNYep2avvU3ARkZzadegEm147zJXnjG1Pmz+zhkXXTqem\nukoz7vuoLpvFQp/Jt9z9c3koj4j0MU9u3gPQ3t8So2avvq3L4OLurWb21nwURkT6jtjOkfUNjZSX\nGZvqD3HW2GGFLpbkSaYd+uvMbCnwC+BYLNHdH8lJqUSkpCXuHNna5nzx189hZqqt9BOZBpdBRJMT\nr4hLc0DBRaQPi9+3vjvDgdMtVKng0j9kFFzc/cO5LoiIFJfE2kdsmXygywCh5Vsko9FiZjbBzH5t\nZnvC61dmNiHXhRORwklX++jKuOpBSdO1fEv/kWmz2H8CP+PUJMd/CGlX5aJQIlJ43a19xDehVSRZ\npFLzWPqXTOe5jHb3/3T3lvD6EaBlhEX6sO4sHpm402Rzq2PAyMGVmsfST2Vac9lvZv8A/Dy8v4Go\ng19E+qgFc6bx+V9u4GSrt6cNqixLWvtI1oTmwOABFaz7yuxcF1WKUKY1l48A7wd2ATuJVhxWJ79I\nHzZ/Zg0XTxnVvkw+wOxzx3aofSxZV8+li59Muk4YqAO/P8t0hv617v53eSiPiBSRPUeaeNtZo3ng\nIxfzjz98mv/eup9jTS0MGVjRaTRZMurA778ynaF/A3BXHsojIkVi39Em/rb7aHtN5VPvPIu/v+fP\nvHnRCo6caKHMIK7FrBN14Pdvmfa5PGVm/w48RMcZ+n/NSalEpOBWbz8AwCVnngbAjgPHKTM4fKIF\nSB9YarT/Sr+XaXCZEb7eHpfmdJyxLyJ9yKpt+xk8oJzpNSOAqNO+LU1AiampruKphfrV0N9l0udS\nBtzj7g/noTwiUiRWbdtP7eRRVJZH434y6ZxXU5jEdDlazN3bgM/noSwiUiRi/S2XnDmqPS1V53y5\nmeaySCeZNov93sw+R+c+lwM5KZWIFFSsv2XWlNPa0xbMmdZpdFhVZbkCiiSVaXC5Lny9JS7NgTOz\nWxwRKQartu2nqrKcCyaMaE+LBZCerJIs/U+mqyJPyXVBRKTw4jf4GlhRxmPP7uwQPLR7pGQqbZ+L\nmX0+7vh9Cee+matCiUj+xa8PBtDU0satj2xkybr6ApdMSlFXHfrXxx3fmnBubpbLIiI5EFuiZcrC\nx7h08ZOdgkXs/KceWt/jJfZFEnXVLGYpjpO9F5Ei09WGX5ks4aL1waQnugounuI42XsRKTKpNvz6\n7MMb+PRD67tcwgW0Ppj0TFfB5UIzO0xUS6kKx4T3ybeaE5GikarW0eoevqb/vCZFSk+lDS7uXp6v\ngohI9o2vrkq5HH5XtD6Y9Eam+7mISAlaMGcaVZXd+zGvqizn7utm8NTCKxRYpMcUXET6sPkza/j8\n3LPb35db8nE4WsJFsi3TGfoiUqLOGB51jz56y6Vs33dMS7hIXuSs5mJmE81spZk9b2abzOyTIX2U\nmS03sxfD15Eh3czsO2a21cyeNbOL4u51Y7j+RTO7MS79jWa2MXzmO2bRn2Wp8hDpjzbUHWJAeRln\njxvG/Jk1LLp2OjXVVaqpSE7lsubSAnzW3f9qZsOAtWa2HPgQsMLdF5vZQmAh8AXgXcDU8JoF3APM\nMrNRwG1ALdHw57VmttTdD4ZrPgo8DTxONLHzt+GeyfIQ6Xc27GjgnPHDGVgRjc/REi6SDzmrubj7\nzthOle5+BHgBqAHmAfeHy+4H5ofjecADHlkFVJvZOGAOsNzdD4SAshyYG84Nd/dV7u7AAwn3SpaH\nSL/S2uZsrD/EhXELUIrkQ176XMxsMjCTqIYx1t13hlO7gLHhuAbYEfexupCWLr0uSTpp8kgs183A\nzQCTJk3q5nclUjxiC04mrla8be9Rjja1cOGE6kIXUfqZnAcXMxsK/Ar4lLsftrjRKu7uZpbTmf7p\n8nD3e4F7AWpra7XigJSkdEu8NLe2AXDhRAUXya+cDkU2s0qiwPJTd38kJO8OTVqEr3tCej0wMe7j\nE0JauvQJSdLT5SFSslItQJlqiZc7l21hQ10DwwZWcObpQwpRZOnHclZzCSO3fgi84O7/FndqKXAj\nsDh8fTQu/eNm9iBRh/4hd99pZsuAb8aN+JoN3OruB8zssJldQtTc9kHg/3aRh0hJid9fxTi1oF99\nQyMLfrGBr/7XJg4eb0762dcaGtmw4xAXTBxBWZnWmZX8ymWz2KXAPwIbzWx9SPsi0S/8h83sJuAV\n4P3h3OPA1cBW4DjwYYi2UjazrwHPhOtuj9te+WPAj4AqolFivw3pqfIQKRmJzV2J7bbNbZ4ysACM\nGzGIF3Ye5ua3a8NYyb+cBRd3/29SL8t/ZZLrnY7bKMefuw+4L0n6GuD8JOn7k+UhUkqSNXdlqtyM\n979pInf//kX1t0hBaPkXkSLV031Uhg6soNWde/+4DYDbHt2k3SQl7xRcRIpUT/ZRqamu4ovXRGuJ\nHT8Z1Xp2HT6h7Yol7xRcRIrUgjnTSFxnMva2uqqSyvKOJ2N7r3z3yZc63UvbFUu+aeFKkSJ1yZmn\n4Q7DB1Vw5ERLh8mRkHri5KcfWp/0ftquWPJJwUWkSC3btAuARz52KW8YM7TT+VRrhKXaIEzbFUs+\nqVlMJEdSTXrM1O+e28UbxgxNGljSiTYI67iJrLYrlnxTzUUkB9ItyTJ/Zk3KJq2Y/UebeHr7fm65\n/A3dzjt2n3T3F8k1i6aXSG1tra9Zs6bQxZASFz+jPplyM1rdO8y2B9rf11RXcfnZo/nNhp00NDYz\nethAvnT1OQoMUrTMbK271yamq+YikiWJtZVkWsMfc4l/0sUv6/KTVa+2p+890tShxiNSKtTnIpIl\nvZlRn46GEUspUnARyZJcDvXVMGIpNWoWE+mBZB3yqYYAw6m+lp7SMGIpNaq5iHRTrG+lvqER59RI\nsMvPHs3Aio4/UlWV5dx93Qz+z/sv7DQ82BK+pqJhxFKKVHMRSSNZDSXV5lwrN+/l2otq+PnqHRgk\nHQKcbHhwYh6Xnz2alZv3ahixlDQNRQ40FFkSJRv9VVVZnrLT3oD3vnECy1/YzV+/fJU26JJ+IdVQ\nZDWLiaSQqoZSnriaZDC+uorVLx/gTZNHKbBIv6fgIpJCqhFayTrmqyrL+ejbp/DK/uPMmjIq10UT\nKXoKLiIpdDVCa9igqMuystxYdO10Rg0ZCMCsKaflvGwixU7BRSSFBXOmpR3J1dLqzJ8xnuZW5/ya\nEazevp+hAys4Z9ywvJVRpFgpuIikcOHEahwYUVWZ9Hxjcyurth2gzODR9fWs3n6AN75uJBXl+rES\n0U+BSAorXtgNwGP/+60pazC7D59g6pihfHflVv62+yjrXj2o7YRF0DwXESD5fJYVL+zh7DOGMWHk\n4JSz70dUVbJt3zHaQh//4RMtWmhSBNVcRJLOuF/4yLOs2rafK84eA6TegMsMmls7jh7TQpMiCi4i\nSeeznGhuw4GH19SxZF0982fWsOja6dRUV2FE+64sunY6Dcebk95TC01Kf6dmMen30gWCfUc77qeS\n2NSVamMwLTQp/Z1qLtLvjRk+MO35dM1c2q9eJDnVXKTfinXi7z7c1OW1qWo32q9eJDkFF+mXki1K\nmbivfbx0zVzJmstE+jsFF+nz4ocZj6iqxAwOJumId6C6qpKmlrZOKyGrmUuke9TnIn1a4jDjhsbm\npIEl5lBjc9JRYaqZiHRPzmouZnYf8G5gj7ufH9JGAQ8Bk4GXgfe7+0EzM+DbwNXAceBD7v7X8Jkb\ngS+H237d3e8P6W8EfgRUAY8Dn3R3T5VHrr5PKS6xWkp9Q2OPthYeX12lZi6RLMhlzeVHwNyEtIXA\nCnefCqwI7wHeBUwNr5uBe6A9GN0GzAIuBm4zs5HhM/cAH4373Nwu8pA+Lr6WAsmXxk9HzV8i2ZOz\n4OLufwQOJCTPA+4Px/cD8+PSH/DIKqDazMYBc4Dl7n4g1D6WA3PDueHuvsqjrTQfSLhXsjykD1iy\nrp5LFz/JlIWPceniJzus45VsMmSm1Pwlkl357tAf6+47w/EuYGw4rgF2xF1XF9LSpdclSU+XRydm\ndjNRTYlJkyZ193uRPEsc4VXf0NhhgmNPZsVXVZYrqIjkQME69EONo3vtFlnOw93vdfdad68dPXp0\nLosiWZBq2+HYBMdxIwZ1eY/qqkpGDq5UZ71IjuW75rLbzMa5+87QtLUnpNcDE+OumxDS6oHLEtL/\nENInJLk+XR5S4lLVTGLps6aM4tfrX0t6jWooIvmV75rLUuDGcHwj8Ghc+gctcglwKDRtLQNmm9nI\n0JE/G1gWzh02s0vCSLMPJtwrWR5SAtL1qaSbyDjj9if49frXKDcYOTja3Kvcol1YVEMRyb9cDkX+\nOVGt43QzqyMa9bUYeNjMbgJeAd4fLn+caBjyVqKhyB8GcPcDZvY14Jlw3e3uHhsk8DFODUX+bXiR\nJg8pcl31qdz4ltfxzcc3d/qcQ/vqxK0erWh893UzFExECihnwcXdb0hx6sok1zpwS4r73AfclyR9\nDXB+kvT9yfKQ4peuT2X+zBqOnGgBoMxo35wrmfjPiEhhaPkXKRqp+lTqGxp5y+IVvNZwgoEVZTS1\ntPX4XiKSH1r+RYpGqj4VA15rOAFAU0tbyv3sM7mXiOSHgosUjX++7MxOaclWKvaQnopm2osUnprF\npODi1wMDGD6ogsMnWqiqLE85496JRoHFr3TccLxZ+6mIFAkFF8mp+OXuk/3iT7avSnOr85bXj2Ld\nq4cYN2IQOw+d6HTfmuoqnlp4RV6+BxHpPgUXyZmuhhYvWVfPZx/e0GmBycbmVp6tO0RjcyuNhzrX\nXNTsJVL8zLu5cmxfVVtb62vWrCl0MfqExGauRLGl8NPt/JhKjZq9RIqKma1199rEdNVcJKuSNXMl\nitVUehJY1BQmUho0WkyyqoBaxUIAAAuMSURBVDfL3ndFc1dESoeCi2TNknX1KZvCMlVu1r42WCLN\nXREpHWoWkx6LHwk2oqqSYydbenW/2MrFQKemNXXii5QWBRfplvjO+vgO+YbG5pSfSddxHzuXrKM+\n3RBmESluCi6SkSXr6vmXpZs6BJFMO+Q/cMkkfrW2vlNfzMjBldz2P85LGjTmz6xRMBEpYQou0qVM\nRoClUlNdxdfnT6f2daNUExHpRxRcpEs9HQEW30+imohI/6LgIh3E96nEJjv2lHZ/FOm/FFykXWLz\nV28CS011lQKLSD+m4NLPxQ8nLutmTWXk4EquuWBcp856DRsWEQWXfirZ6K/uBJb4PerVWS8iiRRc\n+plkQaW7Epu81FkvIokUXPqBrlYp7g41eYlIJhRc+rBs1FLg1BL5Wu5eRDKl4NJHJK7zdbKllePN\nbb26p5a4F5GeUnApMYlBxAwOHm/OeJ2vTKn5S0R6Q8GlRCRr4urJOl/JxIYUr9y8VyO+RCQrFFyK\nVLLl7Jtbs7sltREtKvn1+dOzel8REQWXIvTlJRv56apXs9rMFVNm0Obai15EckvBpUBy0QGfTrrl\n7UVEsk3BJY9SDQ3OZs0kkZq+RKQQFFyyLFtzS3qqsswYOqiChuPN6pgXkYLps8HFzOYC3wbKgR+4\n++Js5zHrG8vZfeRktm+bEQtjj2PDkRVMRKSY9MngYmblwHeBq4A64BkzW+ruz2crj0IGlqrKcu2V\nIiJFrazQBciRi4Gt7r7N3U8CDwLzsplBvgNLuRkQjfJSYBGRYtcnay5ADbAj7n0dMCvxIjO7GbgZ\nYNKkSfkpWQZis+01XFhESlVfDS4Zcfd7gXsBamtrsztDsQc0XFhE+oq+GlzqgYlx7yeEtKwZO2xA\nj5vGFEREpK/rq8HlGWCqmU0hCirXA/8zmxk8/aWruuzUHzKgnG+8R/0jItL/9Mng4u4tZvZxYBnR\nUOT73H1TtvN5+ktXZfuWIiJ9Qp8MLgDu/jjweKHLISLSH/XVocgiIlJACi4iIpJ1Ci4iIpJ1Ci4i\nIpJ15l7wuYNFwcz2Aq/08OOnA/uyWJxsUbkyV4xlApWrO4qxTND3y/U6dx+dmKjgkgVmtsbdawtd\njkQqV+aKsUygcnVHMZYJ+m+51CwmIiJZp+AiIiJZp+CSHfcWugApqFyZK8YygcrVHcVYJuin5VKf\ni4iIZJ1qLiIiknUKLiIiknUKLr1kZnPNbIuZbTWzhTnOa6KZrTSz581sk5l9MqSPMrPlZvZi+Doy\npJuZfSeU7VkzuyjuXjeG6180sxuzULZyM1tnZr8J76eY2dMh74fMbEBIHxjebw3nJ8fd49aQvsXM\n5mShTNVm9ksz22xmL5jZm4vkWX06/Ps9Z2Y/N7NBhXheZnafme0xs+fi0rL2fMzsjWa2MXzmO2Zh\nr+6elevO8O/4rJn92syqu3oOqX42Uz3rnpQr7txnzczN7PR8Pq9UZTKzT4TntcnM7sj3swLA3fXq\n4YtoOf+XgDOBAcAG4Nwc5jcOuCgcDwP+BpwL3AEsDOkLgW+F46uB3xLtnHwJ8HRIHwVsC19HhuOR\nvSzbZ4CfAb8J7x8Grg/H3wP+Vzj+GPC9cHw98FA4Pjc8v4HAlPBcy3tZpvuBfwrHA4DqQj8roi24\ntwNVcc/pQ4V4XsDbgYuA5+LSsvZ8gNXhWguffVcvyjUbqAjH34orV9LnQJqfzVTPuiflCukTibb3\neAU4PZ/PK8Wzuhz4PTAwvB+T72fl7gouvXkBbwaWxb2/Fbg1j/k/ClwFbAHGhbRxwJZw/H3ghrjr\nt4TzNwDfj0vvcF0PyjEBWAFcAfwm/HDsi/tl0P6cwg/hm8NxRbjOEp9d/HU9LNMIol/ilpBe6GdV\nA+wIv1wqwvOaU6jnBUxO+MWUlecTzm2OS+9wXXfLlXDuPcBPw3HS50CKn810/zd7Wi7gl8CFwMuc\nCi55e15J/g0fBt6Z5Lq8Pis1i/VO7BdFTF1Iy7nQPDITeBoY6+47w6ldwNguypftct8NfB5oC+9P\nAxrcvSXJ/dvzDucPheuzXaYpwF7gPy1qrvuBmQ2hwM/K3euBfwVeBXYSff9rKfzzisnW86kJx9ku\nH8BHiP6y70m50v3f7DYzmwfUu/uGhFOFfF5nAW8LzVn/z8ze1MMy9epZKbiUIDMbCvwK+JS7H44/\n59GfGHkbX25m7wb2uPvafOWZoQqi5oJ73H0mcIyomaddvp8VQOjDmEcU/MYDQ4C5+SxDpgrxfLpi\nZl8CWoCfFkFZBgNfBL5S6LIkqCCqGV8CLAAezrS/K5sUXHqnnqi9NWZCSMsZM6skCiw/dfdHQvJu\nMxsXzo8D9nRRvmyW+1Lg78zsZeBBoqaxbwPVZhbb6TT+/u15h/MjgP1ZLhNEf2XVufvT4f0viYJN\nIZ8VwDuB7e6+192bgUeInmGhn1dMtp5PfTjOWvnM7EPAu4EPhMDXk3LtJ/Wz7q7XE/2RsCH8/58A\n/NXMzuhBubL5vOqARzyymqhF4fQelKl3z6q7bbR6dWjDrCDqkJvCqY6w83KYnwEPAHcnpN9Jx07Y\nO8LxNXTsVFwd0kcR9UeMDK/twKgslO8yTnXo/4KOHYEfC8e30LGD+uFwfB4dOxu30fsO/T8B08Lx\nv4TnVNBnBcwCNgGDQ173A58o1POic3t91p4PnTuor+5FueYCzwOjE65L+hxI87OZ6ln3pFwJ517m\nVJ9L3p5Xkmf1z8Dt4fgsoiYvy/uz6s0Pr17to0L+RjTa4ks5zuutRM0UzwLrw+tqorbRFcCLRKNE\nYv9ZDfhuKNtGoDbuXh8BtobXh7NUvss4FVzODD8sW8N/0NjIlUHh/dZw/sy4z38plHULGY4s6qI8\nM4A14XktCT/MBX9WwFeBzcBzwI/DD3venxfwc6J+n2aiv3ZvyubzAWrD9/gS8O8kDK7oZrm2Ev2S\njP2//15Xz4EUP5upnnVPypVw/mVOBZe8PK8Uz2oA8JNwr78CV+T7Wbm7ln8REZHsU5+LiIhknYKL\niIhknYKLiIhknYKLiIhknYKLiIhknYKLSA+FVXB/Eve+wsz2WlgZOs3nZpjZ1WnO15rZd3pZtvFm\n9sve3EOkNxRcRHruGHC+mVWF91eR2QzmGUTzCjoxswp3X+Pu/7s3BXP319z9vb25h0hvKLiI9M7j\nRLOxIVrJ9uexE2Z2sZn9JSyc+Wczmxb2w7gduM7M1pvZdWb2L2b2YzN7CvixmV1mp/bF+baZfSUc\nzzGzP5pZh59bM3tHuNf6kNcwM5sc2+MjLNoZO7/XzG4L6QvM7Jmw38hXc/2gpH9RcBHpnQeB681s\nEHAB0SrVMZuBt3m0cOZXgG+6+8lw/JC7z3D3h8K15xItk35Dwv1vJQpElwPfIZrR3ZZwzeeAW9x9\nBvA2oDH+pLv/Uzg3j2gJ9R+Z2WxgKnAxUU3qjWb29p4/BpGOKrq+RERScfdnw/YHNxDVYuKNAO43\ns6lEy/ZUprnVUndvTEx09+Nm9lHgj8Cn3f2lJJ99Cvg3M/sp0YKFdYmL4Ibg9wvgE+7+ipl9gmgD\nrnXhkqFEweaP6b5fkUwpuIj03lKiPVouI1qbK+ZrwEp3f08IQH9Ic49jac5NJ1qhdnyyk+6+2Mwe\nI+rHeSpsX3si4bLvEQWe34f3Bixy9++nyVekx9QsJtJ79wFfdfeNCekjONXB/6G49CNE21R3ycxe\nB3yWaGO4d5nZrCTXvN7dN7r7t4BngLMTzt8CDHP3xXHJy4CPhL2BMLMaMxuTSZlEMqHgItJL7l7n\n7smGDt8BLDKzdXRsJVgJnBvr0E9137DB0w+Bz7n7a0Qr3v4gNHHF+5SZPWdmzxKtjvvbhPOfA6bH\nder/s7s/AfwM+IuZbSTa7yajgCeSCa2KLCIiWaeai4iIZJ2Ci4iIZJ2Ci4iIZJ2Ci4iIZJ2Ci4iI\nZJ2Ci4iIZJ2Ci4iIZN3/B+h8AdMUDC9QAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "tags": [] + } + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "R9-BvtI_uS0Z", + "colab_type": "text" + }, + "source": [ + "It is apparent that there is a precision error when matrix sizes increase, this could be caused by e.g repeated floating point arithmetic error." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "colab_type": "text", + "id": "k9b0yn9qF_Jq" + }, + "source": [ + "### QR Factorization ###" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "R9We443pu31o", + "colab_type": "text" + }, + "source": [ + "QR Factorization is of use when we are solving systems of equations. Triangular matrices are often important when solving equation systems, QR Factorization aims to take a matrix A and express it as triangular matrix R and matrix Q" + ] + }, + { + "cell_type": "code", + "metadata": { + "colab_type": "code", + "id": "401b1-yLGFaH", + "colab": {} + }, + "source": [ + "# LN-DD2363-part1.pdf, pp. 22\n", + "\n", + "def qrf(A):\n", + "\n", + " assert(np.linalg.det(A) != 0)\n", + "\n", + " Q = np.zeros(A.shape)\n", + " R = np.zeros(A.shape)\n", + "\n", + " V = A[:,[0]]\n", + "\n", + " for i in range(A.shape[0]):\n", + " R[i, i] = np.linalg.norm(V)\n", + " Q[:, [i]] = V / R[i, i]\n", + " for j in range(i+1, A.shape[0]):\n", + " R[i, j] = np.sum(np.dot(Q[:,[i]].T, A[:,[j]]))\n", + " V = A[:,[j]] - R[i, j]*Q[:,[i]]\n", + "\n", + " return [Q, R]" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "fPft7QXJvvST", + "colab_type": "text" + }, + "source": [ + "We formulate a quick test to make sure that QR Factorization is working properly, an important property is that R is triangular:" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "o9pt0xjBsIzW", + "colab_type": "code", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 153 + }, + "outputId": "bd331f37-094b-4870-a3ea-76b264f53309" + }, + "source": [ + "def is_triangular(A):\n", + " for i in range(A.shape[0]):\n", + " for j in range(0, i):\n", + " if (A[i, j] != 0): return False\n", + " return True\n", + "\n", + "def qrf_test():\n", + " A = np.matrix(\"2 -1; -1 2\")\n", + " Q, R = qrf(A)\n", + " print(\"A = \", A)\n", + " print(\"Q = \", Q)\n", + " print(\"R = \", R)\n", + " print(\"is_triangular(R) = \", is_triangular(R))\n", + "\n", + "qrf_test()\n" + ], + "execution_count": 102, + "outputs": [ + { + "output_type": "stream", + "text": [ + "A = [[ 2 -1]\n", + " [-1 2]]\n", + "Q = [[ 0.89442719 0.4472136 ]\n", + " [-0.4472136 0.89442719]]\n", + "R = [[ 2.23606798 -1.78885438]\n", + " [ 0. 1.34164079]]\n", + "is_triangular(R) = True\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "colab_type": "text", + "id": "lVE-mL2SDDZu" + }, + "source": [ + "### Direct Solver ###" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "I64B1TRbwBJQ", + "colab_type": "text" + }, + "source": [ + "From litterature we know that if Ax=b then x=A⁻1b, this is what is refered to as a direct solver technique:" + ] + }, + { + "cell_type": "code", + "metadata": { + "colab_type": "code", + "id": "BBZekOW0DGxN", + "colab": {} + }, + "source": [ + "def ds(A, b):\n", + "\n", + " A_inverse = np.linalg.inv(A)\n", + "\n", + " x = A_inverse*b\n", + "\n", + " return x" + ], + "execution_count": 0, + "outputs": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "KOh_k1pQwQsG", + "colab_type": "text" + }, + "source": [ + "We formulate a quick test to make sure that the direct solver is working correctly, for example, if a solution is correct then Ax-b == 0" + ] + }, + { + "cell_type": "code", + "metadata": { + "id": "ASb9__vvo4uR", + "colab_type": "code", + "colab": { + "base_uri": "https://localhost:8080/", + "height": 153 + }, + "outputId": "99018ee5-e1e2-465c-8b2d-191c4e5abe78" + }, + "source": [ + "def ds_test():\n", + " A = np.matrix(\"1 0; 2 2\")\n", + " b = np.matrix(\"1; 6\")\n", + " x = ds(A, b)\n", + " print(\"A = \", A)\n", + " print(\"b = \", b)\n", + " print(\"x = \", x)\n", + "\n", + " print(\"||Ax-b|| = \", np.linalg.norm(A*x-b))\n", + "\n", + "ds_test()" + ], + "execution_count": 104, + "outputs": [ + { + "output_type": "stream", + "text": [ + "A = [[1 0]\n", + " [2 2]]\n", + "b = [[1]\n", + " [6]]\n", + "x = [[1.]\n", + " [2.]]\n", + "||Ax-b|| = 0.0\n" + ], + "name": "stdout" + } + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "SsQLT38gVbn_", + "colab_type": "text" + }, + "source": [ + "# **Results**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "RLwlnOzuV-Cd", + "colab_type": "text" + }, + "source": [ + "The results appear to be correctly implemented. This is not surprising since the various algorithms follow what is described in LN-DD2363-part3.pdf." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_4GLBv0zWr7m", + "colab_type": "text" + }, + "source": [ + "# **Discussion**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "6bcsDSoRXHZe", + "colab_type": "text" + }, + "source": [ + "Sparse Matrix Vector Product was of particular interest. The problem of large datasets is of real world importance. The comparison between dense and sparse matrix products showed that there may be accuracy loss when matrices approach larger sizes." + ] + } + ] +} diff --git a/Lab-3/lab3_iterativemethods.ipynb b/Lab-3/lab3_iterativemethods.ipynb new file mode 100644 index 0000000..8b416bf --- /dev/null +++ b/Lab-3/lab3_iterativemethods.ipynb @@ -0,0 +1,543 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "[](https://colab.research.google.com/github/johanhoffman/DD2363-VT20/blob/tree/viktorme/Lab-3/lab3_iterativemethods.ipynb)\n", + "# Lab 3: Iterative methods #\n", + "**Viktor Meyer - DD2363 Methods in Scientific Computing**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Abstract**\n", + "This lab is an exercise in iterative methods. The mandatory part includes Jacobi iteration, Gauss-Seidel iteration, Newton's method for scalar nonlinear equation f(x)=0. There is also an extra assignment GMRES method or Newton's method for vector nonlinear equation f(x)=0." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Environment**\n", + "To have access to the neccessary modules you have to run this cell." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 1, + "outputs": [], + "source": [ + "# Load neccessary modules.\n", + "#from google.colab import files\n", + "\n", + "import time\n", + "import numpy as np\n", + "from scipy import optimize\n", + "from scipy.misc import derivative\n", + "\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib import tri\n", + "from matplotlib import axes" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Introduction**\n", + "The methods to be implemented were covered in the during lectures and definitions are available in the lecture notes. The aim of this lab is to approach the problem of solving systems of linear equations using iterative methods. The iterative methods are generally fast with a low memory footprint while generally yielding acceptable approximate solutions to large systems. \n", + "> LN-DD2363-part4.pdf, p. 127" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Methods**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "### Jacobi iteration, Ax=b ###\n", + "\n", + "Jacobi iteration is a fixed point iteration method using matrix splitting. The method typically starts with an initial guess where all components of x are zero. In each iteration, the method updates x with new changes that aim to closer approximate a correct solution.\n", + "> LN-DD2363-part4.pdf, pp. 130-133\n", + "> https://www.maa.org/press/periodicals/loci/joma/iterative-methods-for-solving-iaxi-ibi-jacobis-method" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 2, + "outputs": [], + "source": [ + "def jacobi(A, b, iterations = 1000):\n", + " \n", + " #Diagonal from A, apply reciprocal, this is equal to (D-1)\n", + " di = np.diag(np.array(np.reciprocal(A.diagonal().copy().astype(float))).ravel())\n", + " \n", + " #Clear the diagonal we extracted from A, this is equal to (L+U)\n", + " np.fill_diagonal(A, 0)\n", + " \n", + " #Invert each element in A\n", + " A *= -1\n", + " \n", + " solution = np.matrix(np.zeros(b.shape))\n", + " for i in range(iterations):\n", + " solution = di*(A*solution+b)\n", + " \n", + " return solution" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Gauss-Seidel iteration, Ax=b ###\n", + "Gauss-Seidel iteration shares a lot in common with Jacobi iteration, the main difference is in how changes to x are made. Recall that Jacobi iteration only applies updates to all components of x AFTER each iteration. Gauss-Sediel improves upon this by updating each component of x AS SOON AS POSSIBLE, likely leading to faster convergence rates.\n", + "> LN-DD2363-part4.pdf, pp. 130-133\n", + "> https://www.maa.org/press/periodicals/loci/joma/iterative-methods-for-solving-iaxi-ibi-gauss-seidel-method" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 3, + "outputs": [], + "source": [ + "def gaussseidel(A, b, iterations = 1000):\n", + " \n", + " LD = np.tril(A) \n", + " \n", + " U = A-LD\n", + "\n", + " x = np.matrix(np.zeros(b.shape))\n", + " for i in range(iterations):\n", + " x = np.linalg.inv(LD)*(-1*U*x+b)\n", + " \n", + " return x" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Newton's method for scalar nonlinear equations, f(x)=0###\n", + "\n", + "Newton's method is an approach to solving nonlinear scalar equations. The general idea is to use an initial guess and progressively improve over time. Improvements are made by using the derivative as guide in each iteration.\n", + "\n", + "> LN-DD2363-part4.pdf, pp. 147-149" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 4, + "outputs": [], + "source": [ + "def newtonsmethod(f, x0, iterations):\n", + " x = x0\n", + " for i in range(iterations):\n", + " df = derivative(f, x)\n", + " x -= (f(x)/df) \n", + " return x" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Results**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "### Jacobi iteration, Ax=b ###\n", + "\n", + "Below are tests for Jacobi iteration. Input A and b are given such that there exists an exact solution y for x under the condition of Ax=b. This means that the test explores convergence ||Ax-b|| and the approximation error ||x-y|| at the same time." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 5, + "outputs": [ + { + "name": "stdout", + "text": [ + "Using:\n", + "A=\n", + "[[ 4 -1 -1]\n", + " [-2 6 1]\n", + " [-1 1 7]]\n", + "b=\n", + "[[ 3]\n", + " [ 9]\n", + " [-6]]\n", + "Converged in 18:\n", + "x=\n", + "[[ 0.99999999]\n", + " [ 2.00000001]\n", + " [-0.99999999]]\n", + "|error|=\n", + "9.877244466771629e-09\n" + ], + "output_type": "stream" + }, + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5xcdX3/8dd7L0k2F7ILWSCZXRK0iKJAElbUghW1TQI/hYiKibal1T7ys5Xay6/4CO2j6ANbpaW1VqXa/ChS+7PcL6YVjChYLBTNkoSEWyRGkN0EsuQKySbZy+f3x5xdJpuZ3UkyZ2d25/18POYxZ77nnO98mAzz2XO+N0UEZmZmQ9WUOwAzM6tMThBmZpaXE4SZmeXlBGFmZnk5QZiZWV515Q6glGbMmBFz5swpdxhmZmPGY4899nJENOfbN64SxJw5c2hvby93GGZmY4ak5wvt8y0mMzPLywnCzMzycoIwM7O8nCDMzCwvJwgzM8srtQQhqVXSg5KekvSkpD/Kc4wkfUXSJknrJc3P2Xe5pGeTx+VpxXnP2k7Ou/YBTl3+Xc679gHuWduZ1luZmY0paXZz7QX+T0SskTQNeEzS/RHxVM4xFwKnJY+3AV8H3ibpeOCzQBsQybkrI2JnKQO8Z20nV921ge6ePgA6d3Vz1V0bAFg8L1PKtzIzG3NSu4KIiK0RsSbZfgV4Ghj6q3sJ8K3IehRolDQTWAjcHxE7kqRwP7Co1DFet2rjYHIY0N3Tx3WrNpb6rczMxpxRaYOQNAeYB/xkyK4M8ELO646krFB5vrqXSWqX1N7V1XVEcW3Z1X1E5WZm1ST1BCFpKnAn8McRsafU9UfEiohoi4i25ua8o8ULmtXYcETlZmbVJNUEIamebHL4dkTcleeQTqA153VLUlaovKSuXHg6DfW1h5Q11Ndy5cLTS/1WZmZjTpq9mAT8C/B0RHypwGErgd9OejO9HdgdEVuBVcACSU2SmoAFSVlJLZ6X4YuXnsmUidkkkWmcxBcvPdMN1GZmpNuL6Tzgt4ANktYlZX8OnAIQEd8A7gUuAjYB+4DfTfbtkPR5YHVy3jURsSONIBfPy/Dyqwf4q+8+zXc//U4aJ09I423MzMac1BJERPw3oBGOCeBTBfbdCNyYQmiHaWnKtjl07Ox2gjAzS3gkNZBpnAxkE4SZmWU5QfDaFUSnu7eamQ1yggAaJ9czeUItHTv3lTsUM7OK4QQBSKKlqYFO32IyMxvkBJHINDa4DcLMLIcTRKKlabLbIMzMcjhBJDJNDezu7uGV/T3lDsXMrCI4QSTck8nM7FBOEIlMMkGfG6rNzLKcIBKZnNHUZmbmBDGoeepEJtbV+BaTmVnCCSIhKenq6sFyZmbgBHGIjAfLmZkNcoLI0dLkwXJmZgOcIHK0NE1m+96DdB/sK3coZmZl5wSRY7CrqxuqzcxSXXL0RknbJD1RYP+VktYljyck9Uk6Ptn3nKQNyb72tGIc6rWFg9xQbWaW5hXETcCiQjsj4rqImBsRc4GrgP8asqzou5P9bSnGeIiMR1ObmQ1KLUFExENAsetILwVuTiuWYp04bRJ1NXJDtZkZFdAGIWky2SuNO3OKA/i+pMckLRvh/GWS2iW1d3V1HVMstTViVqO7upqZQQUkCOD9wMNDbi+dHxHzgQuBT0n6tUInR8SKiGiLiLbm5uZjDsaD5czMsiohQSxhyO2liOhMnrcBdwPnjlYwLU0NboMwM6PMCULSdOBdwHdyyqZImjawDSwA8vaESkOmqYGX9hzgQK/HQphZdatLq2JJNwMXADMkdQCfBeoBIuIbyWEfAL4fEXtzTj0JuFvSQHz/HhHfSyvOoVqaJgOwddd+5syYMlpva2ZWcVJLEBGxtIhjbiLbHTa3bDNwdjpRjSx3sJwThJlVs0pog6goHixnZpblBDHEydMnUSOvLGdm5gQxRH1tDScfN8mD5cys6jlB5NHSNJkOd3U1syrnBJGHFw4yM3OCyKulqYEX9+ynt6+/3KGYmZWNE0QemcYG+vqDrbv3lzsUM7OycYLIY2CwnKfcMLNq5gSRx+C6EG6HMLMq5gSRx6zGSQDu6mpmVc0JIo+JdbWcOG0inbs8mtrMqpcTRAGZpgZfQZhZVXOCKKClabIbqc2sqjlBFJBpbGDLrm76+6PcoZiZlYUTRAEtTQ309AXbXjlQ7lDMzMrCCaKAwa6ubqg2syqVWoKQdKOkbZLyLhcq6QJJuyWtSx5X5+xbJGmjpE2SlqcV43BaB9eFcDuEmVWnNK8gbgIWjXDMjyNibvK4BkBSLXA9cCFwBrBU0hkpxpnXrEYnCDOrbqkliIh4CNhxFKeeC2yKiM0RcRC4BbikpMEVYfKEOk6YMsEJwsyqVrnbIN4h6XFJ90l6c1KWAV7IOaYjKctL0jJJ7ZLau7q6ShpcpqnBXV3NrGqVM0GsAWZHxNnAV4F7jqaSiFgREW0R0dbc3FzSADONDV6b2syqVtkSRETsiYhXk+17gXpJM4BOoDXn0JakbNS1JAsHRXgshJlVn7IlCEknS1KyfW4Sy3ZgNXCapFMlTQCWACvLEWOmsYEDvf28/OrBcry9mVlZ1aVVsaSbgQuAGZI6gM8C9QAR8Q3gQ8DvS+oFuoElkf1TvVfSFcAqoBa4MSKeTCvO4eSuC9E8bWI5QjAzK5vUEkRELB1h/9eArxXYdy9wbxpxHYncdSHmtjaWORozs9FV7l5MFS0zOFjODdVmVn2cIIZx3KR6jptU566uZlaVnCBG0NI02YPlzKwqOUGMIJN0dTUzqzZOECMYGCznsRBmVm2cIEbQ0tTA3oN97O7uKXcoZmajygliBC2e9tvMqpQTxAgGBss5QZhZtXGCGEGmcWBlOScIM6suThAjaJxcz5QJtR4sZ2ZVxwliBJLc1dXMqpITRBE8WM7MqpETRBEyjV5ZzsyqjxNEETJNDezu7uGV/R4LYWbVwwmiCANjIXwVYWbVJLUEIelGSdskPVFg/8ckrZe0QdIjks7O2fdcUr5OUntaMRZroKtrxw4nCDOrHmleQdwELBpm/y+Ad0XEmcDngRVD9r87IuZGRFtK8RUtd2U5M7NqkeaKcg9JmjPM/kdyXj4KtKQVy7GaMXUCE+tqnCDMrKpUShvEJ4D7cl4H8H1Jj0laNtyJkpZJapfU3tXVlUpwA2MhPFjOzKpJalcQxZL0brIJ4vyc4vMjolPSicD9kp6JiIfynR8RK0huT7W1taU2J3em0YPlzKy6lPUKQtJZwA3AJRGxfaA8IjqT523A3cC55YnwNS1NDR4sZ2ZVpWwJQtIpwF3Ab0XEz3LKp0iaNrANLADy9oQaTS1Nk9m+9yDdB/vKHYqZ2ahI7RaTpJuBC4AZkjqAzwL1ABHxDeBq4ATgnyQB9CY9lk4C7k7K6oB/j4jvpRVnsV6b1XUfv3LitDJHY2aWvjR7MS0dYf/vAb+Xp3wzcPbhZ5RX7sJBThBmVg0qpRdTxct4NLWZVRkniCKdOG0S9bVyQ7WZVQ0niCLV1oiZ093V1cyqhxPEEWjxYDkzqyJOEEfA60KYWTVxgjgCmaYGXtpzgAO9HgthZuOfE8QRGJjVdeuu/WWOxMwsfU4QR2BwXQg3VJtZFXCCOAKvrSznhmozG/9GTBCSaiX93WgEU+lOnj6JGuGurmZWFUZMEBHRx6FTcVet+toaZk73rK5mVh2KnYtpraSVwO3A3oHCiLgrlagqWKaxgQ53dTWzKlBsgpgEbAfek1MWZKfrriotTQ385Bc7yh2GmVnqikoQEfG7aQcyVmSaGnjx8f309vVTV+s2fjMbv4r6hZPUIuluSduSx52SWtIOrhJlGhvo6w+27vZYCDMb34r9E/ibwEpgVvL4j6Ss6gwMlvOUG2Y23hWbIJoj4psR0Zs8bgKaRzpJ0o3JFUfeJUOV9RVJmyStlzQ/Z9/lkp5NHpcXGWfqMk0eLGdm1aHYBLFd0m8mYyJqJf0m2UbrkdwELBpm/4XAacljGfB1AEnHk12i9G3AucBnJTUVGWuqZjVOAjwWwszGv2ITxMeBy4AXga3Ah4ARG64j4iFguC4/lwDfiqxHgUZJM4GFwP0RsSMidgL3M3yiGTUT62o5cdpEj6Y2s3FvxF5MkmqBSyPi4hTePwO8kPO6IykrVJ4vvmVkrz445ZRTUgjxcNl1IXwFYWbjW7EjqZeOQixHJSJWRERbRLQ1N4/YLFISmabJbqQ2s3Gv2FtMD0v6mqR3Spo/8CjB+3cCrTmvW5KyQuUVoaWpgS27uunvj3KHYmaWmmJHUs9Nnq/JKQsOHVl9NFYCV0i6hWyD9O6I2CppFfCFnIbpBcBVx/heJZNpbKCnL9j2ygFOnj6p3OGYmaWimDaIGuDrEXHbkVYu6WbgAmCGpA6yPZPqASLiG8C9wEXAJmAfScN3ROyQ9HlgdVLVNRFRMfNbvNbVdZ8ThJmNWyMmiIjol/QZ4IgTREQM23YREQF8qsC+G4Ebj/Q9R0Pr4LoQ3bSVORYzs7QU2wbxA0l/JqlV0vEDj1Qjq2CzvLKcmVWBYtsgPpI85/61H8DrShvO2DB5Qh0nTJngBGFm41qxs7memnYgY02mqcFdXc1sXBv2FlPS9jCw/eEh+76QVlBjQXawnEdTm9n4NVIbxJKc7aHdTCti6otyyTQ20Lmzm2w7u5nZ+DNSglCB7Xyvq0pL02QO9Pbz8qsHyx2KmVkqRkoQUWA73+uqkml8raurmdl4NFIj9dmS9pC9WmhItkleV/UIsdzBcnNbG8scjZlZ6Q2bICKidrQCGWsGEoTXhTCz8arYgXI2xHGT6jluUp3HQpjZuOUEcQxaPO23mY1jThDHINPU4FtMZjZuOUEcg4HBch4LYWbjkRPEMcg0NrD3YB+7u3vKHYqZWck5QRyDlqbJgGd1NbPxyQniGLQ0edpvMxu/Uk0QkhZJ2ihpk6Tlefb/g6R1yeNnknbl7OvL2bcyzTiPVqbxtcFyZmbjTbHrQRwxSbXA9cBvAB3AakkrI+KpgWMi4k9yjv9DYF5OFd0RMZcK1ji5nikTat3V1czGpTSvIM4FNkXE5og4CNwCXDLM8UuBm1OMp+QkuaurmY1baSaIDPBCzuuOpOwwkmYDpwIP5BRPktQu6VFJiwu9iaRlyXHtXV1dpYj7iLQ0TXYbhJmNS5XSSL0EuCMi+nLKZkdEG/BR4MuSXp/vxIhYERFtEdHW3Nw8GrEeItPoleXMbHxKM0F0Aq05r1uSsnyWMOT2UkR0Js+bgR9xaPtExWhpamB3dw+v7PdYCDMbX9JMEKuB0ySdKmkC2SRwWG8kSW8EmoD/ySlrkjQx2Z4BnAc8NfTcSjA4q6uvIsxsnEktQUREL3AFsAp4GrgtIp6UdI2ki3MOXQLcEofOV/EmoF3S48CDwLW5vZ8qyeBguR1OEGY2vqTWzRUgIu4F7h1SdvWQ15/Lc94jwJlpxlYqXlnOzMarSmmkHrNmTJ3AxLoaD5Yzs3HHCeIYDY6F8BWEmY0zThAlkGn0YDkzG3+cIErAg+XMbDxygiiBlqYGtu89SPfBvpEPNjMbI5wgSuDF3dmrhzOu/h7nXfsA96wtNB7QzGzscII4Rves7eTW9g4Agmx316vu2uAkYWZjnhPEMbpu1UYO9vYfUtbd08d1qzaWKSIzs9JwgjhGWwp0by1UbmY2VjhBHKNZyUjqYsvNzMYKJ4hjdOXC02morz2krKG+lisXnl6miMzMSiPVuZiqweJ52TWQrlu1kc5d3dRKfGHxWwbLzczGKl9BlMDieRkeXv4evrJ0Hn0RnDh9UrlDMjM7Zk4QJbTgjJM4blIdt7W/MPLBZmYVzgmihCbV17J4Xob7nniR3fu8wpyZjW2pJghJiyRtlLRJ0vI8+39HUpekdcnj93L2XS7p2eRxeZpxltJlba0c7O1n5fot5Q7FzOyYpJYgJNUC1wMXAmcASyWdkefQWyNibvK4ITn3eOCzwNuAc4HPSmpKK9ZSevOs43jTzOO43beZzGyMS/MK4lxgU0RsjoiDwC3AJUWeuxC4PyJ2RMRO4H5gUUpxlpQkLmtrYX3Hbp7euqfc4ZiZHbU0E0QGyP0zuiMpG+qDktZLukNS6xGei6RlktoltXd1dZUi7mO2eG6GCbU1bqw2szGt3I3U/wHMiYizyF4l/OuRVhARKyKiLSLampubSx7g0WiaMoHfOOMk7lnbyYFeTwFuZmNTmgmiE2jNed2SlA2KiO0RcSB5eQNwTrHnVrrL3trKzn09/PDpbeUOxczsqKSZIFYDp0k6VdIEYAmwMvcASTNzXl4MPJ1srwIWSGpKGqcXJGVjxvm/MoOZ0yf5NpOZjVmpJYiI6AWuIPvD/jRwW0Q8KekaSRcnh31a0pOSHgc+DfxOcu4O4PNkk8xq4JqkbMyorREfOqeFh37WxdbdntnVzMYeRUS5YyiZtra2aG9vL3cYg57fvpd3Xfcjrlx4Op9696+UOxwzs8NIeiwi2vLtK3cj9bg2+4QpvP11x3Nb+wuMp0RsZtXBCSJll7W18vz2ffz0F2PqDpmZmRNE2i58y0ymTqzjtmTdajOzscIJImUNE2p5/9mzuHfDVl7Z7wn8zGzscIIYBZe1tdDd08d/rt9a7lDMzIrmBDEK5rY28oaTpnpMhJmNKU4QoyA7gV8ra3+5i2dfeqXc4ZiZFcUJYpQsnpehrkbc/pgbq81sbHCCGCUzpk7kvW86kbvWdNDT11/ucMzMRuQEMYoua2vl5VcP8uAznsDPzCqfE8QoetcbmmmeNtFjIsxsTHCCGEV1tTV8cH4LD27cxrZX9pc7HDOzYTlBjLIPt7XQ1x/cvWZMLW9hZlXICWKUvb55Km2zm7jVE/iZWYVzgiiDy9pa2dy1lzW/3FnuUMzMCnKCKIP/ddZMJk+o5bbVbqw2s8qVaoKQtEjSRkmbJC3Ps/9PJT0lab2kH0qanbOvT9K65LFy6Llj2ZSJdbzvrJn85/ot7D3QW+5wzMzySi1BSKoFrgcuBM4Alko6Y8hha4G2iDgLuAP425x93RExN3lczDhzWVsrew/2ce8GT+BnZpUpzSuIc4FNEbE5Ig4CtwCX5B4QEQ9GxL7k5aNAS4rxVJRzZjfxuhlTuN1jIsysQqWZIDJA7vSlHUlZIZ8A7st5PUlSu6RHJS0udJKkZclx7V1dXccW8SiSxIfbWvnpczvY3PVqucMxMztMRTRSS/pNoA24Lqd4drKQ9keBL0t6fb5zI2JFRLRFRFtzc/MoRFs6H5yfobZG3OEJ/MysAqWZIDqB1pzXLUnZIST9OvAXwMURcWCgPCI6k+fNwI+AeSnGWhYnHjeJC97QzJ1rOuj1BH5mVmHSTBCrgdMknSppArAEOKQ3kqR5wD+TTQ7bcsqbJE1MtmcA5wFPpRhr2Xy4rZWX9hzgoWfHzu0xM6sOqSWIiOgFrgBWAU8Dt0XEk5KukTTQK+k6YCpw+5DurG8C2iU9DjwIXBsR4zJBvOeNJ3LClAkeE2FmFacuzcoj4l7g3iFlV+ds/3qB8x4Bzkwztkoxoa6GS+dn+ObDz7H91QOcMHViuUMyMwMqpJG62n24rZXe/uA9f/8jTl3+Xc679gHuWevJ/MysvFK9grDiPLVlDxLs7s6Oqu7c1c1Vd20AskuVmpmVg68gKsB1qzYydGLX7p4+rlu1sTwBmZnhBFERtuzqPqJyM7PR4ARRAWY1NuQtb5xc7zUjzKxsnCAqwJULT6ehvvaQMgl27uvht2/8KS/s2FfgTDOz9DhBVIDF8zJ88dIzyTQ2ICDT2MCXPnQ2n7/kzax5ficL/uEhbvjxZvr6fTVhZqNH4+kWRltbW7S3t5c7jJLasqubv7znCX74zDbObpnOtR88izfNPK7cYZnZOCHpsWTeu8P4CqLCzWps4IbL2/jq0nl07Ozm/V/9b65b9Qz7e/rKHZqZjXNOEGOAJN5/9ix+8KfvYvG8DNc/+HMu+scf8+jm7eUOzczGMSeIMaRpygT+7sNn8/8+8TZ6+vtZsuJRrrprA3v295Q7NDMbh9wGMUbtO9jLl3/wLDf8eDMzpk7k84vfQvfB7OC6Lbu6mdXYwJULT/dIbDMb1nBtEE4QY9yGjt185s71PL11DzWC3I5ODfW1fPHSM50kzKwgN1KPY2e2TGflFedx3KQ6hvaC7e7p42+/90x5AjOzMc+T9Y0D9bU1vLK/N+++Lbv386GvP8I5s5s4Z3YT82c3McNTiptZEZwgxolZjQ105pm7aerEWvoj+ObDz/HPD20GYM4Jk5mfJIy22cdz2olTqakR96ztPOY2jFLUYWaVIdUEIWkR8I9ALXBDRFw7ZP9E4FvAOcB24CMR8Vyy7yrgE0Af8OmIWJVmrGPdlQtP56q7NtCdMz6iob6Wv1qcbYPY39PHk1t20/7cTh57ficP/ayLu9Zk15yYNqmOWY0N/Hzbq/Qm96mOZsrxe9Z2HhLD0U5bXimJynWMzzoqIYZKqmM4qTVSS6oFfgb8BtBBdo3qpblLh0r6A+CsiPikpCXAByLiI5LOAG4GzgVmAT8A3hARw44Oq8ZG6lxH8mWJCH65Yx+PPb+T9ud3ctvqFwaTQ64awewTpjB5Qm3yqDv0eWItk+vrmDKxlq89sIld3Yd3uW2eOpGbPv5WJtTWUFdbQ32tqK+toa5G1NfVUF+TLautEd9ZtyVvojuSxvahicp1uI5KiqGS6oAy9WKS9A7gcxGxMHl9FUBEfDHnmFXJMf8jqQ54EWgGlucem3vccO9Z7QniWJy6/LsU+ia8/+xZ7DvQy76Dfew7OPCc3d57sI+Dvf2px1cjOH5Ktu1EypbpkG0N7tu2Zz99ef5jamvErMZJg8fn1pM9V4Pbv9yxL2/CrKsRs0+YXFTMz28/sjo0EFCO517eW7COOTOmFBVH2nWcWmQdv6iAOiohhrTryDQ28PDy9xRVBwyfINK8xZQBXsh53QG8rdAxEdEraTdwQlL+6JBz86ZEScuAZQCnnHJKSQKvRoXaMDKNDXx16bxhz+3t62dfTx8LvvQQL+7Zf9j+E6ZM4K8/cCY9ff309vfT0xfZ7eS5py/o7eunp6+frzywKe979AcsePNJOQsrxeD24HNSdvtjHXnr6OsP3jr7eAIGp1EfqC4idzvY/PLe/P+t/cEbi5wL6+ddR1BHgey8adurBes4/aRpRcWRdh2nnTS1qDqerYA6KiGGtOso5ToyY76ROiJWACsgewVR5nDGrEJtGFcuPH3Ec+tqaziutoblF74xbx1/+b4zWPSWk4uK4841nQUT1Rc+cGZRdTzy8+0F6/jSR+YWVcfaXz5QsI7rPzq/qDrWlaKOa4ep42OVUcc/feycouo4rwLqqIQY0q6j0PoyRyPNcRCdQGvO65akLO8xyS2m6WQbq4s510oo35TjR3ovsxR15Fsbo9hE5Tpcx1iIoZLqGEmabRB1ZBup30v2x3018NGIeDLnmE8BZ+Y0Ul8aEZdJejPw77zWSP1D4DQ3UleHSund4TrGZx2VEEMl1VG2qTYkXQR8mWw31xsj4q8lXQO0R8RKSZOAfwPmATuAJRGxOTn3L4CPA73AH0fEfSO9nxOEmdmR8VxMZmaWl+diMjOzI+YEYWZmeTlBmJlZXk4QZmaW17hqpJbUBTx/lKfPAF4uYThpcZylN1ZidZylNVbihHRjnR0Rzfl2jKsEcSwktRdqya8kjrP0xkqsjrO0xkqcUL5YfYvJzMzycoIwM7O8nCBes6LcARTJcZbeWInVcZbWWIkTyhSr2yDMzCwvX0GYmVleThBmZpZX1SUISYskbZS0SdLyPPsnSro12f8TSXPKEGOrpAclPSXpSUl/lOeYCyTtlrQueVw92nEmcTwnaUMSw2EzJSrrK8nnuV5ScSvUlDbG03M+p3WS9kj64yHHlO3zlHSjpG2SnsgpO17S/ZKeTZ6bCpx7eXLMs5IuL0Oc10l6Jvm3vVtSY4Fzh/2ejEKcn5PUmfPve1GBc4f9fRilWG/NifM5SesKnJv+ZxoRVfMgO+34z4HXAROAx4EzhhzzB8A3ku0lwK1liHMmMD/ZnkZ2XY2hcV4A/GcFfKbPATOG2X8RcB/ZpZ/fDvykAr4DL5IdHFQRnyfwa8B84Imcsr8Flifby4G/yXPe8cDm5Lkp2W4a5TgXAHXJ9t/ki7OY78koxPk54M+K+G4M+/swGrEO2f/3wNXl+kyr7QriXGBTRGyOiIPALcAlQ465BPjXZPsO4L3Kt5p8iiJia0SsSbZfAZ6mwJrcY8AlwLci61GgUdLMMsbzXuDnEXG0I+5LLiIeIrseSq7c7+G/AovznLoQuD8idkTETuB+YNFoxhkR34+I3uTlo2RXfyyrAp9nMYr5fSip4WJNfncuA25OM4bhVFuCyAAv5Lzu4PAf3sFjki/+buCEUYkuj+QW1zzgJ3l2v0PS45LuS1bhK4cAvi/pMUnL8uwv5jMfTUso/D9cJXyeA06KiK3J9ovASXmOqbTP9uNkrxbzGel7MhquSG6F3Vjgll2lfZ7vBF6KiGcL7E/9M622BDGmSJoK3El2Rb09Q3avIXub5Gzgq8A9ox1f4vyImA9cCHxK0q+VKY4RSZoAXAzcnmd3pXyeh4ns/YSK7o+erADZC3y7wCHl/p58HXg9MBfYSvbWTaVbyvBXD6l/ptWWIDqB1pzXLUlZ3mOUXVd7OrB9VKLLIamebHL4dkTcNXR/ROyJiFeT7XuBekkzRjlMIqIzed4G3E32Mj1XMZ/5aLkQWBMRLw3dUSmfZ46XBm7FJc/b8hxTEZ+tpN8B3gd8LElmhynie5KqiHgpIvoioh/4vwXevyI+Txj87bkUuLXQMaPxmVZbglgNnCbp1OSvySXAyiHHrAQGeoN8CHig0Jc+Lcm9x38Bno6ILxU45uSBthFJ55L9txzVRCZpiqRpA9tkGyyfGHLYSuC3k95Mbwd259w6GW0F/yKrhM9ziNzv4eXAd/IcswpYIKkpuWWyICkbNZIWAZ8BLo6IfQWOKeZ7kqoh7V4fKPD+xfw+jJZfB56JiI58O0ftM02zBbwSH2R71fyMbG+Fv0jKriH7BQeYRPYWxCbgp0qdAkEAAAL6SURBVMDryhDj+WRvKawH1iWPi4BPAp9MjrkCeJJsT4tHgV8tQ5yvS97/8SSWgc8zN04B1yef9wagrUz/7lPI/uBPzymriM+TbNLaCvSQve/9CbLtXj8EngV+AByfHNsG3JBz7seT7+om4HfLEOcmsvftB76nAz0AZwH3Dvc9GeU4/y35/q0n+6M/c2icyevDfh9GO9ak/KaB72bOsaP+mXqqDTMzy6vabjGZmVmRnCDMzCwvJwgzM8vLCcLMzPJygjAzs7ycIMwSkl5NnudI+miJ6/7zIa8fKWX9ZmlwgjA73BzgiBJEMvJ1OIckiIj41SOMyWzUOUGYHe5a4J3JPPt/Iqk2WfdgdTLZ2/+GwTUkfixpJfBUUnZPMnnakwMTqEm6FmhI6vt2UjZwtaKk7ieSuf0/klP3jyTdoex6C9/OGel9rbJrhayX9Hej/ulY1Rjprx6zarSc7NoB7wNIfuh3R8RbJU0EHpb0/eTY+cBbIuIXyeuPR8QOSQ3Aakl3RsRySVdExNw873Up2QnkzgZmJOc8lOybB7wZ2AI8DJwn6WmyU0W8MSJCBRboMSsFX0GYjWwB2fmk1pGddv0E4LRk309zkgPApyUNTNfRmnNcIecDN0d2IrmXgP8C3ppTd0dkJ5hbR/bW125gP/Avki4F8s5/ZFYKThBmIxPwhxExN3mcGhEDVxB7Bw+SLiA7ydo7Ijtt+Fqyc3sdrQM5231kV27rJTtr5x1kZ1D93jHUbzYsJwizw71CdqnXAauA30+mYEfSG5IZNIeaDuyMiH2S3kh2idUBPQPnD/Fj4CNJO0cz2SUof1oosGSNkOmRnZL8T8jemjJLhdsgzA63HuhLbhXdBPwj2ds7a5KG4i7yLwH6PeCTSTvBRrK3mQasANZLWhMRH8spvxt4B9lZOQP4TES8mCSYfKYB35E0ieyVzZ8e3X+i2cg8m6uZmeXlW0xmZpaXE4SZmeXlBGFmZnk5QZiZWV5OEGZmlpcThJmZ5eUEYWZmef1/GGFTEMsGapwAAAAASUVORK5CYII=\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def jacobi_plot():\n", + " \n", + " iterations = np.array([])\n", + " errors = np.array([])\n", + " \n", + " A = np.matrix(\"4 -1 -1; -2 6 1; -1 1 7\")\n", + " b = np.matrix(\"3; 9; -6\")\n", + " trueresult = np.linalg.solve(A, b)\n", + "\n", + " for n in range(32):\n", + " \n", + " result = jacobi(A.copy(), b, n)\n", + " \n", + " error = np.abs(np.sum(result-trueresult))\n", + " \n", + " iterations = np.append(iterations, n)\n", + " errors = np.append(errors, error)\n", + " \n", + " if np.isclose(error, 0):\n", + " print(\"Using:\")\n", + " print(f\"A=\\n{A}\")\n", + " print(f\"b=\\n{b}\")\n", + " print(f\"Converged in {n}:\")\n", + " print(f\"x=\\n{result}\")\n", + " print(f\"|error|=\\n{error}\")\n", + " break\n", + "\n", + " plt.figure()\n", + " plt.xlabel('Iterations')\n", + " plt.ylabel('Error')\n", + " plt.plot(iterations, errors, 'o-')\n", + "\n", + "jacobi_plot()" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "We find that Jacobi iteration yields an approximate solution in 18 iterations where the solution is very close to zero." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Gauss-Seidel iteration, Ax=b ###\n", + "\n", + "Below are tests for Gauss-Seidel iteration. Input A and b are given such that there exists an exact solution y for x under the condition of Ax=b. This means that the test explores convergence ||Ax-b|| and the approximation error ||x-y|| at the same time." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 6, + "outputs": [ + { + "name": "stdout", + "text": [ + "Using:\n", + "A=\n", + "[[ 4 -1 -1]\n", + " [-2 6 1]\n", + " [-1 1 7]]\n", + "b=\n", + "[[ 3]\n", + " [ 9]\n", + " [-6]]\n", + "Converged in 9:\n", + "x=\n", + "[[ 1.]\n", + " [ 2.]\n", + " [-1.]]\n", + "|error|=\n", + "1.8414112457065812e-09\n" + ], + "output_type": "stream" + }, + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5RddX338fdnbsmQhOSMmSBJJmeiDUEESTIH1GKtvQixj5VI9RGqSKsuapXaKz5g11KLq8IqXa22pWqqFNtH8WkRaahoxKJFUSSThIuBRmLIZSZoBiYJuWcu3+ePsyecTM5kJsk5s8/l81rrrDn7t/c++5uzYD6z9++3f1sRgZmZ2WgNaRdgZmaVyQFhZmZFOSDMzKwoB4SZmRXlgDAzs6Ka0i6glGbPnh2dnZ1pl2FmVjXWrl37XES0F1tXUwHR2dlJd3d32mWYmVUNSVvHWudLTGZmVpQDwszMinJAmJlZUQ4IMzMrygFhZmZFlS0gJHVI+o6kJyVtkPSHRbaRpL+TtEnS45KWFay7RtLTyeuactV5z/peLrnlARbe8HUuueUB7lnfW65DmZlVlXIOcx0E/jQi1kmaAayVdH9EPFmwzZuARcnr1cBngFdLagM+BuSASPZdFRG7SlngPet7ufHuJzg4MARA7+6D3Hj3EwCsWDqvlIcyM6s6ZTuDiIhnI2Jd8n4v8BQw+rfu5cC/RN7DwCxJZwOXAfdHRH8SCvcDy0td462rNx4NhxEHB4a4dfXGUh/KzKzqTEofhKROYCnwo1Gr5gHbC5Z7krax2ot99rWSuiV19/X1nVRdO3YfPKl2M7N6UvaAkDQd+CrwRxHxQqk/PyJWRkQuInLt7UXvFh/T3FmtJ9VuZlZPyhoQkprJh8OXIuLuIpv0Ah0Fy/OTtrHaS+r6yxbT2tx4TFtrcyPXX7a41IcyM6s65RzFJOALwFMR8TdjbLYKeHcymuk1wJ6IeBZYDVwqKSMpA1yatJXUiqXzuPmKC5g9vQWAtmkt3HzFBe6gNjOjvKOYLgGuBp6Q9GjS9hFgAUBEfBa4D/gNYBNwAPjdZF2/pE8Aa5L9boqI/nIUuWLpPJaf/1Iu+Phq3p6b73AwM0uULSAi4vuAxtkmgA+Ose524PYylHacqc2NvHLuTNZtLekoWjOzquY7qRO5bIbHevZweHBo/I3NzOqAAyLRlc1wZHCYDTtKPtDKzKwqOSASXdkMgC8zmZklHBCJOWdOpaOtle4tDggzM3BAHKNrQYa123aR7zs3M6tvDogCXZ1t9O09TM8uT7VhZuaAKNC1IN8P0b21LLdcmJlVFQdEgcUvncH0KU2sdUe1mZkDolBjg1i6YJY7qs3McEAcZ9mCDBt/vpe9hwbSLsXMLFUOiFFynRki4NHtu9MuxcwsVQ6IUZZ0zELCl5nMrO45IEaZMbWZxWfNYN02B4SZ1TcHRBG5zgzrt+1maNg3zJlZ/XJAFNGVzbDv8CAbf7Y37VLMzFLjgCgil20DYK0vM5lZHSvnI0dvl7RT0o/HWH+9pEeT148lDUlqS9ZtkfREsq67XDWOZX6mlfYZU1i7xXdUm1n9KucZxB3A8rFWRsStEbEkIpYANwL/Peqxor+SrM+VscaiJB2duM/MrF6VLSAi4kFgon+CXwXcWa5aTkWuM8P2/oPsfOFQ2qWYmaUi9T4ISWeQP9P4akFzAN+StFbStePsf62kbkndfX19JatrWfIAIc/LZGb1KvWAAH4TeGjU5aXXRcQy4E3AByW9fqydI2JlROQiItfe3l6yol4590xamhocEGZWtyohIK5k1OWliOhNfu4EvgZcPNlFTWlq5ML5M+l2QJhZnUo1ICTNBH4Z+I+CtmmSZoy8By4Fio6EKrdl2Qwbduzh0MBQGoc3M0tVOYe53gn8EFgsqUfSeyW9X9L7CzZ7K/CtiNhf0HYW8H1JjwGPAF+PiG+Wq84TyWXbGBgKnujdk8bhzcxS1VSuD46IqyawzR3kh8MWtm0GLixPVSdn2YJZQH7ivos621KuxsxsclVCH0TFesn0KSycPc0d1WZWlxwQ4+jKZli3bRcRnrjPzOqLA2IcXdkM/fuP8Mxz+8ff2MyshjggxpHzDXNmVqccEON4eft0zpza5IAws7rjgBhHQ4NYls04IMys7jggJiCXzfD0zn3sPnAk7VLMzCaNA2ICRibuW79td8qVmJlNHgfEBCzpmEVjg3yZyczqigNiAs5oaeK8s8+ke6ufMGdm9cMBMUFd2QyPbd/DwNBw2qWYmU0KB8QEdWUzHBwY4qlnX0i7FDOzSeGAmKAu3zBnZnXGATFBc2e1MnfmVAeEmdUNB8RJ8A1zZlZPHBAnoSub4dk9h9ix+2DapZiZlV05nyh3u6Sdkoo+LlTSGyTtkfRo8vpowbrlkjZK2iTphnLVeLJy2fxDg3wWYWb1oJxnEHcAy8fZ5nsRsSR53QQgqRG4DXgTcB5wlaTzyljnhJ179gxamxsdEGZWF8oWEBHxIHAqd5ZdDGyKiM0RcQT4CnB5SYs7Rc2NDSzpmOWAMLO6kHYfxGslPSbpG5JembTNA7YXbNOTtBUl6VpJ3ZK6+/r6ylkrkO+HePLZF9h/eLDsxzIzS1OaAbEOyEbEhcDfA/ecyodExMqIyEVErr29vaQFFtOVzTA0HDzW44n7zKy2pRYQEfFCROxL3t8HNEuaDfQCHQWbzk/aKsKyBfkb5tb5MpOZ1bjUAkLSSyUpeX9xUsvzwBpgkaSFklqAK4FVadU52swzmlk0ZzrdDggzq3FN5fpgSXcCbwBmS+oBPgY0A0TEZ4G3Ab8vaRA4CFwZEQEMSroOWA00ArdHxIZy1Xkqcp0Zvv74swwPBw0NSrscM7OyKFtARMRV46z/B+Afxlh3H3BfOeoqhWULMtz5yHY29e3jnLNmpF2OmVlZpD2KqSp54j4zqwcOiFOwcPY02qa1OCDMrKY5IE6BJJYt8MR9ZlbbHBCnKNeZ4Znn9vP8vsNpl2JmVhYOiFPkfggzq3UOiFN0wbyZNDeKtdscEGZWmxwQp2hqcyPnz5vpO6rNrGY5IE5D14IMj/Xs4fDgUNqlmJmVnAPiNOQ6MxwZHGbDjhfSLsXMrOQcEKfBE/eZWS1zQJyGOWdOpaOtle4tDggzqz0OiNOUy7axdtsu8vMMmpnVDgfEaVqWzdC39zDb+w+mXYqZWUk5IE5TV9IPsXbbqTx+28yscjkgTtPil85g+pQm31FtZjWnbAEh6XZJOyX9eIz175T0uKQnJP1A0oUF67Yk7Y9K6i5XjaXQ2CCWLpjljmozqznlPIO4A1h+gvXPAL8cERcAnwBWjlr/KxGxJCJyZaqvZLqyGTb+fC97Dw2kXYqZWcmULSAi4kFgzAvzEfGDiBj5s/thYH65aim3rmyGCFi/bXfapZiZlUyl9EG8F/hGwXIA35K0VtK1J9pR0rWSuiV19/X1lbXIsSzpmEWDPLOrmdWWsj2TeqIk/Qr5gHhdQfPrIqJX0hzgfkn/k5yRHCciVpJcnsrlcqncjDBjajOLX3om6zyzq5nVkFTPICS9Cvg8cHlEPD/SHhG9yc+dwNeAi9OpcOK6srNYv203Q8O+Yc7MakNqASFpAXA3cHVE/KSgfZqkGSPvgUuBoiOhKkku28a+w4Ns/NnetEsxMyuJsl1iknQn8AZgtqQe4GNAM0BEfBb4KPAS4B8lAQwmI5bOAr6WtDUBX46Ib5arzlJ58Qlz/Zw398yUqzEzO31lC4iIuGqc9e8D3lekfTNw4fF7VLb5mVbaZ0xh7dZdXP3azrTLMTM7bZUyiqnqSSKXzfgRpGZWMxwQJdSVzbC9/yA7XziUdilmZqfNAVFCL/ZD+CzCzKqfA6KEXjl3Ji1NDQ4IM6sJDogSamlq4ML5M+l2QJhZDXBAlFhXto0NO/ZwaGAo7VLMzE6LA6LEurIZBoaCx3v2pF2KmdlpcUCUmDuqzaxWOCBKrG1aCy+bPc0BYWZVb9yAkNQo6a8no5hasSybYd22XUR44j4zq17jBkREDHHsVNw2jlw2Q//+Izzz3P60SzEzO2UTnYtpvaRVwL8DR3/rRcTdZamqyo30Q3Rv3cXL2qenXI2Z2amZaB/EVOB54FeB30xeby5XUdXu5e3TmdnazDr3Q5hZFZvQGURE/G65C6klDQ1i2YJZ7qg2s6o2oTMISfMlfU3SzuT1VUnzy11cNevKZnh65z52HziSdilmZqdkopeY/hlYBcxNXvcmbTaGrmwbAOu37U65EjOzUzPRgGiPiH+OiMHkdQfQPt5Okm5PzjiKPjJUeX8naZOkxyUtK1h3jaSnk9c1E6yzYlzYMZPGBtG9tT/tUszMTslEA+J5Se9K7ololPQu8p3W47kDWH6C9W8CFiWva4HPAEhqI/+I0lcDFwMfk5SZYK0V4YyWJl4590z3Q5hZ1ZpoQLwH+N/Az4BngbcB43ZcR8SDwIn+hL4c+JfIexiYJels4DLg/ojoj4hdwP2cOGgq0rIFGR7bvoeBoeG0SzEzO2kTupMauCIi3hIR7RExJyJWRMS2Ehx/HrC9YLknaRurvVh910rqltTd19dXgpJKpyub4eDAEE89+0LapZiZnbSJ3kl91STUckoiYmVE5CIi194+brfIpMp1euI+M6teE73E9JCkf5D0S5KWjbxKcPxeoKNgeX7SNlZ7VTl7ZitzZ051QJhZVZroVBtLkp83FbQF+TurT8cq4DpJXyHfIb0nIp6VtBr4ZEHH9KXAjad5rFR0dbbRvcUjmcys+owbEJIagM9ExL+d7IdLuhN4AzBbUg/5kUnNABHxWeA+4DeATcABko7viOiX9AlgTfJRN0VEVf6W7Vowi3sf28GO3QeZO6s17XLMzCZs3ICIiGFJHwZOOiAi4oR9F5GfD/uDY6y7Hbj9ZI9ZaUZumOveuou3OCDMrIpMtA/i25L+TFKHpLaRV1krqxGvOHsGrc2NnrjPzKrORPsg3pH8LPxrP4CXlbac2tPU2MCSDk/cZ2bVZ6KzuS4sdyG1LNeZ4R+/+1P2Hx5k2pSJZrKZWbpOeIkp6XsYef/2Ues+Wa6ias2ybIah4eCxHk/cZ2bVY7w+iCsL3o8eZlp1U1+kZVlHcsPcFl9mMrPqMV5AaIz3xZZtDDPPaOacs6azdpsDwsyqx3gBEWO8L7ZsJ9CVzbBu6y6Gh/21mVl1GC8gLpT0gqS9wKuS9yPLF0xCfTVj2YIMLxwaZFPfvrRLMTObkBMOqYmIxskqpNblOvO3jazduotzzpqRcjVmZuOb6I1ydpo6X3IGL5nWQrc7qs2sSjggJokklmUzrHNHtZlVCQfEJOrKZnjmuf08v+9w2qWYmY3LATGJurJ+gJCZVQ8HxCS6YN5Mmhvl+yHMrCo4ICbR1OZGzp830zO7mllVcEBMslw2w2M9ezg8OJR2KWZmJ1TWgJC0XNJGSZsk3VBk/d9KejR5/UTS7oJ1QwXrVpWzzsnUlc1wZHCYDTteSLsUM7MTKtvc05IagduANwI9wBpJqyLiyZFtIuKPC7b/A2BpwUccjIgl1Jhl2Rcn7lu2IDPO1mZm6SnnGcTFwKaI2BwRR4CvAJefYPurgDvLWE9FmDNjKgvazvBIJjOreOUMiHnA9oLlnqTtOJKywELggYLmqZK6JT0sacVYB5F0bbJdd19fXynqLruubIa123aRfyS3mVllqpRO6iuBuyKisOc2GxE54LeBT0l6ebEdI2JlROQiItfe3j4ZtZ62rmyGvr2H2d5/MO1SzMzGVM6A6AU6CpbnJ23FXMmoy0sR0Zv83Ax8l2P7J6ra0RvmtvWnXImZ2djKGRBrgEWSFkpqIR8Cx41GknQukAF+WNCWkTQleT8buAR4cvS+1eqcs2YwY0qTJ+4zs4pWtlFMETEo6TpgNdAI3B4RGyTdBHRHxEhYXAl8JY69IP8K4HOShsmH2C2Fo5+qXWODWLJgljuqzayilS0gACLiPuC+UW0fHbX88SL7/YAafyBRVzbDp//rafYeGmDG1Oa0yzEzO06ldFLXnVy2jQhYv233+BubmaXAAZGSCztm0iDP7GpmlcsBkZIZU5tZ/NIzHRBmVrEcECnKZTOs37aLoWHfMGdmlccBkaKubIb9R4bY+LO9aZdiZnYcB0SKXnzCnG+YM7PK44BI0fxMK3NmTHE/hJlVJAdEiiQdnbjPzKzSOCBS1pXNsL3/IDtfOJR2KWZmx3BApOzFfgifRZhZZXFApOyVc2cypamBbgeEmVUYB0TKWpoauHC+J+4zs8rjgKgAy7IZNuzYw6GBofE3NjObJA6ICpDLZhgYCh7v2ZN2KWZmRzkgKsAyd1SbWQUqa0BIWi5po6RNkm4osv53JPVJejR5va9g3TWSnk5e15SzzrS1TWuhfUYLn/72T1h4w9e55JYHuGf9WE9nNTObHGV7YJCkRuA24I1AD7BG0qoiT4b7fxFx3ah924CPATkggLXJvjX5J/Y963vp3zfAUPJQvd7dB7nx7icAWLF0XpqlmVkdK+cZxMXApojYHBFHgK8Al09w38uA+yOiPwmF+4HlZaozdbeu3ng0HEYcHBji1tUbU6rIzKy8ATEP2F6w3JO0jfZbkh6XdJekjpPcF0nXSuqW1N3X11eKuifdjt0HT6rdzGwypN1JfS/QGRGvIn+W8MWT/YCIWBkRuYjItbe3l7zAyTB3VutJtZuZTYZyBkQv0FGwPD9pOyoino+Iw8ni54Guie5bS66/bDGtzY3HtLU0NnD9ZYtTqsjMrLwBsQZYJGmhpBbgSmBV4QaSzi5YfAvwVPJ+NXCppIykDHBp0laTViydx81XXMC8Wa0IaGoQUnDB/Jlpl2Zmdaxso5giYlDSdeR/sTcCt0fEBkk3Ad0RsQr4kKS3AINAP/A7yb79kj5BPmQAboqImn6qzoql846OWOrZdYAVtz3E+77Yzdc+8IvMOqMl5erMrB4ponaeh5zL5aK7uzvtMkpi7dZ+rlr5I3KdGb74notpbky7u8jMapGktRGRK7bOv3UqVFe2jZuvuIAf/PR5/uLeDWmXY2Z1qGyXmOz0/VbXfH6ycy+f++/NnHPWDN792s60SzKzOuIziAr34cvO5ddfcRZ/ce+TfO/p6rzPw8yqkwOiwjU2iE9duYRFc6bzgS+t46d9+9IuyczqhAOiCkyf0sTnr8nR0tjA+77Yze4DR9IuyczqgAOiSszPnMHnru6id9dBPvjldQwMDaddkpnVOAdEFcl1tvHJKy7goU0e2WRm5edRTFXmbV3zedojm8xsEvgMogrlRzbN8cgmMysrB0QVyo9sWsqiOdP5oEc2mVmZOCCq1PQpTfzTu3M0e2STmZWJA6KKdbR5ZJOZlY8DosoVjmy66d7Rj/s2Mzt1HsVUA97WNZ+nf76Xzz24mUVnTffIJjMrCZ9B1IgPL/fIJjMrLQdEjRg9smmzRzaZ2Wkqa0BIWi5po6RNkm4osv5PJD0p6XFJ/yUpW7BuSNKjyWvV6H3teIUjm977xW72HBhIuyQzq2JlCwhJjcBtwJuA84CrJJ03arP1QC4iXgXcBfxVwbqDEbEkeb2lXHXWmo62M/js1V307DrAB7681iObzOyUlfMM4mJgU0RsjogjwFeAyws3iIjvRMSBZPFhYH4Z66kbF3W28cm3emSTmZ2ecgbEPGB7wXJP0jaW9wLfKFieKqlb0sOSVoy1k6Rrk+26+/rcOTvi7bkOfu/1L+NfH97Kv/5wS9rlmFkVqohhrpLeBeSAXy5ozkZEr6SXAQ9IeiIifjp634hYCawEyOVyMSkFV4kPLz+XTTv38fF7n2Th7Om8btHstEsysypSzjOIXqCjYHl+0nYMSb8O/Dnwlog4PNIeEb3Jz83Ad4GlZay1JjU2iE9ftZRfaJ/OB7601iObzOyklDMg1gCLJC2U1AJcCRwzGknSUuBz5MNhZ0F7RtKU5P1s4BLAF9NPwcjT6EbmbPLIJjObqLIFREQMAtcBq4GngH+LiA2SbpI0MirpVmA68O+jhrO+AuiW9BjwHeCWiHBAnKKRkU3bPbLJzE6CImrnsn0ul4vu7u60y6hY/969nevvepyrX5PlEyvOT7scM6sAktZGRK7YuoropLbJ8fZcB5t27uNzD27mnLOmc7XnbDKzE/BUG3Xmw8vP5dfOncPH732S7z/9XNrlmFkFc0DUGY9sMrOJckDUoZGRTU0e2WRmJ+CAqFMjT6PbvuuAn0ZnZkU5IOrYyJxN39/0HJ/4T48iNrNjeRRTnXt7roOnd+5j5YObWTTHI5vM7EUOCOP/LD+Xn+7cx0f/YwOf+vbT9O8/wtxZrVx/2WJWLD3R/IpmVst8iclobBBvPG8OAM/vP0IAvbsPcuPdT3DP+uOmzzKzOuGAMAD+/oGfMvqe+oMDQ/zV6v9JpR4zS58vMRkAO3YfHKP9EFd/4Udc1NnGRZ1tLOmYRWtL4yRXZ2ZpcEAYAHNntdJbJCSmtTTSt/cwf/vtnxABTQ3i/Hkzuagzw0WdbeQ622ib1pJCxWZWbg4IA+D6yxZz491PcHBg6Ghba3Mjf/nWC1ixdB57Dgywbtsu1mzpZ82Wfr74g6380/eeAeDl7dO4eGEbuWz+LKOjrRVJaf1TzKxEPJurHXXP+l5uXb2RHbsPjjuK6dDAED/u3cMjW/rp3rKL7i39vHBoEICzzpxCrrONi7IZcp1tvOLsM2lscGCYVaITzebqgLCSGB4Ont65LwmMftY808+OPYeA/NQey7IZLspmuGhhvh9jarP7McwqgQPCUtG7+2A+LLb0s+aZXWz8+V4Amhvz/RgXJ30YuWyGTEE/xsmcyZjZ6UktICQtBz4NNAKfj4hbRq2fAvwL0AU8D7wjIrYk624E3gsMAR+KiNXjHc8BUdn2HBhg7bZ+Hnkmf0nq8Z49HEnmgPqFOdO5qLONBgV3re3l8OCLc0O1Njdy8xUXTHpIVEpQuY7Kq6MSaihVHakEhKRG4CfAG4Ee8s+ovqrw0aGSPgC8KiLeL+lK4K0R8Q5J5wF3AhcDc4FvA+dExNDo4xRyQFSXQwNDPNG7h0eeyV+W6t66i71JP8ZoZ7Q0smLpPFoaG2huFM2NDTQ3NtDSNGq5sYHmplHLI/s0vbg8st/R9SOf09BAQ4O4Z31v0U77yQ4q11F5dVRCDaWsI62AeC3w8Yi4LFm+ESAibi7YZnWyzQ8lNQE/A9qBGwq3LdzuRMd0QFS34eHg5R+577gb9ka0z5jCwNAwRwaHGRgaZmCoPP/tNjWIoeEoWkeDYM6MqYwM0irseh89cuvoNgXNSvYYvX/hvip4s+35AwwOH19JU4PonD1tQv+eUtjy3P4x61g4iXU8UwF1VEINJ6pj3qxWHrrhVyf8OWk9cnQesL1guQd49VjbRMSgpD3AS5L2h0ftWzQSJV0LXAuwYMGCkhRu6Who0Jj3YxT7jz4iGBiKJCyGOZKExsDgqOWhYQYGX1w+csz64WT7eHF5aJjbvvPTojUOB7z+nNnJ8QtqOVrTyHIcu+KYbaLoPsW22dy3v2gdg8PB4rNmFF1XDpt2Fn+w1OBwsOis6ZNWx9MVUEcl1HCiOsa66fVUVP19EBGxElgJ+TOIlMux0zTW/RjXX7b4uG0l0dIkWppKP2PMPet3jBlUf/W2C0t+vLGs3/bAmHXc9s5lk1bHo7eMXcc/vrNr0uq4pALqqIQaTlTH3FmtJTtGOedi6gU6CpbnJ21Ft0kuMc0k31k9kX2tBq1YOo+br7iAebNaEfn/6dLooL7+ssW0jhqKO1ZQuY76qqMSapisOsp5BrEGWCRpIflf7lcCvz1qm1XANcAPgbcBD0RESFoFfFnS35DvpF4EPFLGWq2CrFg6L/VhrSPHT3ukiuuovDoqoYbJqqPcw1x/A/gU+WGut0fEX0q6CeiOiFWSpgL/CiwF+oErI2Jzsu+fA+8BBoE/iohvjHc8d1KbmZ0c3yhnZmZFnSgg/DwIMzMrygFhZmZFOSDMzKwoB4SZmRVVU53UkvqArae4+2zguRKWU838XRzL38ex/H28qBa+i2xEtBdbUVMBcTokdY/Vk19v/F0cy9/Hsfx9vKjWvwtfYjIzs6IcEGZmVpQD4kUr0y6ggvi7OJa/j2P5+3hRTX8X7oMwM7OifAZhZmZFOSDMzKyoug8IScslbZS0SdINadeTJkkdkr4j6UlJGyT9Ydo1pU1So6T1kv4z7VrSJmmWpLsk/Y+kp5LHCtctSX+c/H/yY0l3JrNT15S6DghJjcBtwJuA84CrJJ2XblWpGgT+NCLOA14DfLDOvw+APwSeSruICvFp4JsRcS5wIXX8vUiaB3wIyEXE+eQfaXBlulWVXl0HBHAxsCkiNkfEEeArwOUp15SaiHg2ItYl7/eS/wWQ7pN7UiRpPvC/gM+nXUvaJM0EXg98ASAijkTE7nSrSl0T0Jo8DfMMYEfK9ZRcvQfEPGB7wXIPdfwLsZCkTvIPcvpRupWk6lPAh4HhtAupAAuBPuCfk0tun5c0Le2i0hIRvcBfA9uAZ4E9EfGtdKsqvXoPCCtC0nTgq+Sf5PdC2vWkQdKbgZ0RsTbtWipEE7AM+ExELAX2A3XbZycpQ/5qw0Lyj0WeJuld6VZVevUeEL1AR8Hy/KStbklqJh8OX4qIu9OuJ0WXAG+RtIX8pcdflfR/0y0pVT1AT0SMnFHeRT4w6tWvA89ERF9EDAB3A7+Yck0lV+8BsQZYJGmhpBbynUyrUq4pNZJE/hrzUxHxN2nXk6aIuDEi5kdEJ/n/Lh6IiJr7C3GiIuJnwHZJi5OmXwOeTLGktG0DXiPpjOT/m1+jBjvtm9IuIE0RMSjpOmA1+VEIt0fEhpTLStMlwNXAE5IeTdo+EhH3pViTVY4/AL6U/DG1GfjdlOtJTUT8SNJdwDryo//WU4PTbniqDTMzK6reLzGZmdkYHBBmZlaUA8LMzIpyQJiZWVEOCDMzK8oBYZaQtC/52Snpt0v82R8ZtfyDUn6+WTk4IMyO1wmcVEAkE7adyDEBERE1d9et1R4HhNnxbgF+SdKjyZz/jZJulbRG0uOSfg9A0hskfU/SKpK7iiXdI2lt8pyAa5O2W64vF/oAAAHdSURBVMjP+vmopC8lbSNnK0o++8eSnpD0joLP/m7B8xe+lNyxi6Rbkmd2PC7pryf927G6Udd3UpuN4QbgzyLizQDJL/o9EXGRpCnAQ5JGZu5cBpwfEc8ky++JiH5JrcAaSV+NiBskXRcRS4oc6wpgCfnnK8xO9nkwWbcUeCX5aaQfAi6R9BTwVuDciAhJs0r+rzdL+AzCbHyXAu9Oph/5EfASYFGy7pGCcAD4kKTHgIfJTwS5iBN7HXBnRAxFxM+B/wYuKvjsnogYBh4lf+lrD3AI+IKkK4ADp/2vMxuDA8JsfAL+ICKWJK+FBXP/7z+6kfQG8rN8vjYiLiQ/P8/pPIbycMH7IaApIgbJP+jqLuDNwDdP4/PNTsgBYXa8vcCMguXVwO8nU6Ej6ZwxHpYzE9gVEQcknUv+sa0jBkb2H+V7wDuSfo528k9te2SswpJndcxMJlD8Y/KXpszKwn0QZsd7HBhKLhXdQf5ZzJ3AuqSjuA9YUWS/bwLvT/oJNpK/zDRiJfC4pHUR8c6C9q8BrwUeAwL4cET8LAmYYmYA/yFpKvkzmz85tX+i2fg8m6uZmRXlS0xmZlaUA8LMzIpyQJiZWVEOCDMzK8oBYWZmRTkgzMysKAeEmZkV9f8Bx7c5v29RgPIAAAAASUVORK5CYII=\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def gaussseidel_plot():\n", + " \n", + " iterations = np.array([])\n", + " errors = np.array([])\n", + " \n", + " A = np.matrix(\"4 -1 -1; -2 6 1; -1 1 7\")\n", + " b = np.matrix(\"3; 9; -6\")\n", + " trueresult = np.linalg.solve(A, b)\n", + "\n", + " for n in range(32):\n", + " \n", + " result = gaussseidel(A.copy(), b, n)\n", + " \n", + " error = np.abs(np.sum(result-trueresult))\n", + " \n", + " iterations = np.append(iterations, n)\n", + " errors = np.append(errors, error)\n", + " \n", + " if np.isclose(error, 0):\n", + " print(\"Using:\")\n", + " print(f\"A=\\n{A}\")\n", + " print(f\"b=\\n{b}\")\n", + " print(f\"Converged in {n}:\")\n", + " print(f\"x=\\n{result}\")\n", + " print(f\"|error|=\\n{error}\")\n", + " break\n", + "\n", + " plt.figure()\n", + " plt.xlabel('Iterations')\n", + " plt.ylabel('Error')\n", + " plt.plot(iterations, errors, 'o-')\n", + "\n", + "gaussseidel_plot()" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "We find that Gauss-Seidel iteration yields an approximate solution in 18 iterations where the solution is very close to zero. It is important to note the improvement in convergence speed, compared to the Jacobi iteration method, Gauss-Seidel is twice as fast!" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Newton's method, f(x)=0 ###\n", + "\n", + "Below are tests for Newton's method. Input function f(x) is specified such that there exists an exact solution y for f(x)=0. This means that the test explores convergence ||f(x)|| and the approximation error ||x-y|| at the same time." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 7, + "outputs": [ + { + "name": "stdout", + "text": [ + "Using:\n", + "fn=\n", + "lambda\n", + "Converged in 12:\n", + "x=\n", + "-2.0000000062094836\n", + "|error|=\n", + "6.2094835939774384e-09\n" + ], + "output_type": "stream" + }, + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3xV5Z3v8c8vOyEJkAtIghBAQJBwEcHipV4YxVoca1vqtFMv7bTjnMEzvbcz9CWdP+Z0OqPO2Os57XSK1tq+am2dVq11rIwVK2q9cRMIFy+AQAImICQBkpBk/84fewdjzGUn2Wuvffm+X6/9yt5r76znt7h8s/KsZz2PuTsiIpJ98sIuQEREgqGAFxHJUgp4EZEspYAXEclSCngRkSyVH3YB3Y0bN86nTp0adhkiIhlj/fr1h9y9orf30irgp06dyrp168IuQ0QkY5jZG329py4aEZEspYAXEclSCngRkSylgBcRyVIKeBGRLBXoKBozKwfuAuYBDtzk7s8F2WYQHtpYyx2rd1J3tIWJ5cWsWDqLZQurwi5LRKRfQQ+T/B7wmLt/1MxGACMDbi/pHtpYy8oHttDS3glA7dEWVj6wBUAhLyJpLbAuGjMrAxYDPwZw95PufjSo9oJyx+qdp8K9S0t7J3es3hlSRSIiiQmyD34a0AD8xMw2mtldZjaq54fMbLmZrTOzdQ0NDQGWMzR1R1sGtV1EJF0EGfD5wLnAD919IXAcuKXnh9x9lbsvcvdFFRW93m0bqonlxYPaLiKSLoIM+P3Afnd/If7618QCP6OsWDqLooJ3/jEVF0RYsXRWSBWJiCQmsIB394PAPjPrSsIrgG1BtReUZQurWL54+qnXeQa3fmSeLrCKSNoLehz854F7zWwzsAC4NeD2AjFm5AggdjYfdZg5viTkikREBhZowLv7pnj/+nx3X+buR4JsLyg1dU2MG13IXy6ajBms2VEfdkkiIgPSnawJqKlrYu7EUipKCpk/qVwBLyIZQQE/gLaOTl59s5m5E0sBuKK6kpf3H+XQsbaQKxMR6Z8CfgCvHDxGR9SZV1UGwJLqStzhjzvTb8y+iEh3CvgB1NQ1Apw6g587sZTxpYU8qW4aEUlzCvgB1NQ1UVKYz+QxsWl0zIzLZ1Wy9pUG2jujIVcnItI3BfwAttY1MmdiKXl5dmrbkupKmts6eGnPWyFWJiLSPwV8Pzqjzo4DzcydWPaO7RfPGMeISB5rtqubRkTSlwK+H7sPHaOlvfNU/3uXUYX5XDB9LGt2KuBFJH0p4PtRU9cEwNyq0ne9d0V1JbsajrPn0PFUlyUikhAFfD+21jZSmJ/HjIrR73pvSfV4QHe1ikj6UsD3o6auierTS8iPvPuPacppI5lROVoBLyJpSwHfB3enpq6JOT0usHa3pLqSF3Yf5lhbRworExFJjAK+D7VHW2hsaX/XBdbullRX0t7pPPOq7moVkfSjgO/D1trYBdauKQp6854zxlBSlK9uGhFJSwr4PmyraySSZ1Sf3vfc7wWRPP7srArW7GggGvUUViciMjAFfB9q6po4s2IURQWRfj+3pLqSQ8fa2Bqfs0ZEJF0o4Puwta6Ref1cYO1y2axKzOAJ3dUqImlGAd+LQ8faeLOpjTn9XGDtMnbUCBZOLudJ3dUqImlGAd+LU3ewJnAGD3DF7PFs3t9IfVNrkGWJiAyKAr4XXXPAJ3IGD3D5rEpAi4CISHpRwPeipraJKWNHUlZckNDnZ08oYUJZEU/seDPgykREEqeA70VNXWO/Nzj1ZGZcXl3JM68eoq2jM8DKREQSp4Dvobm1nT2HTwwq4CE2u+Txk528uFuLgIhIelDA97D9QDOQ+AXWLhedOY7C/Dzd1SoiaSPQgDezPWa2xcw2mdm6INtKlq218UW2e5kDvj/FIyJcdOZprNlRj7vuahWR8KXiDP5yd1/g7otS0Naw1dQ1UVFSSGVJ0aC/d0l1JW8cPsEuLQIiImlAXTQ9DPYCa3eXV8eGSz6pbhoRSQNBB7wD/2Nm681seW8fMLPlZrbOzNY1NIQ7jry1vZNX648NOeAnjRnJrPElmrZARNJC0AF/ibufC/w58FkzW9zzA+6+yt0XufuiioqKgMvp3ytvNtMZ9YTmoOnL5dWVvLTnLZpa25NYmYjI4AUa8O5eG/9aDzwInB9ke8M12CkKenPF7Eo6os7TrxxKVlkiIkMSWMCb2SgzK+l6Drwf2BpUe8lQU9dISVE+k8cWD3kfCyeXU1ZcoOGSIhK6/AD3PR540My62vmFuz8WYHvDtrW2ibkTS4nXPCT5kTwum1XBH3fWE406eXlD35eIyHAEdgbv7rvc/Zz4Y667/2tQbSVDZ9TZcbBpWN0zXZZUV3L4+Ele3n80CZWJiAyNhknG7Wo4Rmt7dMgjaLr7s7MqyDPUTSMioVLAx3UtuZeMM/jykSN4zxljFPAiEioFfFxNbROF+XmcWTEqKftbUj2emromDjZqERARCYcCPq6mronqCaXkR5LzR7Kk665WLeUnIiFRwAPuPqwpCnpz1vjRVJUX665WEQmNAh7Yf6SFptaOpAa8mbGkupJnXztEa7sWARGR1FPA8/YarMOZoqA3S2ZX0tLeyfO7Did1vyIiiVDAE+t/j+QZs04vSep+3zv9NIoK8jS7pIiEQgFPLOBnVIymqCCS1P0WFUS4ZMY4ntAiICISAgU8sVWcBruCU6Iur65k/5EWXqs/Fsj+RUT6kvMB39DcRn1zW1JucOpN13DJJ9RNIyIplvMBX3PqDtZgzuAnlBUze0Kp7moVkZRTwMfngJ8TUMADXFFdyfo3jtB4QouAiEjqKODrGjnjtJGUFhUE1sbl1ZV0Rp2nXg13SUIRyS0K+LqmwLpnuiyYXM7YUSNYs/3NQNsREekupwO+qbWdNw6fCOwCa5dInnHZWRU89UoDnVENlxSR1MjpgN+Wgv73LpdXV3LkRDub9h0JvC0REcjxgO+6wJrsKQp6s/isCiJ5psnHRCRlcjzgG6ksKaSipDDwtsqKC1ikRUBEJIVyOuC3peACa3dXzK5kx8Fmao+2pKxNEcldORvwre2dvFp/jHlVwXfPdDm1CIjO4kUkBXI24HcebKYz6ik9gz+zYjRTxo5UN42IpETOBnzXBdagh0h2130RkJaTWgRERIKVswG/ta6R0qJ8Jo0pTmm7S6oraeuI8tyuQyltV0RyT+ABb2YRM9toZo8E3dZgxO5gLcPMUtruBdPHMnJERN00IhK4VJzBfxHYnoJ2EtbRGWXHgdSOoOlSmB9bBGTNdi0CIiLBCjTgzWwS8AHgriDbGaxdh47T1hENbJGPgSyprqSusZWdbzaH0r6I5Iagz+C/C3wViPb1ATNbbmbrzGxdQ0NqZlvcWts1B3zqLrB2d3nXIiC6q1VEAhRYwJvZNUC9u6/v73PuvsrdF7n7ooqKiqDKeYeauiaKCvKYPm5UStrraXxpEfOqSjUeXkQCFeQZ/MXAh8xsD/BLYImZ/TzA9hJWU9dI9eml5EfCG0S0pHo8G/Ye4cjxk6HVICLZLbCEc/eV7j7J3acC1wFr3P0TQbWXKHdPyRzwA1lSXUnU4alXtAiIiAQj58bB73urhebWjpROUdCb+VVljBs9Qotxi0hgUhLw7v5Hd78mFW0NJOhFthOVl2dcNquSp3bW09HZ5zVoEZEhy7kz+Jq6JiJ5xlnjS8IuhSXVlTS1drD+DS0CIiLJl3MBv7WukZmVoykqiIRdCpfOHEd+nrFmp7ppRCT5ci7gu6YoSAclRQWcP20sazQeXkQCkFMBX9/cSkNzW+j9790tqa7k1fpj7HvrRNiliEiWyamAf3uK4PQKeECTj4lI0uVWwMenKJiTRgE/vWI008aNUsCLSNLlVsDXNTH1tJGUFBWEXco7XD6rkud2HebEyY6wSxGRLJJzAZ8uF1i7u2J2JSc7ojz72uGwSxGRLJIzAd/Y0s7et06kVfdMl/OmjmV0Yb66aUQkqXIm4Lel4QXWLiPy87h05jie3KFFQEQkeXIm4N+eoiD9umggNkf8waZWth1oCrsUEckS+WEXkCrb6poYX1pIRUlh2KX06vJZseGS1616nmOtHUwsL2bF0lksW1gVcmUikqlyJuC31jWm7dk7wLOvHcIMmltjI2lqj7aw8oEtAAp5ERmSnOiiaW3v5PWG48xLw/73Lnes3knP7veW9k7uWL0znIJEJOPlRMDvONhMZ9SZk8Zn8HVHWwa1XURkIDkR8OkyB3x/JpYXD2q7iMhAciLgt9Y2UVZcwKQx6RuWK5bOorjHFMbFBRFWLJ0VUkUikukGDHgzi5jZN1NRTFC21TUyd2IpZhZ2KX1atrCK2649m4nlRQDk5xm3fmSeLrCKyJANGPDu3glckoJaAtHeGWX7wea07p7psmxhFX+65Qr+/aPz6Yg6p41OzyGdIpIZEu2i2WhmD5vZJ83s2q5HoJUlyesNxzjZEU3rIZI9fXjBRCpLClm1dlfYpYhIBks04IuAw8AS4IPxR1osoj2Qmtr0naKgL4X5Ef764mk889ohtsanOBYRGayEbnRy978OupCg1NQ1UVSQx/SK0WGXMig3XDCF7695lTuf3sX3rlsYdjkikoESOoM3s0lm9qCZ1ccfvzGzSUEXlww1dY3MnlBKJC99L7D2pqy4gOvPn8Ijmw+w/4iW8xORwUu0i+YnwMPAxPjjd/FtaS0adbbVNWVU90x3N10yDQPufmZP2KWISAZKNOAr3P0n7t4Rf9wDVPT3DWZWZGYvmtnLZlZjZl8fdrWDtO/ICZrbOpiXQRdYu5tYXswHz5nIL1/aS+OJ9rDLEZEMk2jAHzazT8THxEfM7BPELrr2pw1Y4u7nAAuAq8zswuEUO1hvL7KdmQEP8LeXTufEyU5+/sIbYZciIhkm0YC/CfhL4CBwAPgo0O+FV485Fn9ZEH+kdDWLrbWN5OcZZ52eWRdYu5szsZRLZ47jnj/toa2jM+xyRCSDJHQnK3Ctu3/I3SvcvdLdl7n73kS+18w2AfXA4+7+Qi+fWW5m68xsXUNDw5AOoi81dU3MqBxNYX5k4A+nsZsXn0lDcxsPbawNuxQRySCJ3sl6/VB27u6d7r4AmAScb2bzevnMKndf5O6LKir67dYftJq6JuZVZW73TJeLZ5zGnAmlrFq7i2hUS/qJSGIS7aJ51sy+b2aXmtm5XY9EG3H3o8CTwFVDqnII6ptaOXSsLWNH0HRnZtz8Z9N5veG4FuYWkYQlGvALgLnAPwPfij/6nYDMzCrMrDz+vBi4Etgx9FIHZ2uar8E6WFefPYGq8mJNXyAiCRvwTlYzywN+6O73D3LfE4Cfxvvw84D73f2RIdQ4JF1TFMyeUJKqJgNVEMnjpkum8Y1HtrFx7xEWThkTdkkikuYS6YOPAl8d7I7dfbO7L3T3+e4+z93/eUgVDlFNXRPTxo2ipKgglc0G6rrzJlNalK+zeBFJSKJdNH8ws38ws8lmNrbrEWhlw1RzoJE5WdD/3t2ownw+ceEZPFZzkD2HjoddjoikuUQD/uPAZ4G1wPr4Y11QRQ1X44l29r3VkhUXWHv69EVTKcjL465ndBYvIv1LKODdfVovj+lBFzdUNQdiF1gzdYqC/lSWFvGRhVX817r9HD7WFnY5IpLG+g14M/tqt+cf6/HerUEVNVzb6jJvDvjB+NvF02jriPKz5zR9gYj0baAz+Ou6PV/Z472UjWkfrK21jZxeWpS1S97NqCzhfbMr+dlze2g5qekLRKR3AwW89fG8t9dpoyaDpwhO1PLFZ3LkRDu/Xr8v7FJEJE0NFPDex/PeXqeFlpOdvN5wjLlZMEVBf86bOoYFk8u565nddGr6AhHpxUABf46ZNZlZMzA//rzr9dkpqG/QdhxsIurZ2//excy4efF03jh8gtU1B8MuR0TSUL8B7+4Rdy919xJ3z48/73qdlncQbc3yC6zdvX/u6Uw9bSQ/WrsLd53Fi8g7JToOPmNsq2ukrLiAqvLisEsJXCTP+JtLp/PyvqO8tOdI2OWISJrJuoCPTRFcilnaXgNOqo+eO4mxo0awau3rYZciImkmqwK+vTPKjoPNWTODZCKKR0T45IVn8Ift9bxW3xx2OSKSRrIq4F+rP8bJjmhO9L9391fvPYPC/DzuXLs77FJEJI1kVcBnwyLbQ3Ha6EI+tmgSD26spb6pNexyRCRNZFnAN1JcEGHauFFhl5Jy/+uS6bRHo9zzpz1hlyIiaSK7Ar62idkTSojk5cYF1u6mjhvFVXNP5+fPv8Gxto6wyxGRNJA1AR+NOtsONOVc90x3yxdPp6m1g1+9pOkLRCSLAn7vWyc41tbBvKrcusDa3cIpYzh/6ljufmY37Z3RsMsRkZBlTcDn6gXWnpYvnk7t0RYe3XIg7FJEJGRZE/Bb6xrJzzNmjh8ddimhWlJdyZkVo/jRU5q+QCTXZU3A19Q1MXN8CYX5kbBLCVVenrF88XS2HWji2dcOh12OiIQoKwLe3ampbWRejt3g1JdlC6uoKCnkR5q+QCSnZUXAv9nUxuHjJ3PuDta+FOZH+PRFU3n61UOnli8UkdyTFQFfUxdbZDvbF/kYjE9ccAYjR0S48+ldYZciIiEJLODNbLKZPWlm28ysxsy+GEQ7D22s5Uu/2gTAF+7byEMba4NoJuOUjSzguvOm8LuX66g72hJ2OSISgiDP4DuAv3f3OcCFwGfNbE4yG3hoYy0rH9hCc2vszs0Dja2sfGCLQj7upkum4sDdz2gSMpFcFFjAu/sBd98Qf94MbAeqktnGHat30tLe+Y5tLe2d3LF6ZzKbyViTxozkmvkTuO/FvTS2tIddjoikWEr64M1sKrAQeKGX95ab2TozW9fQ0DCo/fbV9aAuibctXzyd4yc7+cULe8MuRURSLPCAN7PRwG+AL7n7u4Z0uPsqd1/k7osqKioGte+JfSzL19f2XDR3YhmXzBjHT57dTVtH58DfICJZI9CAN7MCYuF+r7s/kOz9r1g6i+KCd97YVFwQYcXSWcluKqMtXzyd+uY2frupLuxSRCSFghxFY8CPge3u/u0g2li2sIrbrj2bqvJiDKgqL+a2a89m2cKkdvVnvEtnjmP2hFLuXLuLaFTTF4jkivwA930x8Elgi5ltim/7mrs/msxGli2sUqAPwMxYvngaX/7Vy/zxlXqWVI8PuyQRSYEgR9E84+7m7vPdfUH8kdRwl8RdM38iE8uK+NFTuvFJJFdkxZ2sMrCCSB43XTKNF3a/xcv7joZdjoikgAI+h1x3/hQK843rVj3HtFv+m4tvX6ObwkSyWJB98JJm/rDtTTqj0BaNrfZUe7SFlQ9sAdB1DJEspDP4HHLH6p109BhFozt/RbKXAj6H6M5fkdyigM8huvNXJLco4HNIb3f+GvCFJTPCKUhEAqWAzyE97/wdN3oEDjy367AW6BbJQhpFk2N63vn7vT+8ynf+8AoXTD+N68+fEmJlIpJsOoPPcZ9bMoNLZ47jnx6uObX0oYhkBwV8jovkGd/5+ALGjCzgs/duoKlVC4OIZAsFvDBudCH/7/pz2XekhVt+s1n98SJZQgEvAJw/bSwrls7i0S0H+emf9oRdjogkgQJeTll+6XTeN7uSf310O5s0IZlIxlPAyyl5ecY3P3YOlSVFfPbeDRw9cTLskkRkGBTw8g7lI0fwgxvPpb65lb+//2WtACWSwRTw8i4LJpfzj1fP5okd9ax6WguEiGQqBbz06lMXTeXqs0/njtU7eXH3W2GXIyJDoICXXpkZt//FfCaPKebz923g0LG2sEsSkUFSwEufSosK+MGN53LkRDtf+uUmOtUfL5JRFPDSr7kTy/j6h+byzGuH+P6a18IuR0QGQQEvA7ruvMl8ZGEV333iFZ597VDY5YhIghTwMiAz41+WzePMitF88ZcbebOpNeySRCQBCnhJyKjCfH5447kcb+vk8/dtpKMzGnZJIjKAwALezO42s3oz2xpUG5JaM8eXcOu183hx91t8+/FXwi5HRAYQ5Bn8PcBVAe5fQvCRhZO4/vzJ/McfX+fJHfVhlyMi/Qgs4N19LaA7ZLLQP31wLrMnlPLl+zdRe7Ql7HJEpA+h98Gb2XIzW2dm6xoaGsIuRxJQVBDhP248l45O53O/2MDJDvXHi6Sj0APe3Ve5+yJ3X1RRURF2OZKgaeNG8W9/MZ+Ne49y++93hF2OiPQi9ICXzPWB+RP49EVTufvZ3Ty29UDY5YhIDwp4GZaVV1dzzqQyVvzXZt44fDzsckSkmyCHSd4HPAfMMrP9ZvY3QbUl4SnMj/D9G87FDD5z7wZa2zvDLklE4oIcRXO9u09w9wJ3n+TuPw6qLQnX5LEj+dZfLqCmrolvPLIt7HJEJE5dNJIUV84Zz82Lp3PvC3v57abasMsREcDc02cK2EWLFvm6devCLkOGqL0zyg13Ps+mfUcZM3IEDc1tTCwvZsXSWSxbWBV2eSJZyczWu/ui3t7TGbwkTUEkjw/Mn0B7p1Pf3IYDtUdbWPnAFh7aqLN6kVRTwEtS3bl297u2tbR3csfqnSFUI5LbFPCSVHV9TF3Q13YRCY4CXpJqYnlxr9tHFeZrCKVIiingJalWLJ1FcUHkHdsiecaxtg4+8H+fZv0bmn9OJFUU8JJUyxZWcdu1Z1NVXowBVeXFfOtj5/Czm86ntT3KR//zOb7+uxpOnOwIu1SRrKdhkpIyx9o6+PfHdvCz595g8thi/u3a+Vw0Y1zYZYlkNA2TlLQwujCff/7wPH61/EIiZtxw1wusfGALTa3tYZcmkpUU8JJyF0w/jce+tJibF0/nVy/t5f3fXsuaHW+GXZZI1lHASyiKCiKsvHo2D3zmYkqL87npnnV85VebOHL8ZNiliWQNBbyEasHkcn73+Uv4whUzefjlOq78zlM8ukVzy4skgwJeQleYH+ErV57Fw5+7hNPLivjMvRv4u5+vp765NezSRDKaAl7SxpyJpTz0mYv56lWzeGJHPVd+ey0PbNhPOo30EskkCnhJK/mRPD5z2Qwe/cKlzKgczVfuf5mb7nlJUx2IDIECXtLSjMrR3H/ze/mnD87h+V1v8f7vrOUXL+zV2bzIICjgJW1F8oy/vngaq7+0mPmTyvjag1u44c4XtParSIJ0J6tkBHfnly/t49b/3k57NMqKpdWMKS7gW4+/Qt3RFi0sIjmrvztZ81NdjMhQmBnXnz+Fy2ZV8I8PbuUbj2zDDLrOT7oWFgEU8iJx6qKRjDKhrJgff2oRY0YW0POXz5b2Tm77/Xb104vE6QxeMo6ZcfRE7/PXvNnUxnv+5Q+cXVXG/EllzJ9UzvxJZYwvLUpxlSLhU8BLRppYXkxtL0Mny4oLeN/sSjbvb+QHTzYQjZ/Mjy8tjIV9VRnzJ8e+jhk1IsVVi6SWAl4y0oqls1j5wBZauq0SVVwQ4esfmnuqD77lZCc1dY1s3t/I5v1H2by/kce3vT2p2eSxxW+H/qRy5lWVUlJU8I52HtpYyx2rd+pCrmSkQAPezK4CvgdEgLvc/fYg25Pc0RWy/YVv8YgIi6aOZdHUsae2NbW2s3V/I5trY6G/ae9R/ntzbO4bM5g+btSpbp2jJ9r50drXaW2PAsFdyE3VD5FUtJMtbaSqnaDbCGyYpJlFgFeAK4H9wEvA9e6+ra/v0TBJCcPhY21srm1kS/xM/+X9jTQ0t/X5+ZEjIlx33hQKC/IozM+jMD9CYX4eI/Ljrwsi8e3x93p8LvY69vz3Ww7wtQe3vus3kduuPTvpP0R6+40nme1kSxupaidZbfQ3TDLIgH8v8H/cfWn89UoAd7+tr+9RwEs6cHfebGrjwtue6PMzowvzaevopL0zmP8/eQanjS7EiP1mAWBYt+exi81dzPr/3L63TtARfXet+XnGGaeNTErNbxzuu42p40YlpY09h4732ca0JLUBsDsF7fTVRlV5Mc/esiTh/YQ1Dr4K2Nft9X7ggp4fMrPlwHKAKVOmBFiOSGLMjNPLiqjq40Ju9/+A0ahzsjNKW3uUto5O2jpiX1vbo6eet3W8/f7Jjq7tsdf//tjOXmuIOrxvduWpoaDu4Hi35z22d30OTg0T9W6f3X2o97t/O6JO9YTSofwxvcvrDX23MWt8SVLaeK3+WJ9tzBw/OiltALyagnb6aiOZ8y6FfpHV3VcBqyB2Bh9yOSKn9HUhd8XSWade5+UZRXkRigoiQEEve+nfvc/v7fOHyG3Xzh9S3b3Z8MaRPtv5wQ3nJqWNTXvX9N3GjUlq4/a+2/iPG9+TlDYALk5BO321MbG8OCn7h2BvdKoFJnd7PSm+TSQjLFtYxW3Xnk1VeTFGV+gmt693xdJZFBdE3rGt5w+RTGknW9pIVTupaCPIM/iXgJlmNo1YsF8H3BBgeyJJt2xhVaDDIhMZDZQp7WRLG6lqJxVtBDrZmJldDXyX2DDJu939X/v7vC6yiogMTmiTjbn7o8CjQbYhIiK902RjIiJZSgEvIpKlFPAiIllKAS8ikqXSask+M2sA3hjit48DDiWxnDBly7Fky3GAjiUdZctxwPCO5Qx3r+jtjbQK+OEws3V9DRXKNNlyLNlyHKBjSUfZchwQ3LGoi0ZEJEsp4EVEslQ2BfyqsAtIomw5lmw5DtCxpKNsOQ4I6Fiypg9eRETeKZvO4EVEpBsFvIhIlsr4gDezq8xsp5m9Zma3hF3PUJnZZDN70sy2mVmNmX0x7JqGy8wiZrbRzB4Ju5bhMLNyM/u1me0ws+3x5Sgzjpl9Of5va6uZ3WdmRWHXlCgzu9vM6s1sa7dtY83scTN7Nf51TJg1JqqPY7kj/u9rs5k9aGblyWgrowM+vrD3D4A/B+YA15vZnHCrGrIO4O/dfQ5wIfDZDD6WLl8EtoddRBJ8D3jM3auBc8jAYzKzKuALwCJ3n0dsCu/rwq1qUO4Bruqx7RbgCXefCTwRf50J7uHdx/I4MM/d5wOvACuT0VBGBzxwPvCau+9y95PAL4EPh1zTkLj7AXffEH/eTCxEgltpImBmNgn4AHBX2LUMh5mVAYuBHwO4+0l3PxpuVUOWDxSbWT4wEqgLuZ6Eufta4K0emz8M/DT+/KfAspKVZUkAAARESURBVJQWNUS9HYu7/4+7d8RfPk9sBbxhy/SA721h74wNxS5mNhVYCLwQbiXD8l3gq0A07EKGaRrQAPwk3t10l5mNCruowXL3WuCbwF7gANDo7v8TblXDNt7dD8SfHwTGh1lMEt0E/D4ZO8r0gM86ZjYa+A3wJXdvCrueoTCza4B6d18fdi1JkA+cC/zQ3RcCx8mcroBT4v3THyb2A2siMMrMPhFuVcnjsfHeGT/m28z+kVh37b3J2F+mB3xWLextZgXEwv1ed38g7HqG4WLgQ2a2h1i32RIz+3m4JQ3ZfmC/u3f9NvVrYoGfad4H7Hb3BndvBx4ALgq5puF608wmAMS/1odcz7CY2aeBa4AbPUk3KGV6wJ9a2NvMRhC7aPRwyDUNiZkZsX7e7e7+7bDrGQ53X+nuk9x9KrG/kzXunpFni+5+ENhnZl1L3V8BbAuxpKHaC1xoZiPj/9auIAMvFvfwMPCp+PNPAb8NsZZhMbOriHVpfsjdTyRrvxkd8PGLEp8DVhP7x3q/u9eEW9WQXQx8ktjZ7qb44+qwixIAPg/ca2abgQXArSHXM2jx30B+DWwAthD7v58xt/qb2X3Ac8AsM9tvZn8D3A5caWavEvsN5fYwa0xUH8fyfaAEeDz+f/8/k9KWpioQEclOGX0GLyIifVPAi4hkKQW8iEiWUsCLiGQpBbyISJZSwEvWMLNj8a9TzeyGJO/7az1e/ymZ+xcJggJestFUYFABH5+Aqz/vCHh3z/S7QCUHKOAlG90OXBq/YeTL8Xnp7zCzl+Lzbd8MYGaXmdnTZvYw8btTzewhM1sfnzd9eXzb7cRmYdxkZvfGt3X9tmDxfW81sy1m9vFu+/5jt3nk743fQYqZ3R6f93+zmX0z5X86kjMGOmsRyUS3AP/g7tcAxIO60d3PM7NC4Fkz65pJ8Vxi83Dvjr++yd3fMrNi4CUz+42732Jmn3P3Bb20dS2xu1vPAcbFv2dt/L2FwFxi0/I+C1xsZtuBjwDV7u7JWthBpDc6g5dc8H7gr8xsE7EpmE8DZsbfe7FbuAN8wcxeJjYn9+Run+vLJcB97t7p7m8CTwHnddv3fnePApuIdR01Aq3Aj83sWiBp846I9KSAl1xgwOfdfUH8Ma3bXOjHT33I7DJic5q8193PATYCw1nWrq3b804gPz5/0vnE5oW5BnhsGPsX6ZcCXrJRM7GJm7qsBv4uPh0zZnZWH4t2lAFH3P2EmVUTWzqxS3vX9/fwNPDxeD9/BbHVn17sq7D4fP9l7v4o8GViXTsigVAfvGSjzUBnvKvlHmJrqk4FNsQvdDbQ+/JujwH/O95PvpNYN02XVcBmM9vg7jd22/4g8F7gZWILTnzV3Q/Gf0D0pgT4rcUWvDbgK0M7RJGBaTZJEZEspS4aEZEspYAXEclSCngRkSylgBcRyVIKeBGRLKWAFxHJUgp4EZEs9f8Bri2/z7W4zikAAAAASUVORK5CYII=\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def newtonsmethod_plot():\n", + " \n", + " iterations = np.array([])\n", + " errors = np.array([])\n", + "\n", + " x0 = 0.0\n", + " fn = lambda x : x**3+8\n", + " trueresult = optimize.root_scalar(fn, bracket=[-10, 10]).root\n", + "\n", + " for n in range(128):\n", + " \n", + " result = newtonsmethod(fn, x0, n)\n", + " \n", + " error = np.abs(result-trueresult)\n", + " \n", + " iterations = np.append(iterations, n)\n", + " errors = np.append(errors, error)\n", + " \n", + " if np.isclose(error, 0):\n", + " print(\"Using:\")\n", + " print(f\"fn=\\nlambda\")\n", + " print(f\"Converged in {n}:\")\n", + " print(f\"x=\\n{result}\")\n", + " print(f\"|error|=\\n{error}\")\n", + " break\n", + "\n", + " plt.figure()\n", + " plt.xlabel('Iterations')\n", + " plt.ylabel('Error')\n", + " plt.plot(iterations, errors, 'o-')\n", + "\n", + "newtonsmethod_plot()" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "We find that Newton's method yields an approximate solution in 12 iterations where the solution is very close to zero. It is interesting to see that the approximation error actually increases in the beginning and then rapidly decreases." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Discussion**\n", + "\n", + "This lab showed that iterative methods can be highly effective in solving equations. They are definitely an alternative to consider in comparison with the direct methods previously explored." + ], + "metadata": { + "collapsed": false + } + } + ], + "metadata": { + "kernelspec": { + "name": "python3", + "language": "python", + "display_name": "Python 3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + }, + "pycharm": { + "stem_cell": { + "cell_type": "raw", + "source": [], + "metadata": { + "collapsed": false + } + } + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/Lab-4/lab4_approximation.ipynb b/Lab-4/lab4_approximation.ipynb new file mode 100644 index 0000000..eeba629 --- /dev/null +++ b/Lab-4/lab4_approximation.ipynb @@ -0,0 +1,230 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "[](https://colab.research.google.com/github/johanhoffman/DD2363-VT20/blob/tree/viktorme/Lab-4/lab4_approximation.ipynb)\n", + "# Lab 4: Approximation #\n", + "**Viktor Meyer - DD2363 Methods in Scientific Computing**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Abstract**\n", + "This lab is an exercise in approximation methods. The mandatory part includes L2 projection to piece-wise linear approximation over mesh in 1D. There is also an extra assignment L2 projection to piece-wise linear approximation over triangular mesh in 2D." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Environment**\n", + "To have access to the neccessary modules you have to run this cell." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 208, + "outputs": [], + "source": [ + "# Load neccessary modules.\n", + "\n", + "import numpy as np\n", + "import scipy.integrate as integrate\n", + "import math\n", + "\n", + "from matplotlib import pyplot as plt" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Introduction**\n", + "The methods to be implemented were covered in the during lectures and definitions are available in the lecture notes. The aim of this lab is to approach the problem of solving systems of linear equations using approximation. This may be of use when e.g solutions do not exist and we seek a best possible solution or are constrained by computational work/memory footprint. \n", + "> LN-DD2363-part5.pdf, p. 161" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Methods**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "### L2 projection ###\n", + "\n", + "> LN-DD2363-part5.pdf, pp. 167-171" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 209, + "outputs": [], + "source": [ + "def l2(f, x):\n", + " \n", + " n = x.size\n", + " a = np.zeros(shape=(n,n))\n", + " b = np.zeros(n)\n", + " \n", + " for i in range(1, n):\n", + " \n", + " # build mass matrix\n", + " \n", + " h = x[i] - x[i - 1]\n", + " \n", + " # a_{i, i} = h_{i}/3\n", + " a[i, i] += h/3\n", + " \n", + " # a_{i, i} = h_{i+1}/3 \n", + " a[i - 1, i - 1] += h/3\n", + " \n", + " # a_{i, i+1} = h_{i+1}/6 \n", + " a[i - 1, i] += h/6\n", + " \n", + " # a_{i, i-1} = h_{i}/6 \n", + " a[i, i - 1] += h/6\n", + " \n", + " # build load vector\n", + " \n", + " l = x[i-1]\n", + " c = x[i]\n", + " r = x[i+1 if i+1 < n else i]\n", + " b[i] += integrate.quad(lambda x: (f(x)*((x-l)/h)), l, c)[0]\n", + " b[i] += integrate.quad(lambda x: (f(x)*((r-x)/h)), c, r)[0]\n", + " \n", + " return np.linalg.solve(a,b)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Results**\n", + "\n", + "### Convergence and Accuracy###\n", + "We check that the error converges towards zero as the number of nodes increases. This test also ensures that accuracy is sufficient since we evaluate the real function value at each point and compare it to the approximation" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 210, + "outputs": [ + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAeJUlEQVR4nO3deZSddZ3n8ff37rVnqYWsJISQGEG2miAqGDKOB7SFllabOLh0K6hHZvTonNZeju3Q3TPtMjqnbUYF0XZ0FFGxJ7bYoBKkYSBQYTVAsAgBErJU1tqXW/c7f9ynKjdFVVKp1FPPvfV8Xufcc5/t3vt9zoX65Pk99/f7mbsjIiLxlYi6ABERiZaCQEQk5hQEIiIxpyAQEYk5BYGISMyloi7gZDU2NvqyZcuiLkNEpKJs2bJlv7s3jbev4oJg2bJltLW1RV2GiEhFMbMXJ9qnpiERkZhTEIiIxJyCQEQk5hQEIiIxpyAQEYk5BYGISMwpCEREYi42QdC24yBf+Ndn0bDbIiLHik0QPLnzCF+/93kO9w5FXYqISFmJTRC01OcA2NvVH3ElIiLlJUZBkAVgb+dAxJWIiJSXGAVBcEXQqSsCEZFSsQmCprriFcE+BYGIyDFiEwS5dJI51Wk1DYmIjBGbIABoqcupaUhEZIxYBUFzfZa9XboiEBEpFasgaKnP6R6BiMgYMQuCLPu6BigU1LtYRGREzIIgx3DBOdAzGHUpIiJlI1ZB0FynvgQiImPFKghGehfv0zATIiKjYhYEI1cE+uWQiMiIWAXBSO9iNQ2JiBwVqyBIJxM01mZ0RSAiUiLUIDCzy81sm5m1m9lnj3PcH5mZm1lrmPVA8Yax+hKIiBwVWhCYWRK4CbgCWANsMLM14xxXB3wC2BxWLaVa6rOak0BEpESYVwRrgXZ33+7ug8BtwFXjHPc3wBeAGfnr3FKfU9OQiEiJMINgEfByyfrOYNsoM7sAWOLuvzjeG5nZ9WbWZmZtHR0dp1RUc32O/d0D5IcLp/Q+IiKzRWQ3i80sAXwF+PSJjnX3m9291d1bm5qaTulzW+qzuMP+bvUuFhGBcINgF7CkZH1xsG1EHXA2cK+Z7QBeD2wM+4Zxi3oXi4gcI8wgeARYaWbLzSwDXANsHNnp7kfcvdHdl7n7MuAh4Ep3bwuxJk1ZKSIyRmhB4O554AbgLuAZ4HZ332pmN5rZlWF97omMTmKveQlERABIhfnm7n4ncOeYbZ+b4Nh1YdYyYn5tloRp7mIRkRGx6lkMkEwYTXVZNQ2JiARiFwSgvgQiIqViGQTNmsReRGRULINgZMpKERGJbRDkONgzyEB+OOpSREQiF9MgCGYq030CEZF4BkFz0KlMU1aKiMQ0CI4OM6ErAhGReAZBvaasFBEZEcsgmFudIZ00XRGIiBDTIEgkTFNWiogEYhkEAM2aslJEBIhxELTUaZgJERGIcRCc1qBhJkREIMZB0Fyfpas/T+9gPupSREQiFdsgGOlLoN7FIhJ38Q0CTVkpIgLEOgg0ZaWICMQ4CEbHG9IVgYjEXGyDoD6XIpdOqGlIRGIvtkFgZpqyUkSEGAcBjHQq0xWBiMRbrIOgWVNWiojEOwiKTUP9uHvUpYiIRCbmQZCld3CY7gH1LhaR+Ip5EGimMhGRWAdBc536EoiIxDoIjvYuVhCISHzFOgia1TQkIhLvIKjNpqjNptSXQERiLdZBAEFfAl0RiEiMxT4I1LtYROJOQaBJ7EUk5hQEwcBz6l0sInEV+yBors8xmC9wpG8o6lJERCIR+yAY6UuwR/cJRCSmFATqSyAiMacgqNMk9iISb6EGgZldbmbbzKzdzD47zv6PmtlTZva4md1vZmvCrGc8zUHTkMYbEpG4Ci0IzCwJ3ARcAawBNozzh/4H7n6Ou58HfBH4Slj1TCSXTtJQlVbTkIjEVphXBGuBdnff7u6DwG3AVaUHuHtnyWoNEMlvOFvqs2oaEpHYSoX43ouAl0vWdwIXjT3IzD4OfArIAOvHeyMzux64HmDp0qXTXmhLfY69mrJSRGIq8pvF7n6Tu68APgP81QTH3Ozure7e2tTUNO01NNfldI9ARGIrzCDYBSwpWV8cbJvIbcAfhljPhFqCSewLBfUuFpH4CTMIHgFWmtlyM8sA1wAbSw8ws5Ulq28Hfh9iPRNqqc8xXHAO9AxG8fEiIpEK7R6Bu+fN7AbgLiAJfNvdt5rZjUCbu28EbjCztwBDwCHgA2HVczyjM5V19tNUl42iBBGRyIR5sxh3vxO4c8y2z5UsfyLMz5+skZnK9nX1Aw3RFiMiMsMiv1lcDjTMhIjEmYIAaKo92jQkIhI3CgIgk0owvyajKwIRiSUFQaC5Xn0JRCSeFAQBTVkpInGlIAgUJ7FX05CIxI+CINDSkGN/9wD54ULUpYiIzCgFQaClPos77O9W72IRiRcFQUAzlYlIXCkIAkc7lSkIRCReFASB0fGGNC+BiMSMgiAwvzZLwjR3sYjEj4IgkEwYTXWaslJE4kdBUKKlXn0JRCR+FAQlmutyuiIQkdhREJQYmbJSRCROFAQlWupzHOwZZCA/HHUpIiIz5oRBYGZJM/vyTBQTtZGfkHboqkBEYuSEQeDuw8CbZqCWyDVrpjIRiaHJzln8mJltBH4M9IxsdPc7QqkqIiPDTKgvgYjEyWSDIAccANaXbHNgdgVB0DS0R0EgIjEyqSBw9z8Ju5ByMLc6QzppahoSkViZ1K+GzGyxmf3MzPYFj5+a2eKwi5tpiYTRXKcpK0UkXib789HvABuBhcHj58G2WadZU1aKSMxMNgia3P077p4PHv8ENIVYV2Q0ZaWIxM1kg+CAmV0b9ClImtm1FG8ezzot9Rp4TkTiZbJB8KfAe4A9wG7gXcCsvIHcXJ+jqz9P72A+6lJERGbECX81ZGZJ4Gp3v3IG6olcS0mnsuWNk/11rYhI5Zpsz+INM1BLWTirpRaAx146FHElIiIzY7JNQw+Y2T+a2SVmdsHII9TKInL2wgYaa7Ns2tYRdSkiIjNism0f5wXPN5Zsc47taTwrJBLGulVN3L11D/nhAqmkBmgVkdltMvcIEsDX3f32GainLKxf3cxPtuzk0ZcOs3b5vKjLEREJ1WTuERSAP5uBWsrGm1Y2kkoY9zy7L+pSRERCN9l2j1+b2X8xsyVmNm/kEWplEarPpWldNpdNCgIRiYHJBsEfAx8H7gO2BI+2sIoqB+tXN7Ntbxe7DvdFXYqISKgmFQTuvnycxxlhFxel9aubAXRVICKz3nGDwMz+rGT53WP2/bewiioHK5pqWTKvSkEgIrPeia4IrilZ/vMx+y4/0Zub2eVmts3M2s3ss+Ps/5SZPW1mT5rZb8zs9EnUPCPMjPWrmnng+f30D2kyexGZvU4UBDbB8njrx+4sDk1xE3AFsAbYYGZrxhz2GNDq7q8DfgJ88YQVz6B1q5vpHyrw4PZZOb6eiAhw4iDwCZbHWx9rLdDu7tvdfRC4DbjqmDdw3+TuvcHqQ0BZTXZz8RnzyaUTah4SkVntREFwrpl1mlkX8LpgeWT9nBO8dhHwcsn6zmDbRD4E/HK8HWZ2vZm1mVlbR8fMDf2QSyd544pG7nl2H+4nyj0Rkcp03CBw96S717t7nbunguWR9fR0FRHMb9AKfGmCOm5291Z3b21qmtn5cC5b3czOQ30839E9o58rIjJTwhxIZxewpGR9cbDtGGb2FuAvgSvdveymBrss+BmpehmLyGwVZhA8Aqw0s+VmlqH4C6SNpQeY2fnANymGQFn+pV00p4rVp9UpCERk1gotCNw9D9wA3AU8A9zu7lvN7EYzG5nk5ktALfBjM3vczDZO8HaRWreqmbYdh+jsH4q6FBGRaRfqFFzufidw55htnytZfkuYnz9d1q9u5hu/fZ77f7+ft52zIOpyRESmlQbbn4QLls6hoSqt5iERmZUUBJOQSia49Kwm7t22j0JBPyMVkdlFQTBJ61c3sb97kKd2HYm6FBGRaaUgmKRLVzZhpp+RisjsoyCYpPm1Wc5bMod7tykIRGR2URCchPWrmnli5xE6usqu35uIyJQpCE7CSC9jXRWIyGyiIDgJr11YT0t9lk0KAhGZRRQEJ8HMuGxVM//23H6GhgtRlyMiMi0UBCdp3apmugbyPLLjYNSliIhMCwXBSXrTykbSSePebTM3L4KISJgUBCepNpviouXz1Z9ARGYNBcEUXLa6mfZ93bx8sPfEB4uIlDkFwRSs12Q1IjKLKAimYHljDcvmVysIRGRWUBBM0WWrm3lw+wH6BoejLkVE5JQoCKZo/epmBvMF7n56T9SliIicEgXBFF18xnxWtdTxpbu20T+kqwIRqVwKgilKJRP89TvWsPNQH7fctz3qckREpkxBcArecGYjV5x9Gv/r3ud55XBf1OWIiEyJguAU/cXbXkPBnf/+y2ejLkVEZEoUBKdoybxqPvLmFfz8iVd4+AWNPyQilUdBMA0+9uYVLGzI8fmNWxnW5PYiUmEUBNOgKpPkz9/2Gp7e3cltj7wUdTkiIidFQTBN/uB1C1i7fB5fvmsbR3qHoi5HRGTSFATTxMz4/Dtey5G+Ib766+eiLkdEZNIUBNNozcJ6NqxdyvceepFte7qiLkdEZFIUBNPs029dRU0myY3/shV33TgWkfKnIJhm82oyfPqtq3ig/QB3bd0bdTkiIiekIAjBf7xoKata6vjbXzytcYhEpOwpCEKgcYhEpJIoCEJSOg7R7iMah0hEypeCIESj4xDdqXGIRKR8KQhCtGReNR+59Aw2PvEKD7Tvj7ocEZFxKQhC9rF1Z3JGYw0f/d4Wntx5OOpyREReRUEQsqpMku9/+CIaqtO879aHefqVzqhLEhE5hoJgBiycU8UPr3s91Zkk1966md/vVa9jESkfoQaBmV1uZtvMrN3MPjvO/kvN7FEzy5vZu8KsJWpL5lXzg+teTzJhvPdbm9ne0R11SSIiQIhBYGZJ4CbgCmANsMHM1ow57CXgg8APwqqjnCxvrOEHH76IQsF57y2beelAb9QliYiEekWwFmh39+3uPgjcBlxVeoC773D3J4FCiHWUlZUtdXz/wxfRnx9mwy0PsUtzHYtIxMIMgkXAyyXrO4NtJ83MrjezNjNr6+jomJbiovSaBfV8708vorNviPfe8hB7O/ujLklEYqwibha7+83u3ururU1NTVGXMy3OWdzAdz+0lv1dA7z3lofo6BqIuiQRiakwg2AXsKRkfXGwTQIXLJ3Ld/5kLa8c7ufab23mYM9g1CWJSAyFGQSPACvNbLmZZYBrgI0hfl5FWrt8Ht/6QCsvHOjhfbdu1jSXIjLjQgsCd88DNwB3Ac8At7v7VjO70cyuBDCzf2dmO4F3A980s61h1VPO3nhmI99834U8t7eLq266n0dfOhR1SSISI1Zps2i1trZ6W1tb1GWEYvP2A3zq9ifYfaSPj192Jv9p/UoyqYq4jSMiZc7Mtrh763j79FemjFx0xnx++clLuPqCxXztnnau/voD6oUsIqFTEJSZ+lyaL7/7XL5x7YW8crift3/tfm69/wUKhcq6chORyqEgKFOXn30ad33yUi5d2cjf/MvTXHvrZnU+E5FQKAjKWFNdllve38rfX30OT7x8mMu/eh93PLqTSruvIyLlTUFQ5syMa9Yu5ZefuJRVp9Xxqduf4GPff1R9DkRk2igIKsTS+dX86CMX85nLV/ObZ/fy5i9u4it3b+NwrwJBRE6NgqCCJBPGx9at4Bf/+RLeeGYj/3BPO2/6wia+fNc2DukKQUSmSP0IKtgzuzv5x3vaufN3u6lOJ3n/G5Zx3SVnMK8mE3VpIlJmjtePQEEwCzy3t4t/+M3v+cVTu6lKJ3nfxadz3SVn0Fibjbo0ESkTCoKYaN/XxdfuaefnT7xCNpXk2tcv5bpLz6C5Lhd1aSISMQVBzDzf0c1N97Tzz4/vImHGZaub+aMLFrN+dbOGrBCJKQVBTL2wv4cfPvwSP3tsFx1dA8ypTnPluQu5+oLFnLu4ATOLukQRmSEKgpjLDxe4v30/P310F3dv3cNAvsCKphquvmAx7zx/EQvnVEVdooiETEEgozr7h7jzyd3c8eguHt5xEDN4w4r5vPP8xVy2qon5usEsMispCGRcLx7o4Y5Hd3HHYzt5+WAfZnDu4jmsW9XEZauaOWdRA4mEmo9EZgMFgRyXu/PkziPcu62De5/bx+MvH8Yd5tdkeHMQCpeubKKhOh11qSIyRQoCOSkHewa577kONm3bx2+f6+Bw7xAJK86xvG5VE2uXz+d1ixvIpZNRlyoik6QgkCkbLjhP7DzMvc/uY9O2Dp7adQSAdNI4e1EDrafP5cLT59G6bK46sImUMQWBTJuDPYNsefEQbS8eZMuOQzy58wiDwwUAljfWcOHpc2k9fS6ty+ZyRmOt7jGIlAkFgYRmID/M73Yd4ZEdh2jbcYgtLx7kUO8QANWZJK9ZUM+aBfWsWVjPaxfWc1ZLnZqURCKgIJAZ4+5s39/DlhcP8fQrncXH7k66B/JAcQTVFU01vHZhw2hAnNVSR2NtRh3cREJ0vCBIzXQxMruZGSuaalnRVDu6rVBwXj7UOxoKW1/p5MHnD/Czx3aNHlOfS7GiuZYzm2pZ0VwbvEcNS+dVk0pqWAyRMOmKQCJzoHuAp3d30r6vm+c7uoPnHjq6BkaPSSeNZfNrWNFUy/IgGEYeCxpyCgmRSdIVgZSl+bVZLlnZxCUrm47ZfqRviO0dxVAYCYnn9nbxm2f3MjR89B8uyYSxaE4VS+dVs6QkIBbPrWLBnByNNVndrBaZBAWBlJ2GqjTnL53L+UvnHrN9uODs6eznpQO9vHywl5dKHndv3cOBMbO0ZZIJTmvIsaAhx6I5xXBY0FB1dLm+ivqqlO5NSOwpCKRijFwBLJpTxcUr5r9qf/dAnpcO9LLrcB+7j/QVnw/3s/tIH5tfOMiezn6GC8c2hWZTCZrrszTX5WiuyxYf9bljn+uyzK3O6OpCZi0FgcwatdkUaxYWf4k0nuGC09E1MBoUe470s69rgH2d/eztHOC5vV3c376frv78q16bMJhXk6WxNsP82gyNtVnm12SD5czo8ryaDHNrMtRldaUhlUNBILGRTBinNeQ4rSEHzJ3wuL7BYfZ1FUNib2c/HV0DHOge5EDPAPu7B9nfPcBjLx3mQPcAPYPD475HKmHMqc4wtzrN3Jri87yaDHOqM8yrztBQnaahKs2cqnTJcoZcOqEAkRmnIBAZoyqT5PT5NZw+v+aEx/YNDnOgpxgU+7sHONQ7xOHeQQ72DHKod4hDPYMc6h3khf09bHnxMId7B8kXJv6lXiaZOCYk6qvS1OVS1OeC53HW63Mp6nJparMpqjNJBYmcNAWByCmoyiRZnKlm8dzqSR3v7nQN5DnSO8SRvqOPw8H64b5BOkvW93X183xHns6+ITr786+6xzFWwqAmm6Ium6I2l6I2m6I2ly6uZ1PUZFPUZpPUZFNUjyxnivuqS/dliqGS1s9zY0FBIDKDzIz6XJr6XJolJ/lad6dvaJjOvjxd/UN09hfDobNviO6BPN39eboH8nQFzyPrR/qG2HWod3TbRM1Z48kkE1Rnk1Snk1RliiFRlQ6eM8Xt1ZkkVZni9upMklzm6PFVmeTo9qp0klzwqMokyaUS6gdSJhQEIhXCzIJ/qaeC+xxTUygUA6VnoBgKPQPFwOgdzNM9MExvsN43OEzP4DB9g/ngeZiewTy9wT2U3sFhegeG6Rsq7hsZfPBkpJN2NByCRzadIJcKnkfCIzWynCCbKj7n0kmyqeJ6Nl36fPSYbGrkmOJyJpUgk0qQ1C/AjqEgEImZRMKoCZqJplN+uFAMhSAYegePXe4fOvroGxqmf6gwun8gX3we2d4/NExnf56OrgEG8oWS1xbozw9zqgMipBJWDId0kkwyQTadIJNMjAbFyHJ2zHpxOQiUpJFJJUiX7Esni69JJ4uvSacSpJM2+vp0snSfHbueNJIJi+Qej4JARKZFKpmgLpmgLhfuTHbuzuBwgYF8gYGhAgP54dHl/vzwMdv6h4YZzBePHXkeyI/dduz6yHt39ec5EKwPluwrfQ7DSCikRkNiZNn45FvO4h3nLpz2z1QQiEhFMbOgyScJU28hO2XuTr7gDOYLDI0JiKFhH10fKnkM5v3Y9eC4fMn60HCBoXyh+N7B8tBwgaGCMyek6WIVBCIiU2BmpJM2K35ZVflnICIip0RBICISc6EGgZldbmbbzKzdzD47zv6smf0o2L/ZzJaFWY+IiLxaaEFgZkngJuAKYA2wwczWjDnsQ8Ahdz8T+CrwhbDqERGR8YV5RbAWaHf37e4+CNwGXDXmmKuA7wbLPwH+vWmgFBGRGRVmECwCXi5Z3xlsG/cYd88DR4BXDTRvZtebWZuZtXV0dIRUrohIPFXEzWJ3v9ndW929tamp6cQvEBGRSQszCHbBMeNqLQ62jXuMmaWABuBAiDWJiMgYYXYoewRYaWbLKf7BvwZ475hjNgIfAB4E3gXc4378UUS2bNmy38xenGJNjcD+Kb623Ohcys9sOQ/QuZSrUzmX0yfaEVoQuHvezG4A7gKSwLfdfauZ3Qi0uftG4Fbge2bWDhykGBYnet8ptw2ZWZu7t0719eVE51J+Zst5gM6lXIV1LqEOMeHudwJ3jtn2uZLlfuDdYdYgIiLHVxE3i0VEJDxxC4Kboy5gGulcys9sOQ/QuZSrUM7FTnBvVkREZrm4XRGIiMgYCgIRkZiLTRCcaCTUSmJmO8zsKTN73Mzaoq7nZJjZt81sn5n9rmTbPDP7lZn9PnieG2WNkzHBeXzezHYF38vjZva2KGucLDNbYmabzOxpM9tqZp8ItlfU93Kc86i478XMcmb2sJk9EZzLfw22Lw9Gam4PRm7OTMvnxeEeQTAS6nPAf6A45tEjwAZ3fzrSwqbIzHYAre5ecZ1kzOxSoBv43+5+drDti8BBd//7IKTnuvtnoqzzRCY4j88D3e7+5ShrO1lmtgBY4O6PmlkdsAX4Q+CDVND3cpzzeA8V9r0Eg2/WuHu3maWB+4FPAJ8C7nD328zsG8AT7v71U/28uFwRTGYkVJkB7n4fxc6DpUpHof0uxf95y9oE51GR3H23uz8aLHcBz1AcELKivpfjnEfF8aLuYDUdPBxYT3GkZpjG7yQuQTCZkVAriQN3m9kWM7s+6mKmQYu77w6W9wAtURZzim4wsyeDpqOybkoZTzA51PnAZir4exlzHlCB34uZJc3scWAf8CvgeeBwMFIzTOPfsbgEwWzzJne/gOKkPx8PmilmhWCsqUptr/w6sAI4D9gN/I9oyzk5ZlYL/BT4pLt3lu6rpO9lnPOoyO/F3Yfd/TyKA3auBVaH9VlxCYLJjIRaMdx9V/C8D/gZxf9IKtneoH13pJ13X8T1TIm77w3+5y0At1BB30vQDv1T4P+4+x3B5or7XsY7j0r+XgDc/TCwCbgYmBOM1AzT+HcsLkEwOhJqcJf9Goojn1YcM6sJboRhZjXAW4HfHf9VZW9kFFqC5/8bYS1TNvJHM/BOKuR7CW5M3go84+5fKdlVUd/LROdRid+LmTWZ2ZxguYriD12eoRgI7woOm7bvJBa/GgIIfjL2Pzk6EurfRVzSlJjZGRSvAqA4aOAPKulczOyHwDqKw+nuBf4a+GfgdmAp8CLwHncv6xuxE5zHOorNDw7sAD5S0sZetszsTcC/AU8BhWDzX1BsX6+Y7+U457GBCvtezOx1FG8GJyn+g/12d78x+P//NmAe8BhwrbsPnPLnxSUIRERkfHFpGhIRkQkoCEREYk5BICIScwoCEZGYUxCIiMScgkBkAmb2l8HIj08Go1ZeZGafNLPqqGsTmU76+ajIOMzsYuArwDp3HzCzRiAD/D8qdORXkYnoikBkfAuA/SOddYI//O8CFgKbzGwTgJm91cweNLNHzezHwTg3I3NGfDGYN+JhMzsz2P5uM/tdMM78fdGcmsixdEUgMo7gD/r9QDXwa+BH7v7b0rkggquEO4Ar3L3HzD4DZIMeoDuAW9z978zs/RR75f6BmT0FXO7uu8xsTjCOjEikdEUgMo5gLPgLgeuBDuBHZvbBMYe9HlgDPBAMF/wB4PSS/T8seb44WH4A+Cczu47i8AEikUud+BCReHL3YeBe4N7gX/IfGHOIAb9y9w0TvcXYZXf/qJldBLwd2GJmF7r7gemtXOTk6IpAZBxmtsrMVpZsOo/iwGtdQF2w7SHgjSXt/zVmdlbJa/645PnB4JgV7r7Z3T9H8UqjdHh0kUjoikBkfLXA14KhgPNAO8Vmog3Av5rZK+5+WdBc9EMzywav+yuK82MDzDWzJ4GB4HUAXwoCxoDfAE/MyNmIHIduFouEoPSmctS1iJyImoZERGJOVwQiIjGnKwIRkZhTEIiIxJyCQEQk5hQEIiIxpyAQEYm5/w+1CakYHrRb7wAAAABJRU5ErkJggg==\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + " f = lambda x : math.sin(x)\n", + " errs = np.array([])\n", + " steps = np.array([1/x for x in range(1, iterations)])\n", + " \n", + " for i in range(steps.size):\n", + "\n", + " points = np.arange(0, 2*math.pi, steps[i])\n", + " apx = l2(f, points)\n", + " \n", + " err = 0\n", + " for j in range(points.size): err += apx[j]-f(points[j])\n", + " errs = np.append(errs, np.abs(err))\n", + "\n", + " plt.plot(errs)\n", + " plt.xlabel(\"Steps\")\n", + " plt.ylabel(\"Error\")\n", + " plt.show()\n", + " " + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Discussion**\n", + "\n", + "This lab showed that approximation methods can yield tight solutions. With many iteration points the approximation error quickly approaches zero when compared to the actual function value." + ], + "metadata": { + "collapsed": false + } + } + ], + "metadata": { + "kernelspec": { + "name": "python3", + "language": "python", + "display_name": "Python 3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + }, + "pycharm": { + "stem_cell": { + "cell_type": "raw", + "source": [], + "metadata": { + "collapsed": false + } + } + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/Lab-5/lab5_integration.ipynb b/Lab-5/lab5_integration.ipynb new file mode 100644 index 0000000..d33fb14 --- /dev/null +++ b/Lab-5/lab5_integration.ipynb @@ -0,0 +1,414 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "[](https://colab.research.google.com/github/johanhoffman/DD2363-VT20/blob/tree/viktorme/Lab-5/lab5_integration.ipynb)\n", + "# Lab 5: Quadrature #\n", + "**Viktor Meyer - DD2363 Methods in Scientific Computing**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Abstract**\n", + "This lab is an exercise in approximation methods of integrals. Multiple methods were implemented with different characteristics. Two of the methods provided very close approximations in a single evaluation while the third method had an iterative approach with the goal of eventually converging." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Environment**\n", + "To have access to the necessary modules you have to run this cell." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 201, + "outputs": [], + "source": [ + "# Load necessary modules.\n", + "\n", + "import numpy as np\n", + "import scipy.integrate as integrate\n", + "import math\n", + "import random\n", + "\n", + "from matplotlib import pyplot as plt" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Introduction**\n", + "The methods to be implemented were covered in the during lectures and definitions are available in the lecture notes. The aim of this lab is to approach the problem of integrating functions using approximation. The mandatory part includes 2-point Gauss quadrature over a unit interval, 3-point edge midpoint quadrature over a reference triangle, Monte Carlo quadrature over a unit interval. There is also an extra assignment Monte Carlo quadrature over a reference triangle.\n", + "> [1] LN-DD2363-part6.pdf, p. 161" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Methods**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "### 2-point Gauss quadrature over a unit interval ###\n", + "Code is commented with references to the theoretical background.\n", + "> [2] LN-DD2363-part6.pdf, pp. 205-206\n", + "\n", + "> [3] http://mathforcollege.com/nm/mws/gen/07int/mws_gen_int_txt_gaussquadrature.pdf" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 202, + "outputs": [], + "source": [ + "def randpoly():\n", + " a = random.uniform(1, 10)\n", + " b = random.uniform(1, 10)\n", + " c = random.uniform(1, 10)\n", + " d = random.uniform(1, 10)\n", + " \n", + " f = lambda x : a*math.pow(x,3) + b*math.pow(x,2) + c*math.pow(x,1) + d\n", + " F = lambda x : a*math.pow(x,4)/4 + b*math.pow(x,3)/3 + c*math.pow(x,2)/2 + d*math.pow(x,1)/1\n", + " \n", + " return [f, F(1) - F(0)]\n", + "\n", + "def gaussquadunit(f):\n", + " \n", + " # c1 = (b-a)/2. [3]\n", + " c1 = (1-0)/2\n", + " \n", + " # c2 = ((b-a)/2). [3]\n", + " c2 = (1-0)/2\n", + " \n", + " # x1 = ((b-a)/2)*(-1/sqrt(3))+((b-a)/2). [3]\n", + " x1 = (1-0)/2*(-1/math.sqrt(3))+(1-0)/2 \n", + " \n", + " # x2 = ((b-a)/2)*(-1/sqrt(3))+((b-a)/2). [3]\n", + " x2 = (1-0)/2*(1/math.sqrt(3))+(1-0)/2\n", + " \n", + " return f(x1)*c1 + f(x2)*c2" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### 3-point edge midpoint quadrature over a reference triangle ###\n", + "Code is commented with references to the theoretical background.\n", + "> [4] LN-DD2363-part6.pdf, pp. 206\n", + "\n", + "> [5] http://www.cs.rpi.edu/~flaherje/pdf/fea6.pdf" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 203, + "outputs": [], + "source": [ + "def randtet():\n", + " a = random.uniform(1, 10)\n", + " b = random.uniform(1, 10)\n", + " c = random.uniform(1, 10)\n", + " d = random.uniform(1, 10)\n", + " e = random.uniform(1, 10)\n", + " f = random.uniform(1, 10)\n", + " \n", + " fn = lambda x,y : a*math.pow(x,2) + b*math.pow(y,2) + c*x*y + d*x + e*y + f\n", + " Fn = (1/24)*(2*a+2*b+c+4*d+12*f+4*e)\n", + " \n", + " return [fn, Fn]\n", + "\n", + "def gaussquadtriangle(f):\n", + " \n", + " # An 3-point rule which is exact for quadratic integrands is obtained by choosing\n", + " # the quadrature points as the midpoints of the three edges, with weights\n", + " # w0=w1=w2=1/6,motivated by the symmetry of the quadrature points. [4]\n", + " w0 = 1/6\n", + " w1 = 1/6\n", + " w2 = 1/6\n", + " \n", + " # As expected the optimal evaluation point is the centroid of the triangle. [5]\n", + " x0 = [0, 0.5]\n", + " x1 = [0.5, 0.5]\n", + " x2 = [0.5, 0]\n", + " \n", + " return f(x0[0],x0[1])*w0 + f(x1[0],x1[1])*w1+f(x2[0],x2[1])*w2" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Monte Carlo quadrature over a unit interval ###\n", + "Monte Carlo integration samples a function at random. That is, the algorithm evaluates n samples at random positions over its interval. With a large enough samples size (law of large numbers) the result will converge to the exact solution.\n", + "> [6] LN-DD2363-part6.pdf, pp. 224-225" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 204, + "outputs": [], + "source": [ + "def montecarlounit(f, n):\n", + " \n", + " result = 0\n", + " \n", + " for i in range(n):\n", + " result += f(random.uniform(0,1))\n", + " \n", + " result /= n\n", + " \n", + " return result" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Results**" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "### 2-point Gauss quadrature over a unit interval ###\n", + "We check that the approximation is close enough to the exact solution." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 205, + "outputs": [ + { + "name": "stdout", + "text": [ + "Test Completed Successfully!\n" + ], + "output_type": "stream" + } + ], + "source": [ + "def test_gaussquadunit():\n", + " \n", + " for i in range(1, 1000, 1):\n", + " poly = randpoly()\n", + " apx = gaussquadunit(poly[0])\n", + " assert np.isclose(apx, poly[1])\n", + " \n", + " print(\"Test Completed Successfully!\")\n", + " \n", + "test_gaussquadunit()" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### 3-point edge midpoint quadrature over a reference triangle ###\n", + "We check that the approximation is close enough to the exact solution." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 206, + "outputs": [ + { + "name": "stdout", + "text": [ + "Test Completed Successfully!\n" + ], + "output_type": "stream" + } + ], + "source": [ + "def test_gaussquadtriangle():\n", + " \n", + " for i in range(1, 1000, 1):\n", + " tet = randtet()\n", + " apx = gaussquadtriangle(tet[0])\n", + " assert np.isclose(apx, tet[1])\n", + " \n", + " print(\"Test Completed Successfully!\")\n", + " \n", + "test_gaussquadtriangle()" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "### Monte Carlo quadrature over a unit interval ###\n", + "We check that the error converges towards zero as the number of samples increases. By inspection it is apparent that the convergence rate resembles 1/sqrt(samples)." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": 207, + "outputs": [ + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5wU9fkH8M9zhQ4CiogUDywoYgFPBFssEYkaNWqMGo0aE1KIKSYxGFOMGqNGY4/KD0vsWIkBRQHp/QDp7ehHuzuuc3V3n98fM7M7O9vm9navDJ/363Wv292Znflue+Y7z7eMqCqIiMh7Mlq6AERElB4M8EREHsUAT0TkUQzwREQexQBPRORRWS1dALujjjpKc3JyWroYRERtxvLly4tVtVe0Za0qwOfk5CAvL6+li0FE1GaIyM5Yy5iiISLyKAZ4IiKPYoAnIvIoBngiIo9igCci8igGeCIij2KAJyLyKE8E+GdnbsGczUUtXQwiolbFEwH+37PzsSC/uKWLQUTUqngiwAsEvHAJEVE4bwR4ARjfiYjCeSPAt3QBiIhaIU8EeABgBZ6IKJwnAryIMEVDROSQ1umCRWQHgEoAfgA+Vc1Ny34AKOvwRERhmmM++ItVNb19GNnISkQUwRspmpYuABFRK5TuAK8AvhSR5SIyNtoKIjJWRPJEJK+oKLnRqL6A4quNhU0pJxGR56Q7wJ+vqsMBfAvAOBG50LmCqk5Q1VxVze3VK+plBROqrvdjV0k18naUNLG4RETekdYAr6p7zP+FAD4BMCKd+yurbkjn5omI2pS0BXgR6SwiXa3bAEYDWJuu/Rn7SefWiYjalnT2oukN4BMxom4WgHdUdVoa90dERDZpC/Cqug3AGenafjSswRMRhXiimyQREUXyVIAX9ognIgryVIAnIqIQbwV4VuCJiIK8FeCJiCjIUwGeFXgiohBPBXgiIgrxVIAXdoQnIgryVIAnIqIQTwV41t+JiEI8FeCJiCjEUwGeKXgiohBPBXgiIgrxVIDnXDRERCGeCvBERBTiqQDPHDwRUYinAjwREYV4KsCzAk9EFOKpAE9ERCHeCvCswhMRBXkrwBMRUZCnAjz7wRMRhXgqwBMRUYinAjz7wRMRhXgqwBMRUYinAjwr8EREIZ4K8EREFOKpAM9rshIRhaQ9wItIpoisFJEp6d4XERGFNEcN/lcANjTDftiLhojIJq0BXkT6AbgSwMR07oeIiCKluwb/NIB7AQTSvB8A7EVDRGSXtgAvIlcBKFTV5QnWGysieSKSV1RUlK7iEBEddtJZgz8PwNUisgPAewAuEZG3nCup6gRVzVXV3F69ejVph8zBExGFpC3Aq+p9qtpPVXMA3ATgK1W9NV37IyKicJ7qB6/a0iUgImo9sppjJ6o6G8Ds5tgXEREZvFWDb+kCEBG1Ip4K8EREFOKpAM8cPBFRiKcCPBERhXgqwCur8EREQZ4K8EREFOKpAM/6OxFRiKcCPBERhXgqwDMFT0QU4qkAT0REIZ4K8MosPBFRkKcCPBERhXgrwLMCT0QU5K0AT0REQZ4K8KzAExGFeCrAExFRiKcCPPvBExGFeCrAExFRiKcCPPvBExGFeCrAExFRiKcCPHPwREQhngrwREQU4qkAzwo8EVGIpwI8ERGFeCrA85qsREQhngrwREQU4qkAz/o7EVGIpwI8ERGFeCvAswpPRBSUtgAvIh1EZKmIrBKRdSLyt3Tti4iIImWlcdt1AC5R1SoRyQYwX0Q+V9XF6doh56IhIgpJW4BXo89ilXk32/xjBCYiaiZpzcGLSKaIfA2gEMB0VV0SZZ2xIpInInlFRUVN2h+7wRMRhaQ1wKuqX1XPBNAPwAgRGRplnQmqmququb169UpncYiIDivN0otGVcsAzAIwJr37SefWiYjalnT2ouklIt3N2x0BXAZgY7r2R0RE4dLZi6YPgP+ISCaMA8n7qjoljftjCy4RkU3CGrzZUPpEYzesqqtVdZiqnq6qQ1X1weSKmNifrjwlXZsmImqzEgZ4VfUDOL8ZypK0kYOOBMDZJImI7NymaFaKyKcAPgBwyHpQVT9OS6mIiKjJ3DaydgBwEMAlAL5t/l2VrkIla+yby1u6CERErYarGryq3pnugjSFSOh2bYMfHbIzW64wRESthKsavIj0E5FPRKTQ/PtIRPqlu3DJqG3wt3QRiIhaBbcpmtcAfArgWPPvf+ZjrYIgVIVnOysRkcFtgO+lqq+pqs/8ex1Aq5xXgPGdiMjgNsAfFJFbzT7xmSJyK4xG11bBnoMPsApPRATAfYD/IYAbAewHsA/ADQBaZcMr4zsRkSFhLxpzqoHrVPXqZihPUuw1eF70g4jI4HYk683NUJaUsNfgaxv8+OcXG9mzhogOS25Hsi4QkecBTEL4SNYVaSlVI8XqRfPagh14YdZWdGqXhXEXn9ACJSMiajluA/yZ5n/7hGEKY2Rrq2JP0dT5jJp7HWvwRHQYcpODzwDwoqq+3wzlSUpYDp4peCIiAO5y8AEA9zZDWVIiajdJ+xGAiOgw4bab5AwR+Z2I9BeRntZfWkvWCPbwzRo8EZHBbQ7+e+b/cbbHFMCg1BaHiIhSxe1skgPTXZCm4EhWIqJIcVM0InKv7fZ3HcseSVehmoLxnYjIkCgHf5Pt9n2OZWNSXJYmsPWDb8FSEBG1JokCvMS4He1+q8DrshIRGRIFeI1xO9r9FhOeg2+5chARtSaJGlnPEJEKGLX1juZtmPc7pLVkSYsS4VmrJ6LDUNwAr6pt4uKmsfrBS+vMIhERNQu3A53aDKZoiIgMngjwIvZeNIzwRESARwK8HdPtREQGTwR4zkVDRBTJEwHejlMVEBEZ0hbgzZknZ4nIehFZJyK/St++Qrf/+MkalFXXp2tXRERtRjpr8D4Av1XVIQBGAhgnIkPSuD8AwOqCcjwzc0u6d0NE1OqlLcCr6j7rmq2qWglgA4C+6dhXrP7u7FFDRIezZsnBi0gOgGEAlkRZNlZE8kQkr6ioKCX7y+AVnIiI0h/gRaQLgI8A/FpVK5zLVXWCquaqam6vXr2S3Ef4/QzzPkeyEtHhLK0BXkSyYQT3t1X143Tuy441eCKi9PaiEQCvANigqv9K135i7Ls5d0dE1CqlswZ/HoDbAFwiIl+bf1ekcX9BGYzvRESuL7rdaKo6H810UZDIHDwjPBGR50ayAqzBExEBHgnwzpw7c/BERB4J8E5M0RAReSTAO8N5pideFRFR03gyFDpTNJywgIgOR54I8LF60TBTQ0SHM08EeCf2oiEi8kiAd845w0ZWIiKPBHgiIorkiQDvrLD7zcv2ldc0AACe+yofu0uqm7tYREQtyhMB3sm6Lusr87cHH1tdUN5SxSEiahGeCPDOjHsgENkxkg2vRHS48USAdwoo4HcEeba7EtHhxhsB3pmDDyjOfPDL+CsREXmcNwK8g6qistYX9pi9Br98Zwmq630gIvIyTwR4Zz94qxdN+DqG4qo6XP/iItwzaVUzlIyIqOV4IsA7vTBra8Rj1uCnmno/AGDNnsb1qqlt8De9YEREzcgTAd5NA2pTGlkX5Bfj5D9Pw+JtB5PfCBFRM/NEgG+MZAL9wq3FAIBl20tSXBoiovTxRIB3E7ObMj+NldLPYGd6ImpDPBHgXXHEZo3SEBuL1aWefemJqC3xRIBvzDVYrXUbcxEQ62DAWSqJqC3xRIB3oynB2ZrbhuGdiNoSTwR4N4H39leXQlWDtfFGZGhCOXjW4ImoDfFEgHfr4KH6YLDeX1EbfHxvWQ0mr9wT83nMwRNRW+SJAO828AZUg+kWAME54r/70iL8etLXaPAHYj4PYA2eiNoWTwR4twKBUG0cQHDg0r7ymrjPs9I6TYnvj03biJzxU+GLcRAhIkq1rJYuQCo456KJJWDLwQOh3HpAQ8ujsR5NtgZ/97sr8b9VewEAvoAiKzOpzRARNYqnavDdOsQ/XvkDiijXAgmK1fAaStEkVy4ruMfbBxFRqqUtwIvIqyJSKCJr07UPi7rs1a4aPsBpxa5S5IyfGrY8muBBIQU5eLdlJSJqqnTW4F8HMCaN24+QaMCT0cgauv/est0Ry6PRJtbgw/fR9G0QEbmRtgCvqnMBNMvsXG7THn5HLxqn2AHe+O821x+3DIzwRNRMWjwHLyJjRSRPRPKKioqauK34y42BTrGXx4q9Tc3BO8vg1OAPsHcNEaVciwd4VZ2gqrmqmturV6/ktuFyPX8gdi3dLEvUx63An4p+8NEOIkP/+gVG/uOrJm+biMjOE90kLYnCbyBBDT7WsmCKJgU1+GgpmjpfAHVVdU3fOBGRTYvX4FPB7dS/S7eX4NvPz4+5PHEja3iE//V7K5EzfirqfO4v59eYaYqJiJoind0k3wWwCMBgESkQkbvSta9QL8b4VeyJ87fFXZ4oB2/ffFFlHSZ/bfRvt67zCgBbDlTitAe+wH8W7kBZdb3rfVimrz+A3Ien8xqwRNRk6exFc7Oq9lHVbFXtp6qvpG9fxv9EGZSCUndTEjhFC8rlNaHgbR1Y/rdqL+76Tx4qa33466frMPbN5RHP8yeowf9p8hoUV9Xj4KHIg0NL2Li/Auv2Nu4C5UTUOngiRWP1bunSIQs9O7eLuV6i7Eis2rVGWW7v9LJpfyUAY0qCXeYEZgCwo/hQ5D4SVOFLDzUAALJayeUBxzw9D1c+GzutRUStlycC/JFd2uP+K07Bmz88BzPu+UbS27GPMv3T5DWYunofgFCKxl7Dt+frb3x5UdTtFVbWRZwV2O8GAhrRPbLevM/+8kTUVJ4I8ADw4wsHYcCRndCzczv07d4xqW28NHsr/vDhagDAW4t3Ydw7KwCEAruG1eDdBeDCyvDeMfYDwx2vL8MJ938e9XnpCPAf5O3G92IcjIjIezwT4FPhP4t2YlLe7ojHA2Yl217DjzV3vJOz3deeg5+7OfbArnj99ZP1+w9XY8n2pg0u3l1SjT1lNfjxG3lxL5LSFgUCyjMn8hRP9YNPFyso23/79b7wAP/h8gJX23LbTdJtoFm09SDOOq4H2mXFP1Ynyv27dcHjs4K3p68/gGuH9U16W9X1PnRq13q+gne/txJTV+/DjkevbOmiEKUEa/AuWMF2dUGoN0m9owb/uw9WxX2uxW2cdVODX7e3HDf/32L84/MNCdedl1/sbsfNJL+wEkP+8gU+XuHuwNgcrDYXr5mx/gA2H6hs6WJQC2CAT2DWpsJgkH536a7g484afCw+f3igdlsz97lY72CV0ZVyy4GqhOu6LW8sh+p8TXq+00az59GMDQdSul2K9KM38jD6qbktXQxqAZ4M8KkcLXrna8uiBmW3ATOyBp/aFI3Tyl2l2HWwOuJxe6/LZNI15z/WuLly8naUYMzTc2MO2ErFzJxEFJ8nA3wsbueScR4gogZ4l42szpq422NPIMrm63x+VNY2oLCyNubzvvPvhbjwn7MiHre/9kSDraIprW7AU9M3u17/b/9bj437KxOmBlrDzA31vgAHc5EnHVYB3vW88Y6gHC0g1jUkV4N3WzOPts8LHpuF0x74EiP+PjNsO24OXPYac7I9dJ6ZucXVerUN/mBK5+rnF+Dud1dGlscsTrIB/oO83bj/kzXJPdnhwSnrOJjLlF9YFbd3F7Utngzw9phx2ZDejX6+M7jag6k/oDhQUeu60dLnqIonk6KZtbEQdT5/WJ96XyAQfJ2xNllUWYd5W4qwML84bB6HaGcH0by9ZKe7FR2ufn4+ttlG8dqvSWuximN1PT1U52vU/Du//3A13l6yK/GKLuTtKE3Jdlo7VcV7S3fFbU/55r/m4AevLm3GUlE6eTLA2yWTby6sCB+cZA+2z8zYjEufnBM1aEXT1F40y3aU4M7Xl+FfX4anR6KdCby5ODwg3zRhEW57ZSlumbgEczaFamVuUzT3f5Lc5XQ3u2j0dZ51nPrXL/CNKKklN8qrG/Db91ehorYhYtmczUVYv7ci7vOdabRUdSlNpMEfQGmcOYcWbi1GzvipWLvHffpoYX5xzAC+aOtBjP94DR783/pGl5XaJs8H+GTyzde+sCB8G7Yf/LNf5aOqET1KnD1c4jUAL98Zqkla+9xbZkyQtqUwfDtD/vIFdpvz3ljB8u9Tw3+4W4tCtWh78HOTJmpsQ3WDP4Cvd5c16jnGfkK3D1Q0fk784Q9NxxkPfomPVhTgzUWRZxy3v7oUVzw7L+423KTk0uHeD1dj2EPTY77XM9YXAgAWbzvoant7y2pwy8Ql+P2H0bvsVprf29YykR2lnycDvL1ymExlzPkDaMroxt86+sfH29aEuVuDtwMBxeSVe/Cr974GAHy1sTBi/Smrw88i4l1x6uMVoVGnbmqobhuRLU98uQnXvrAAt05c4mp9K6bFKsr24kP41MVZUonts3IGymijjRv8AbyftzvsPUi2p1NTWf3uK2pT0wW1ut7YjtUF1cl6fzJd/Op/GmUmVIquqs6HvB3NcvnpRvNmgLcFusG9uzR5e6kcvh5vU/Y+837VsH730TjL5bbjoZsaal0j+82v22OkQea7bpuIX4ZLnpyNX0ZpnI3HeT2AAxWRvY0mzN2Gez9cjf+uCh3wnO/j1sJDeGzaxrADxsR521Ie9KyDaLw0TWNE+1jDJ8gz/gskePYXy7R1+1NSJq+prvfhjUU7MHHeNqzcVYqDVXW47t8LcMNLi1BeE5kibGmtZ5x4CmWYh617LjsJt4/Kwf/N296k7aXylD1e6sPnqFUm2muDYxBVvx6dsMnFiEU3NVS3vYSSUVZdb2tQjXUVrabvx5lKU1XsKzdSXpW2WrOzIfzO15fiQEUdbhkxAP17dkJRZR0enpp4tHBlbQMCChzRMTvhuit2hdJxpdX1yEFnvLl4JyYt24VPx52PsiSCRfCsyz5jadht4860dfsxbd1+TB53Hs7s3z3m9mob/Kiq8+GoLu0bXRav+vvUDcHG/SM7t0NJdX3wu+p2fqrm5M0Ab9bkzjquBzIzmz6gJpU1+JW7y3Byn25R5623B5qAasJUijMwndq3GzYdqMRpfY+I+7xYvWgOVtVhT1kNhvTp5uoyhPYKs5uumhc/MRtD+nTD1DWhKQH8AY170AsEFBlJzo3f4Attd9mOEox7e0WwJ1J2ZgYqahuQnZEB5+/SukLXwUP1wQBvUdWwM4ULHv8KWRkZmPW7i3DWQzNQ7w9Encumut4Hf0DRtYMR/JdsC53Sl1UbwfzPk41G7V9P+hqfrtqLzu0yI7bjDyge+HQdfnj+QAw8qnPYsmhnXfbv7v7y8DOa/MKquAH+jteWYvG2Evz5qiH44Xk5Ca+YdjiwRo8Dkanc5mqcbwxPpmisAJ8hkpJRrakM8P/8YhPueC16N7SwFE0gcU3bvv5na/YFLzCSqLz2M5KF+cX4ZKUxH8xFT8zG1c8vwB8/WeNqpK4qcMkTsxOuZ9lefCgsuAPArE1F+N0Hq2M+p8Ftn06EArOqImf8VDw6LVTr/u5Li8K6mWZlCE5/4EuMfnoO/BFdWY3/176wwJiz37Z83pZiXPXcvOABcHdJDbab77tVg/7X9M0RPXrOffQrnPbAl8H79tN558HUans4ZL4e42LxRqE27KvAm4t3Rk1fRTvrsn+HnGchztfttNg8CD00ZT3mbomeevvFOytww4sL426npdX5/GFtNclYkF+MpQlmYh3xyMwm7SMdPBngrYpGZobETHNccOJRrreX6ilkVxeUR+0S6Oxv70+wWyulM29LMX7+9gqs2FVmPh7/h2uvadwycQl+M8loCLbSFh8uL3Cdg98W5apVjfWRY8Ixe28i51w+8RwyGxmt1NWC/Ni9T7LMM7vdJTUR+7C/fyXV9WGpsN99sApr91TEzWE/O3MLHvt8Y/B+IKDBWrrFXvGo8wUiatd2j3y2EVeYA7HsAfulOVvDBnvVmgcK+6uJ990tOeQ+DeS8MI1lyup9yNsZexyBastPwfzjN5Zj+EPTXa27p6wGC/OLUVxVFzZi/PsTl+DGlxchUeI0UYXSH1Cs2FXabLV9TwZ4qwafmQG0jzGN7s8uOh4z7rnQ1fbSkVvbGWW+GHvOeNHW4oRfllg/ukQNmImCd0DDT0VToTFf6MenhYLjPz7fgJJD9VjgovHW6v9d6yK9lG3rSuJ8v2ptNeH95bVh73ND8Ipb8bdv9YxZXVCGQX/8LGK5PejV+wIY+Y/4tb8N+yrCnpchwKOfbwwb7GXV4PeV16DcPKDEaz96zPY+J5KZIE1mlau2wR8sKwA8Om0jjv/jZzG/q40RCCjmb0n8u3CKNjJ31sZC5Iyfins/XIWrn5+PP09ei61FVRjz1FzcMnEJch+egRF/j/xMEu3aqgxsOVAZtZyvzt+O6/69EF+ub55J9jwZ4K2vYoYI2mdlYsH4SyLW6dI+Cycc3dXV9hrbZTBZ9lzvfxbtTDiy09nIaklU6/3xG3kAEFYL/dyROvnle+57sKzdU54wP1tU5b6P+9k5PYO331q8C7dOXILvT1wCnz+ANxftCI4NcHo/rwDjP1qNbUWNO6uIV8P8eMUefG/C4uD9UjNwPjV9czC1BUQewBrMg2isswh74H1x9tao60QT7AkT5f22Uj21DQGMfnpO1HJFc9Vz85Azfiqe/HJTzHWyMjJQWduA52ZuiVrhsVIgv/tgFb71zLzgAeb1BTvMsgVQXtOAoX/9Agu3hh+sd5dU4ydv5gVTbLG8tWQnbn1lCb5YFxkcJy3bhWdmbEFhlJ5TFvt78eoCo+PF+3kFWF1QjjcX78SlT84JjhWIJdG7WefzI29HCS57ai7eXLwz4gxmW7ExnmXnwaaf+brh6UZWq9bRs1N4g+ZLt56F0/sZjUt9u3fEnhgBw+I8vU7k6K7tIy7V99dvD8HfEowgdDba7Itz2g7EDkyJaktWzvg3k74OPvazt1eErdOYnOVVz81Hx+zIBkG7RN3yLHU+f0QAWW/WCNfurcCf/7sO7+fFnkP+vWW746YMLPY2hngH8EnLondVtXqiWNY4RpuGavrRtx3+o3f3Y/ebNVggMhVwz6SvUVAa+h5bg8binc1lZwrqfH6sNbu4PvdVfsx1MwTBNoQz+nfHhSf1Clt+qM6HXl3bB/P2c7YUhbUT1PsC2LC/AlV1Pjw7cwvOPd5IkRaUVgcvIjNncyHGDO0TfM7eshoc061DsJHdGkg3Ye5W+AOKK08PrfuHj4xU1esLt2PlX0ZHfQ31/gA6ZMT/niaSqAZfXe8PVjDWFJTjx28sx4Z9FcFKZr0vfABjunkywFuVGyvQZzt60owZekzwdlYKetk4ZUcZSTJsQI9Gb6cywQCYWLn2veW1WLunHF07ZMXdhptUhls1Cc42bnjJ3bVgB/9pWsxlawqMH7gzmDq5Sam5nd++R+d2OFSf+Md4jWP0c4MZWGMF2GTm53957lY8NcOYssL+ufoDio+jXD5RE/TEOqJjO9z48uKYy+3sKYVo6ZqqOh8enrIexeaZ2msLwrsm1/sDwQnvrCCpqnjSNgWH9bt5ZsYWZGcJHp+2Cc/cdCauOTP8qmErdpVhxTsr8MhnHfHmXSMwxXahltI4lbFLn5yDh68ditmbIgcNuhc/wp/zyEyMOdWIL5kZErzeQXlNA47omI3SaqPiNHNjIf52TROK4ZInA7xzRGe8/GFmCrp+ffOUozFjQ+hL4zygAMDJx7hLB516bDesSzB3iqU4Tp78qufmY0DPTnEDfCpee3N6fFrsFIJdvBG9FrejR+214saYu7kIL83ZGtmA6w9gX3kt3lsWee3fRHYWh86C7L1wznzwy2ir41C9P24OvriqLhiQE3l94Y7g7ZW7StG5fRba2Soyn67ai4nzQ0G93BFo632B4DUJrCK9OGcrPrEdmNplZUBVgwcxwDjzc3ZNtewpq8H3Jy6Jeaa7aX8lKm29mfaU1eDO15cBAEYO6hn1OU4n3f85bsjtF7zvJv1vndnZu/euKSjH0L7dgmfGBaU1aPAHolYGU8mTAd45FW28/HCyfazjifahZbnczzVnHus6wCdSXtOAH4w6Dm9EmaMFgPsJ8luJRPlRi/NVndS7S8QEaBXNMOrw8WkbMfbC48MeO+H+z5Penv1s094NM9ZBfMXOUgzq1TnqsqZ44svNeMIx+d10R6OhM+10weOzMKBnJwDA0h0luPr5+WGXwASM9gVnBwBrX4vvuzRqN9B4aczLn459FasalwP56v0BvJPkrKX25936yhL079kxrPJRcqgevbt1SGrbbnmykdV6ExOPBQ2/0lHywjeSFSXAJ+qFYOmYwotQl9c0oHunyAFVgDHwpirK7Itt1VW2fKzzoB1tdsuX525Le5kCavSmSBV7xSFWA7vdD15d2qiJ8ZrCzX522dphnMEdAOoa/DE7FqwqKHM9nfQ9tralWFYlMTEekLiRNZ7dJTXYebAafY4wgvqU1fsiznRSzaMB3vjvpmdezpGpr+G0i5KicTsK8IITjsIPRh2XsrJ0jzFsfs6mIlTU+nDziP4p21cyRg06MqnndWkffiDsYGvkbcy88sm4+oxjE65zRj9jNPHMKJPEJWv9vsaf2d352rKU7T+eVJwRzd5chCUxBhPtKa1x3Wb08co9Sde6E0lFn/4Texvp2oemrMcfP1mTcABVU3gywFvB1M2H8eSNZ+B7udGDnLOnQDyjbRcWiZVX+9lFx6NTu0y886Nzoi6/6vQ+yDmqMx68ZmjY45ef2hs32vKAlkljRyYsV4/O2ejaIfKs4Gdvr0BRZR26dczGn648JeF2mmLshYNiLrvitGNiLovnlnMGhN23Hz6defNfXnJCUvuI5fR+saeC6Nu9I64981j8YFROSvcJIKlAkKgnVqo0dnK6aN5Zsgs/iTGhW1FVXdj4hET+mKKrfTkl6kzghr09buqafbjx5UVpO9PyZIAPNeYkDvBdO2Tj4pOPjrqsIcqXNtoIWBGj66XFCvBHdQlPj/xhzMlY/+AYnHtC9FG03WLUtl++LRf3jjk54vFzotR+nZmgY7p1RKcoc5pYOmVn4Y5zc2IuT4WffeP4mMu6RDn4uNHZkcrKyszAOQOjN5yNPjX2QSRWjvrK0yQwrNoAAA75SURBVEIpH2f7Sa+uxuRb0RrTP/rZuXj6pmH4VoID168uPTHu8tbi7JzYvb9+dP7AZivHW4t3ho1wbill1ZEdGxqb5h11fOTv9okvNqVkWhWntAZ4ERkjIptEJF9ExqdzX3b//O4ZuOr0PjgjykRKb941IuKxWBNrWd3tvmGryT/yndMi1hMYed8PfjoKt5wzIJhvv3Wku1TL4N5dce+YwfiDLYh/8NNRYeu0izEi1+n8E3uFnXn0OaIDOsXJ61fX+6K2GQDAN08JnZV89LNREcvjHRheu+Ps4G17+sR5NnH+CdHPkqIFT7vO7cMPWvW+ACb9JLKMk8edh+6doh84f3jeQBzVOfpMifYpAZztJz3Mdo0hfbrhv+POCx5Y/nHdaTjGzK92apcVdtC/4azwM7BYlYrWwpoMz5lavPO8nOBta2qIxnKm19xI1GW4KcbEqQA4RWvPaWzWZuTAyAD/+sIdaZnMLW0BXkQyAbwA4FsAhgC4WUSGpGt/dsf36oLnbxkeNVViDbCws2ZfvHWkcdr/+8sH47U7zsZPzZrnszcPAwCMGNgT/Xt2wpa/fyvs+Vajydk5PfHId07Dd810yvXDI9MqFnuevXP7TPz8ohPCppk9o1/4wSk7I/y1fPLzcwFEnlGcdHQXvPHDEfjk5+fi6jOORb8ekTX48044Er+97CQA8QfCTLw9N3g62SHKQCZn0LK7+OSjMf03F+LeMYPRsV0mrh/eD9cP7xcMlo9dfxq++u03grVhu4XjL8G8eyNHH9tZ2xkxsCeyMiT42Tl17ZCFru2jB/jsTEG3jkawefS608Ley92lRoPg2Tk9IrpdVpsjLo/q0h5n9O8enCHSOa3u8OOMz/D1O8/G49efHrbszP7dMebUY/C93P64deQAbHxoDJ656Ux88NNRuP+K9KbMAEQc9EYP6Y3fXz4YI3J64vEbTsd4s7LhrJ3+5arQT/jorsb3ftzFsc/QLIN7h9IS4y42UmZu2jLcuHiwUUn40fkDseaB0XjtzrMj2rEuG9I7aipyyt3nu6qIRTtz//lFx+PUY7vhxxe4P5P5/jkD0LFdJj746Sh0yE5/AiWd3SRHAMhX1W0AICLvAbgGQIteEDJab5ZBvboEp3h9+NrwGrr1+NL7L0U384dsnbJ3bZ+FR647LeLC3tec2Tc4OGPqL8/H3rLIPOgD3z4V1w/vh2teWBA112/V2K3TYHuN9svfXIiTzB/MhNtyUVxVFxwN+NvRgwEYA6uswVUnHN0l2PXyH9edhhtz+6PBH8DBQ/X4xcWh/PSQPt0wtG83vJ9XEAwA1lw+ndpl4ZZzBuCzNftwTLcOuH54PwztewT+ecPpuPjko3H3OyuxuqAMuTk9cd8VRnA4sXfXYIPSkzeeAQC48aVFWLqjBKf17Y5BvYyLsSz546XwBxTnPvoVAODY7h0BABseHIOPVhTguCM7YcLcbSiqrMPxR3dBeXVD8P25Mbc/3rfV3O++5ATk7SjFpgOVGHhUZ/Tt3hHtMjNw1nE9gqf4/Xt2xO6SGmRmCH5z2UnYfKAK3xzSG9X1fswzR4r+bvRg5Ob0RLvMDPz4jTzMsc1nYr031qn2XecPxBEdsyOCwNFdO4RNHfy93P74cEUBXrjFqDC8dNtZYetb35nT+x2BKav34tJTeuNf08O7Iz54zan4y3/XBe8/f8sw7C+vxSl9uiFDBBW1Dbj81GOQM34qAOC5m4fh7ndXYuyFgzDB7Dk0IqcnHv7OUIx+KtSNcMIPcgGEgm/JoXqcvKArfn/5YFz/4iJ0bpeJ9386CiKCV27PRXZmBkYOOhJDju2G0UN6463Fu8L65n9r6DEYPqAH5m4pwn3fOgUrd5fi/k/WYvK487DfnJP/2mHH4s7zctCzczss21GKIzpm459fbHR1TV/LRYN74envDcOsTYUYM/QYdMjOxMWDj8bIgUeib/eOKDlUj5fnbsNzNw9Dh+zMsBk1x144CEPNyt32f1yBYQ9NR1l1AzIzBCce3QWVtT7sKavB2r9djk37K4PfjZ6d22HqL89HnyM64t4xJ0NVMeDIzujeMRt3O2b5vG3kcRh74SDsr6jF36duwL2XG7+Ns3N6YviAHli49SCyMgQ/uiB2O1WTqDkVaar/ANwAYKLt/m0Ano+y3lgAeQDyBgwYoOkyfd1+XZBflLLtvbtkp24trGzydrYXVWlNvc/VumWH6rXB54+6bPamQp28siDqsqraBv2/uVv1o+W7Y267pt6n9ea2q+t8wTIVlFbrC7O2aCAQcFXGRPaX1+jklQVRt7diZ4m+u2Snq+3U+/w6Y/3+RpVrb1m1fr5mr9b7/Pri7HzdXXIobLnPH9DJKwvU7w/fZlVtg64pKNMDFTWab37mq3aXpuw9iefrXaX6zpKdeqC8RhfmF6uq6to9ZfrR8t36xqIdMZ8XCAR0f3lNzGWWg1V1un5vedwyTJy3TXcWH4q7zqG6Bp04b5tuLazUySsLtLou/Dvt9weC2wgEAnG35/cHdMuBCv18zT59evpmXbS1WIsra/WFWVu0vKZe95XV6Lo95Tpx3ra4ZbK2daiuIXh/9e4y3by/QgtKqyN+d1sOVOi/Z+Xrwao6bfD5Nb+wMvh9DAQC+u6SnfrVhgO6t6w66r4CgYAWlFbrK/O2aVFlbfCxWA6U1+h/Fm5v8vcIQJ7GiMOiabr+pIjcAGCMqv7IvH8bgHNU9RexnpObm6t5eXlpKQ8RkReJyHJVzY22LJ1JoD0A7P0P+5mPERFRM0hngF8G4EQRGSgi7QDcBODTNO6PiIhs0tbIqqo+EfkFgC8AZAJ4VVXXJXgaERGlSFonG1PVzwBEXs6GiIjSzpMjWYmIiAGeiMizGOCJiDyKAZ6IyKPSNtApGSJSBCDG5YfiOgpAccK1vIWv+fDA13x4aMprPk5Vo87a16oCfLJEJC/WSC6v4ms+PPA1Hx7S9ZqZoiEi8igGeCIij/JKgJ/Q0gVoAXzNhwe+5sNDWl6zJ3LwREQUySs1eCIicmCAJyLyqDYf4Fvqwt7pJiL9RWSWiKwXkXUi8ivz8Z4iMl1Etpj/e5iPi4g8a74Pq0VkeMu+guSJSKaIrBSRKeb9gSKyxHxtk8zppyEi7c37+ebynJYsd7JEpLuIfCgiG0Vkg4iM8vrnLCK/Mb/Xa0XkXRHp4LXPWUReFZFCEVlre6zRn6uI3G6uv0VEbm9MGdp0gG/JC3s3Ax+A36rqEAAjAYwzX9t4ADNV9UQAM837gPEenGj+jQXwYvMXOWV+BWCD7f5jAJ5S1RMAlAK4y3z8LgCl5uNPmeu1Rc8AmKaqJwM4A8Zr9+znLCJ9AfwSQK6qDoUxnfhN8N7n/DqAMY7HGvW5ikhPAH8FcA6M61z/1ToouBLrWn5t4Q/AKABf2O7fB+C+li5Xml7rfwFcBmATgD7mY30AbDJvvwzgZtv6wfXa0h+MK3/NBHAJgCkABMYIvyznZw7jWgOjzNtZ5nrS0q+hka/3CADbneX28ucMoC+A3QB6mp/bFACXe/FzBpADYG2ynyuAmwG8bHs8bL1Ef226Bo/QF8VSYD7mKeYp6TAASwD0VtV95qL9AHqbt73yXjwN4F4AAfP+kQDKVNVn3re/ruBrNpeXm+u3JQMBFAF4zUxLTRSRzvDw56yqewA8AWAXgH0wPrfl8PbnbGns59qkz7utB3jPE5EuAD4C8GtVrbAvU+OQ7pl+riJyFYBCVV3e0mVpRlkAhgN4UVWHATiE0Gk7AE9+zj0AXAPj4HYsgM6ITGV4XnN8rm09wHv6wt4ikg0juL+tqh+bDx8QkT7m8j4ACs3HvfBenAfgahHZAeA9GGmaZwB0FxHr6mP21xV8zebyIwAcbM4Cp0ABgAJVXWLe/xBGwPfy5/xNANtVtUhVGwB8DOOz9/LnbGns59qkz7utB3jPXthbRATAKwA2qOq/bIs+BWC1pN8OIzdvPf4DszV+JIBy26lgm6Cq96lqP1XNgfFZfqWq3wcwC8AN5mrO12y9FzeY67epmq6q7gewW0QGmw9dCmA9PPw5w0jNjBSRTub33HrNnv2cbRr7uX4BYLSI9DDPfEabj7nT0o0QKWjEuALAZgBbAdzf0uVJ4es6H8bp22oAX5t/V8DIPc4EsAXADAA9zfUFRo+irQDWwOih0OKvowmv/yIAU8zbgwAsBZAP4AMA7c3HO5j3883lg1q63Em+1jMB5Jmf9WQAPbz+OQP4G4CNANYCeBNAe699zgDehdHG0ADjTO2uZD5XAD80X3s+gDsbUwZOVUBE5FFtPUVDREQxMMATEXkUAzwRkUcxwBMReRQDPBGRRzHAkyeJyP3mbIWrReRrETknjfuaLSKH1UWiqW3ISrwKUdsiIqMAXAVguKrWichRANq1cLGImh1r8ORFfQAUq2odAKhqsaruFZG/iMgycw7yCeYoSqsG/pSI5JnzsZ8tIh+b828/bK6TY87X/ra5zoci0sm5YxEZLSKLRGSFiHxgziUEEXlUjLn9V4vIE834XtBhjAGevOhLAP1FZLOI/FtEvmE+/ryqnq3GHOQdYdTyLfWqmgvgJRjDx8cBGArgDhGxZi4cDODfqnoKgAoAP7fv1DxT+BOAb6rqcBijU+8xn/8dAKeq6ukAHk7DayaKwABPnqOqVQDOgnHhhCIAk0TkDgAXm1cEWgNjIrNTbU+z5jBaA2Cdqu4zzwC2ITTZ025VXWDefgvGdBJ2I2FceGaBiHwNY66R42BMb1sL4BURuQ5AdcpeLFEczMGTJ6mqH8BsALPNgP4TAKfDmONjt4g8AGOOE0ud+T9gu23dt34nznk9nPcFwHRVvdlZHhEZAWNSrRsA/ALGAYYorViDJ88RkcEicqLtoTNhXCEHAIrNvPgNkc9MaIDZgAsAtwCY71i+GMB5InKCWY7OInKSub8jVPUzAL+BcVk+orRjDZ68qAuA50SkO4xr2+bDSNeUwZi9cD+MqaYbaxOMa+O+CmN627DroapqkZkKeldE2psP/wlAJYD/ikgHGLX8e5LYN1GjcTZJIhfMyyZOMRtoidoEpmiIiDyKNXgiIo9iDZ6IyKMY4ImIPIoBnojIoxjgiYg8igGeiMij/h/5cCXy23BvwwAAAABJRU5ErkJggg==\n" + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def test_montecarlo():\n", + " \n", + " samples = np.array([])\n", + " errs = np.array([])\n", + " poly = randpoly()\n", + " \n", + " for i in range(1, 1000, 1):\n", + " apx = montecarlounit(poly[0], i)\n", + " samples = np.append(samples, i)\n", + " err = apx-poly[1]\n", + " errs = np.append(errs, np.abs(err))\n", + "\n", + " plt.plot(samples, errs)\n", + " plt.xlabel(\"Samples\")\n", + " plt.ylabel(\"Error\")\n", + " \n", + " plt.show()\n", + " \n", + "test_montecarlo()" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n", + "is_executing": false + } + } + }, + { + "cell_type": "markdown", + "source": [ + "# **Discussion**\n", + "\n", + "This lab showed that approximation methods can yield tight solutions. Monte Carlo iteration was very interesting since it is very simple while still being able to yield accurate results. That being said, there is of course the downside of randomness and the need for large sample sizes." + ], + "metadata": { + "collapsed": false + } + } + ], + "metadata": { + "kernelspec": { + "name": "python3", + "language": "python", + "display_name": "Python 3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + }, + "pycharm": { + "stem_cell": { + "cell_type": "raw", + "source": [], + "metadata": { + "collapsed": false + } + } + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file