cvcyl25.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def coeff():\n",
    "    for i in range(1,n):\n",
    "        for j in range(1,m):\n",
    "            aP[i,j]=((r[i]+r[i+1])/2)*k/(r[i]*h)+((r[i]+r[i-1])/2)*k/(r[i]*h)+(h/k)+(h/k)\n",
    "            aE[i,j]=((r[i]+r[i+1])/2)*k/(r[i]*h)\n",
    "            aW[i,j]=((r[i]+r[i-1])/2)*k/(r[i]*h)                                                                                                                                                                                                                                                                                                                                                                                                                              \n",
    "            aN[i,j]=h/k\n",
    "            aS[i,j]=h/k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss():\n",
    "    for k in range(1,p):\n",
    "        for i in range(1,n):\n",
    "            for j in range(1,m):\n",
    "                U[i,j]=(aW[i,j]*U[i-1,j]+aE[i,j]*U[i+1,j]+aS[i,j]*U[i,j-1]+aN[i,j]*U[i,j+1])/aP[i,j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h= 0.1\n",
      "k= 0.2\n",
      "r= [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.  1.1 1.2 1.3 1.4 1.5 1.6 1.7\n",
      " 1.8 1.9 2. ]\n",
      "z= [0.  0.2 0.4 0.6 0.8 1.  1.2 1.4 1.6 1.8 2.  2.2 2.4 2.6 2.8 3.  3.2 3.4\n",
      " 3.6 3.8 4. ]\n",
      "max1= 0.0009710719687276992\n",
      "iteration= 153\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAHwCAYAAAC7cCafAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYZHV97/H3l56NGdaBEWQd9qgEkIzIoqhhEXHBJV5B3K4maK7E7bqgJop5jBdRMbkm0aAiqIhGkUgUEdQYrqjIADMsAooCMqzDomwDzDDf+8c5zdT0VFVX93T16fr1+/U89ZyqU6d+53vmdM+nf+f86pzITCRJUjk2aLoASZI0sQx3SZIKY7hLklQYw12SpMIY7pIkFcZwlySpMIa7NGAi4sSI+Gof2r0pIg6d6HYlTT7DXZpE/QrmUkTE7Ij4YkTcHBEPRMQVEfGCEcscEhHXRcTDEfFfEbFjj22/JyKurtu9MSLeM+L9hXV7D9ft+4eOBpbhLk0hUZnOv5czgFuA5wCbAn8H/HtELASIiC2Bb9fz5wOLgW/02HYArwM2B44Ajo+Io1vePwu4AtgC+CDwrYhYsH6bIzVjOv8nInUVEdtHxLcjYnlE3BMR/1zP3yAi/rbuXd4VEV+OiE3r9xZGREbE6yPi9xFxd0R8sH7vCOADwKsi4sGIWFrP/0lE/ENEXAw8DOwcEdtExLkRcW9E3BARf9VjzddGxItaXs+oa9i3fv2SiLgmIv5Qr/cpHdo5PSI+2vL6uRGxrOX1TXVP+MqIeKjubW8VEd+ve8Y/jIjNW5bfPyJ+Vq93aUQ8t916M/OhzDwxM2/KzNWZ+V3gRuDP6kVeDlyTmd/MzEeAE4G9I+JPImJ+RCyLiBfX69yo/rd7Xd32yZl5eWauyszrge8AB9XL7g7sC3w4M1dk5tnAVcArevl3l6Yaw11qIyKGgO8CNwMLgW2Br9dvv6F+PA/YGdgI+OcRTTwL2AM4BPhQRDwlM88HPgZ8IzM3ysy9W5Z/LXAcsHG9zrOAZcA2wF8AH4uIQ3oo/SzgmJbXzwfuzszL6wA7C3gHsAA4D/jPiJjVQ7vtvAI4DNgdeDHwfao/Xrak+r/lbQARsS3wPeCjVL3tdwNn99Irjoit6vavqWc9DVg6/H5mPgT8FnhaZt4LvBH4fEQ8Cfg0sCQzv9ym3QCePaLd32XmAy2LLa3nSwPHcJfa248qWN9T9yYfycyf1u8dC5ySmb/LzAeB9wNHR8SMls9/pO4BLqUKib3p7vTMvCYzVwFbU/1x8L56vUuAL1D9ATCarwEviYi59etX1/MAXgV8LzMvzMyVwCeBDYEDe2i3nc9k5p2ZeSvw/4BLMvOKzHwUOAd4er3ca4DzMvO8ujd+IdXh9CO7NR4RM4EzgTMy87p69kbAH0cs+keqP4rIzAuAbwI/Al4IvLlD8ydS/f/3pV7alQaN4S61tz1wcx22I21D1bsedjPVueKtWubd0fL8Yarw6OaWEe3fO6IXeTPV0YOuMvMG4FrgxXXAv4Q14b5W3Zm5ul7vqO12cGfL8xVtXg9v847AK+tD8n+IiD9Q/fHy5E4N1+MOvgI8Bhzf8taDwCYjFt8EaP23OhXYE/hSZt7Tpu3jqc69v7D+Q6TXdqWBYbhL7d0C7DCiNz7sNqrAGrYDsIq1w62TTrdhbJ1/GzA/Ilp7jTsAt/bQPqw5NH8U8Ks68IfbfaLu+tD09h3afQiY2/J66x7X3c4twFcyc7OWx7zMPKndwnVdX6T6Y+kV9VGGYdfQchQkIuYBu9Tzh0+n/BvwZeCvI2LXEW2/ETgBOCQzl7W8dQ3VWIfWf/O9WXPYXhoohrvU3i+B24GTImJeRMyJiIPq984C3hkRO0XERqw5j96ulz/SncDCbiPiM/MW4GfA/6nXuxfwJqpD1L34OnA48Nes6bUD/Dvwwqi+SjYT+N/Ao/W6RloCHFkPUtua6jz9eH2V6kjC8yNiqN6m50bEdh2W/yzwFODFmblixHvnAHtGxCsiYg7wIeDKlsP2H6inb6Q67fDlOvCJiGOp9tVhmfm71kYz89f1Nn+4ru9lwF7A2eux3VJjDHepjcx8nGqQ2K7A76kGt72qfvs0qkPGF1GN5H4E+Jsem/5mPb0nIi7vstwxVAP5bqMKtA/X56p7qf124OdU59K/0TL/eqrz358B7qbavhdn5mNtmvkK1ViBm4AL6P3rZu3quYXqKMIHgOVUPfn30Ob/n6i+s/5mYB/gjvpbBQ/WwUxmLqcayPcPwH3AM4Gj68/+GfAu4HX1/vs41RGRE+rmP0r1NbdLW9r9XMvqjwYW1e2eBPxFvT5p4ERmp6OEkiRpENlzlySpMIa7JEmFMdwlSSqM4S5JUmEMd0mSCtPuAh1TzpZbbpkLFy5sugxJkibFZZdddndmjvuuhAMR7gsXLmTx4sVNlyFJ0qSIiJtHX6ozD8tLklQYw12SpMIY7pIkFcZwlySpMIa7JEmFMdwlSSqM4S5JUmEMd0mSCmO4S5JUGMNdkqTCGO6SJBXGcJckqTCGuyRJhTHcJUkqjOEuSVJhDHdJkgpjuEuSVBjDXZKkwhjukiQVxnCXJKkwhrskSYUx3CVJKozhLklSYfoW7hFxWkTcFRFXt3nv3RGREbFlv9YvSdJ01c+e++nAESNnRsT2wGHA7/u4bkmSpq2+hXtmXgTc2+atTwPvBbJf65YkaTqb1HPuEfES4NbMXDqZ65UkaTqZMVkrioi5wAeBw3tc/jjgOIAddtihj5VJklSWyey57wLsBCyNiJuA7YDLI2Lrdgtn5qmZuSgzFy1YsGB45iSVKknS4Jq0cM/MqzLzSZm5MDMXAsuAfTPzjp4aeNvbYIst+lmiJElF6OdX4c4Cfg7sERHLIuJN69kgPP74hNQmSVLJ+nbOPTOPGeX9hWNqcGgIVq9en5IkSZoWBucKdUND9twlSeqB4S5JUmEMd0mSCmO4S5JUmMEKd3BQnSRJoxi8cLf3LklSV4a7JEmFMdwlSSqM4S5JUmEMd0mSCmO4S5JUGMNdkqTCGO6SJBXGcJckqTCGuyRJhTHcJUkqjOEuSVJhDHdJkgpjuEuSVBjDXZKkwhjukiQVxnCXJKkwhrskSYUx3CVJKozhLklSYQx3SZIKY7hLklQYw12SpMIY7pIkFcZwlySpMIa7JEmFMdwlSSqM4S5JUmEMd0mSCmO4S5JUGMNdkqTCGO6SJBXGcJckqTCGuyRJhTHcJUkqjOEuSVJhDHdJkgpjuEuSVBjDXZKkwhjukiQVxnCXJKkwhrskSYUZnHCPqB6GuyRJXQ1OuEPVezfcJUnqqm/hHhGnRcRdEXF1y7xPRMR1EXFlRJwTEZuNqVHDXZKkUfWz5346cMSIeRcCe2bmXsCvgfePqUXDXZKkUfUt3DPzIuDeEfMuyMxV9ctfANuNqVHDXZKkUTV5zv2NwPfH9AnDXZKkUTUS7hHxQWAVcGaXZY6LiMURsXj58uXVTMNdkqRRTXq4R8TrgRcBx2ZmdlouM0/NzEWZuWjBggXVTMNdkqRRzZjMlUXEEcD7gOdk5sNjbsBwlyRpVP38KtxZwM+BPSJiWUS8CfhnYGPgwohYEhGfG1OjhrskSaPqW889M49pM/uL69Wo4S5J0qi8Qp0kSYUx3CVJKozhLklSYQx3SZIKY7hLklQYw12SpMIY7pIkFcZwlySpMIa7JEmFMdwlSSqM4S5JUmEMd0mSCmO4S5JUGMNdkqTCGO6SJBXGcJckqTCGuyRJhTHcJUkqjOEuSVJhDHdJkgpjuEuSVBjDXZKkwhjukiQVxnCXJKkwhrskSYUx3CVJKozhLklSYQx3SZIKY7hLklQYw12SpMIY7pIkFWbwwj2zekiSpLYGL9zB3rskSV0Y7pIkFcZwlySpMIa7JEmFMdwlSSqM4S5JUmEMd0mSCmO4S5JUGMNdkqTCGO6SJBXGcJckqTCGuyRJhTHcJUkqjOEuSVJhDHdJkgpjuEuSVBjDXZKkwvQt3CPitIi4KyKubpk3PyIujIjf1NPNx9So4S5J0qj62XM/HThixLwTgB9l5m7Aj+rXvTPcJUkaVd/CPTMvAu4dMfso4Iz6+RnAS8fUqOEuSdKoJvuc+1aZeTtAPX1SpwUj4riIWBwRi5cvX17NNNwlSRrVlB1Ql5mnZuaizFy0YMGCaqbhLknSqCY73O+MiCcD1NO7xvRpw12SpFFNdrifC7y+fv564Dtj+rThLknSqPr5VbizgJ8De0TEsoh4E3AScFhE/AY4rH7dO8NdkqRRzehXw5l5TIe3Dhl3o4a7JEmjmrID6traoC7XcJckqaPBCnd77pIkjcpwlySpMIa7JEmFMdwlSSqM4S5JUmEMd0mSCmO4S5JUmMEM99Wrm61DkqQpbDDD3Z67JEkdGe6SJBXGcJckqTCGuyRJhTHcJUkqjOEuSVJhDHdJkgozWOHu/dwlSRrVYIV7RBXwhrskSR0NVrhDdWjecJckqSPDXZKkwhjukiQVxnCXJKkwhrskSYUx3CVJKozhLklSYQx3SZIKY7hLklQYw12SpMIY7pIkFcZwlySpMIa7JEmFMdwlSSqM4S5JUmEMd0mSCmO4S5JUGMNdkqTCGO6SJBXGcJckqTCGuyRJhTHcJUkqjOEuSVJhDHdJkgpjuEuSVBjDXZKkwhjukiQVxnCXJKkwhrskSYUx3CVJKkwj4R4R74yIayLi6og4KyLm9Pxhw12SpK4mPdwjYlvgbcCizNwTGAKO7rkBw12SpK6aOiw/A9gwImYAc4Hbev6k4S5JUlc9hXtEzI2Iv4uIz9evd4uIF41nhZl5K/BJ4PfA7cAfM/OCnhsw3CVJ6qrXnvuXgEeBA+rXy4CPjmeFEbE5cBSwE7ANMC8iXtNmueMiYnFELF6+fPmaNwx3SZK66jXcd8nMk4GVAJm5AohxrvNQ4MbMXJ6ZK4FvAweOXCgzT83MRZm5aMGCBWveMNwlSeqq13B/LCI2BBIgInah6smPx++B/etD/QEcAlzb86cNd0mSuprR43IfBs4Hto+IM4GDgDeMZ4WZeUlEfAu4HFgFXAGc2nMDhrskSV31FO6ZeWFEXA7sT3U4/u2Zefd4V5qZH6b6g2HsDHdJkrrqKdwj4uD66QP19KkRQWZe1J+yujDcJUnqqtfD8u9peT4H2A+4DPjzCa9oNIa7JEld9XpY/sWtryNie+DkvlQ0mqGharp6NWwweJfGlySp38abjsuAPSeykJ4Nh7u9d0mS2ur1nPtnqL8GR/UHwT7A0n4V1VVruM+c2UgJkiRNZb2ec1/c8nwVcFZmXtyHekZnz12SpK56Ped+Rr8L6ZnhLklSV13DPSKuYs3h+HVk5l4TXtFoDHdJkroarec+fOe3t9bTr9TTY4GH+1LRaAx3SZK66hrumXkzQEQclJkHtbx1QkRcDPx9P4try3CXJKmrXr8KNy8injX8IiIOBOb1p6RRGO6SJHXV62j5NwGnRcSm9es/AG/sT0mjMNwlSeqq19HylwF7R8QmQGTmH/tbVheGuyRJXY02Wv41mfnViHjXiPkAZOYpfaytPcNdkqSuRuu5D59X37jfhfTMcJckqavRRsv/Wz39yOSU0wPDXZKkrnoaLR8RJ0fEJhExMyJ+FBF3R8Rr+l1cW8PXk1+5spHVS5I01fX6VbjDM/N+qovaLAN2Z+17vE+e2bOr6aOPNrJ6SZKmul7Dffj2a0dS3TTm3j7VMzrDXZKkrnr9nvt/RsR1wArgf0XEAuCR/pXVheEuSVJXPfXcM/ME4ABgUWaupLqu/FH9LKwjw12SpK56HVA3l+rmMZ+tZ20DLOpXUV3NmlVNH3uskdVLkjTV9XrO/UvAY8CB9etlwEf7UtFo7LlLktRVr+G+S2aeDKwEyMwVQPStqm4Md0mSuuo13B+LiA2BBIiIXYBm0tVwlySpq1FHy0d1IfnPAecD20fEmcBBwBv6W1oHw+HuOXdJktoaNdwzMyPi7cDhwP5Uh+Pfnpl397u4toYH1NlzlySprV6/5/4LYOfM/F4/i+mJh+UlSeqq13B/HvDmiLgZeIiq956ZuVffKuvEnrskSV31Gu4v6GsVY7HBBtXNYwx3SZLa6incM/PmfhcyJrNnO6BOkqQOev0q3NQya5Y9d0mSOhjMcJ8923CXJKkDw12SpMIY7pIkFWZww90BdZIktTWY4e6AOkmSOhrMcPewvCRJHRnukiQVxnCXJKkwgxvuDqiTJKmtwQx3B9RJktTRYIa7h+UlSerIcJckqTCGuyRJhRnccHdAnSRJbQ1muDugTpKkjgYz3GfPhpUrYfXqpiuRJGnKGdxwBw/NS5LURiPhHhGbRcS3IuK6iLg2Ig4YUwOGuyRJHc1oaL3/BJyfmX8REbOAuWP69HC4e95dkqR1THq4R8QmwMHAGwAy8zFgbF3wWbOqqeEuSdI6mjgsvzOwHPhSRFwREV+IiHljasGeuyRJHTUR7jOAfYHPZubTgYeAE0YuFBHHRcTiiFi8fPnytd803CVJ6qiJcF8GLMvMS+rX36IK+7Vk5qmZuSgzFy1YsGDtNx1QJ0lSR5Me7pl5B3BLROxRzzoE+NWYGrHnLklSR02Nlv8b4Mx6pPzvgP85pk87oE6SpI4aCffMXAIsGncD9twlSeposK9QZ7hLkrSOwQ53B9RJkrSOwQ53e+6SJK1jMMPdAXWSJHU0mOFuz12SpI4Md0mSCjPY4e6AOkmS1jHY4W7PXZKkdQxmuM+or71juEuStI7BDPeIqvduuEuStI7BDHeowt1z7pIkrWOww92euyRJ6zDcJUkqzOCG+6xZhrskSW0Mbrjbc5ckqa3BDncH1EmStI7BDnd77pIkrcNwlySpMIMb7g6okySprcENd3vukiS1Ndjh7oA6SZLWMdjhbs9dkqR1GO6SJBVmcMPdAXWSJLU1uOFuz12SpLYGO9wdUCdJ0joGO9ztuUuStI7BDvfHH68ekiTpCYMb7rNmVVN775IkrWVww3327GpquEuStJbBD3cH1UmStJbBD3d77pIkrcVwlySpMIMb7g6okySprcENd8+5S5LU1uCHuz13SZLWYrhLklQYw12SpMIMbrg7oE6SpLYGN9wdUCdJUluDH+723CVJWovhLklSYQY33DfaqJo+8ECzdUiSNMUMbrhvuilEwL33Nl2JJElTyuCG+9AQbLYZ3Hdf05VIkjSlDG64A8yfb89dkqQRBjvcN9/ccJckaYTBDnd77pIkrWPww91z7pIkraWxcI+IoYi4IiK+O+5G7LlLkrSOJnvubweuXa8Whs+5r149MRVJklSARsI9IrYDXgh8Yb0amj+/CnYvZCNJ0hOa6rn/I/BeYP263PPnV1MPzUuS9IRJD/eIeBFwV2ZeNspyx0XE4ohYvHz58vYLDYe7g+okSXpCEz33g4CXRMRNwNeBP4+Ir45cKDNPzcxFmblowYIF7VvafPNqas9dkqQnTHq4Z+b7M3O7zFwIHA38ODNfM67GPCwvSdI6Bv977mC4S5LUYkaTK8/MnwA/GXcDw4flPecuSdITBrvnPmcOzJ1rz12SpBaDHe7gzWMkSRph8MPdS9BKkrQWw12SpMKUEe4OqJMk6QmDH+6ec5ckaS2DH+4elpckaS1lhPuKFfDII01XIknSlFBGuIPn3SVJqg1+uHvzGEmS1jL44e715SVJWovhLklSYQx3SZIKM/jh7p3hJElay+CH+yabwNCQPXdJkmqDH+4RXqVOkqQWgx/uAFttBbfd1nQVkiRNCWWE+y67wG9/23QVkiRNCWWE+667VuGe2XQlkiQ1roxw32WX6vryt9/edCWSJDWunHAHD81LkkQp4b7rrtX0hhuarUOSpCmgjHDfYYfqu+723CVJKiTcZ86EHXe05y5JEqWEO6wZMS9J0jRXTrjvsos9d0mSKCncd90V/vAHL0MrSZr2ygl3vw4nSRJQYrh7aF6SNM2VE+4771xN7blLkqa5csJ97lzYZht77pKkaa+ccAfvDidJEqWF+267wa9/3XQVkiQ1qqxwf9rT4K674O67m65EkqTGlBfuANdc02wdkiQ1yHCXJKkwZYX7ttvCJpvA1Vc3XYkkSY0pK9wjqt67PXdJ0jRWVrjDmnDPbLoSSZIaUWa433NPNWpekqRpqLxw33PPauqheUnSNFVeuDtiXpI0zZUX7ltvDZtvbrhLkqat8sLdEfOSpGmuvHAHR8xLkqa1MsN9zz3hvvvg1lubrkSSpElXZrgvWlRNL7202TokSWpAmeG+zz4wYwb88pdNVyJJ0qQrM9znzIG99zbcJUnT0qSHe0RsHxH/FRHXRsQ1EfH2vqzomc+sDsuvXt2X5iVJmqqa6LmvAv53Zj4F2B94a0Q8dcLXst9+8MADcN11E960JElT2aSHe2benpmX188fAK4Ftp3wFe23XzX10LwkaZpp9Jx7RCwEng5c0ua94yJicUQsXr58+dgb32MP2Hhjw12SNO00Fu4RsRFwNvCOzLx/5PuZeWpmLsrMRQsWLBj7CjbYAJ7xDLhknb8bJEkqWiPhHhEzqYL9zMz8dt9W9MxnwpVXwooVfVuFJElTTROj5QP4InBtZp7S15Xttx+sWgVXXNHX1UiSNJU00XM/CHgt8OcRsaR+HNmXNR1wQDW96KK+NC9J0lTUxGj5n2ZmZOZemblP/TivLyvbaiv40z+FCy/sS/OSJE1FZV6hrtVhh8FPfwoPP9x0JZIkTYryw/3QQ+Gxx+Dii5uuRJKkSVF+uB98MMyc6aF5SdK0UX64z5sHBx4IP/xh05VIkjQpyg93qM67X3EFjOdKd5IkDZjpEe6HHlpNf/zjZuuQJGkSTI9wX7QINt8czj236UokSeq76RHuQ0PwylfCf/xHdRtYSZIKNj3CHeB1r6u+6/7t/l3KXpKkqWD6hPuBB8LOO8NXvtJ0JZIk9dX0CfcIeO1rq0F1y5Y1XY0kSX0zfcId4DWvgUw488ymK5EkqW+mV7jvuiscdBB87nPwyCNNVyNJUl9Mr3AH+MhH4Kab4FOfaroSSZL6YvqF+yGHwCteAR/7GNxyS9PVSJI04aZfuAN88pOwejW85z1NVyJJ0oSbnuG+cCG8733wjW/ABRc0XY0kSRNqeoY7wAknwO67w1veAg891HQ1kiRNmOkb7nPmwKmnwo03woknNl2NJEkTZvqGO8BzngN/9Vdwyinwy182XY0kSRNieoc7wMknw7bbwqtf7U1lJElFMNw32wy++tXq8Pzb3tZ0NZIkrTfDHeDgg+EDH4DTT4evfa3paiRJWi+G+7APfQie9Sz4y7+EJUuarkaSpHEz3IfNnAnf/CbMnw8vfSncfXfTFUmSNC6Ge6utt4ZzzoE77qguUevNZSRJA8hwH+kZz4AvfQkuuqi6//vjjzddkSRJY2K4t3PMMdVd4771LXjHO6p7wEuSNCBmNF3AlPWud8Ftt1UhP3NmNY1ouipJkkZluHfziU/AypXw6U9XvfdTTjHgJUlTnuHeTQT84z+umd5zD3zhCzBrVtOVSZLUkeE+moiq577llvB3fwe33FKdi99ii6YrkySpLQfU9SIC/vZvq8vUXnwx7LUXXHhh01VJktSW4T4Wxx4Lv/gFbLopHH44HHdc9Z14SZKmEMN9rPbdFy67DN75zur78LvuCn//9/DQQ01XJkkSYLiPz4YbViPnr7kGnv98+PCHq5D//Oer0fWSJDXIcF8fu+8OZ58NP/0p7LRTdZh+t93gX//VS9dKkhpjuE+Egw6qBtp997uwzTbw1rdWYf+pT8EDDzRdnSRpmjHcJ0oEvPCFVcj/+MfwtKfBu98N224Lxx9fHcKXJGkSGO4TLQKe9zz44Q/hkkvgZS+rLnyz557w3OfC178OK1Y0XaUkqWCGez/ttx+ccQYsWwYf/zjcfHN1U5qttoI3vAEuuABWrWq6SklSYQz3ybDllvDe98INN1Q9+le+srpv/POfD9ttVw3EO+ccuP/+piuVJBUgcgBuZ7po0aJcvHhx02VMrEcegfPOg7POgh/8oBp4N2MGPOtZ1QVyDjwQFi2CefOarlSSNMki4rLMXDTuzxvuU8DKlfCzn8H3vw/nnw9Ll1bzh4aqS93uv3/1OOCA6vv03plOkopmuJfonnuqwXg//3l1udtLLlnzlbpNN60Cf++9q8dee1WD9ebObbZmSdKEMdyng8cfh2uvrcL+8surnv1VV8GDD1bvb7BBdfGcpzwF9tijeuy+ezXdcstma5ckjdn6hru3fB0EQ0NV73zPPdfMW70abryxCvqlS+HKK+G66+B731v7Erjz51dBv+OO1WPhwjXPd9wRNtpo0jdHktRf9txLs2oV3HQTXH/9mscNN1Rfw/v979e99v2mm8LWW1ePrbZq/3zBguqPhLlzPd8vSZNgIHvuEXEE8E/AEPCFzDypiTqKNGNGNehu112rK+a1Wr26ukXtzTevedx2WzXvjjtgyZJq2ukrebNnVyG/xRbVY/h563TzzWGTTdZ9zJtXnT6QJPXdpId7RAwB/wIcBiwDLo2IczPzV5Ndy7SzwQbVte+32aYaed/JihVw551V0N9+O9x9N9x7bzXQr3X661+vef3YY6Ovf+ON1w78jTaqjgbMnVvdaa912su8OXOqPzhmzaoew89neLZJ0vTWxP+C+wE3ZObvACLi68BRgOE+VWy4YXVufuHC3pbPhIcfrkL+vvuqkf33379mOvxo9/ree6vPrlhRTYefr8/pog02WDvsW5+PnLY+nzmz+sOg9TE0tO68TvPHMm9oqKqzqYenV6SiNRHu2wK3tLxeBjyzgTo0USKqw+7z5sEOO6x/e5nw6KPtQ3/4+cMPV8s89ti603bzOr33wANrnq9cWY1ZePzxatr6GJ43cszCIGsN+pEPaD+/ifcmqt1hI/+w6fZ6LMuuz2ebWPYzn6lO36lITYR7uy7DOt20iDgOOA5gh4kIDA2OiOqQ+5w5TVfS3urVo/8RMNq81aun1iNzzQPWft3kexPV7rCRR4S6vR7LsuP5bGuN/ahptGW9r0XRmgj3ZcD2La+3A24buVBmngqcCtVo+ckpTerB8GHUAEhbAAAHkklEQVR/SZqiJv2rcBExA/g1cAhwK3Ap8OrM7HjD84h4ALh+cipsxJbA3U0X0Uclb1/J2wZu36Bz+wbXHpm58Xg/POk998xcFRHHAz+g+ircad2CvXb9+nzfb6qLiMVu32AqedvA7Rt0bt/gioj1urhLI98ZyszzgPOaWLckSaXzqiKSJBVmUML91KYL6DO3b3CVvG3g9g06t29wrde2DcS15SVJUu8GpecuSZJ6NKXDPSKOiIjrI+KGiDih6XrWV0RsHxH/FRHXRsQ1EfH2ev6JEXFrRCypH0c2Xet4RcRNEXFVvR2L63nzI+LCiPhNPd286TrHIyL2aNlHSyLi/oh4xyDvv4g4LSLuioirW+a13V9R+b/17+OVEbFvc5X3psP2fSIirqu34ZyI2KyevzAiVrTsx881V/noOmxbx5/FiHh/ve+uj4jnN1N17zps3zdatu2miFhSzx+ofQdd82Bifv8yc0o+qL4m91tgZ2AWsBR4atN1rec2PRnYt36+MdX3/Z8KnAi8u+n6JmgbbwK2HDHvZOCE+vkJwMebrnMCtnMIuAPYcZD3H3AwsC9w9Wj7CzgS+D7VVSb3By5puv5xbt/hwIz6+cdbtm9h63JT/dFh29r+LNb/zywFZgM71f+3DjW9DWPdvhHvfwr40CDuu7rmTnkwIb9/U7nn/sQNZjLzMWD4BjMDKzNvz8zL6+cPANdSXWu/dEcBZ9TPzwBe2mAtE+UQ4LeZeXPThayPzLwIuHfE7E776yjgy1n5BbBZRDx5ciodn3bbl5kXZObwtVd/QXWVzIHTYd91chTw9cx8NDNvBG6g+j92yuq2fRERwP8AzprUoiZQlzyYkN+/qRzu7W4wU0wQRsRC4OnAJfWs4+tDLacN6mHrWgIXRMRl9f0BALbKzNuh+oEGntRYdRPnaNb+j6WU/Qed91eJv5NvpOoNDdspIq6IiP+OiGc3VdR6avezWNq+ezZwZ2b+pmXewO67EXkwIb9/Uznce7rBzCCKiI2As4F3ZOb9wGeBXYB9gNupDjcNqoMyc1/gBcBbI+LgpguaaBExC3gJ8M16Vkn7r5uificj4oPAKuDMetbtwA6Z+XTgXcDXImKTpuobp04/i0XtO+AY1v7jemD3XZs86Lhom3kd9+FUDveebjAzaCJiJtWOPDMzvw2QmXdm5uOZuRr4PFP8cFk3mXlbPb0LOIdqW+4cPnxUT+9qrsIJ8QLg8sy8E8raf7VO+6uY38mIeD3wIuDYrE9o1oes76mfX0Z1Xnr35qocuy4/iyXtuxnAy4FvDM8b1H3XLg+YoN+/qRzulwK7RcROdU/paODchmtaL/V5oi8C12bmKS3zW8+bvAy4euRnB0FEzIuIjYefUw1cuppqv72+Xuz1wHeaqXDCrNVrKGX/tei0v84FXleP2t0f+OPw4cNBEhFHAO8DXpKZD7fMXxARQ/XznYHdgN81U+X4dPlZPBc4OiJmR8ROVNv2y8mub4IcClyXmcuGZwzivuuUB0zU71/TIwZHGU14JNUIwt8CH2y6ngnYnmdRHUa5ElhSP44EvgJcVc8/F3hy07WOc/t2phqRuxS4ZnifAVsAPwJ+U0/nN13remzjXOAeYNOWeQO7/6j+SLkdWEnVM3hTp/1FdVjwX+rfx6uARU3XP87tu4Hq3OXw7+Dn6mVfUf/cLgUuB17cdP3j2LaOP4vAB+t9dz3wgqbrH8/21fNPB94yYtmB2nd1zZ3yYEJ+/7xCnSRJhZnKh+UlSdI4GO6SJBXGcJckqTCGuyRJhTHcJUkqjOEuDbiI+Fk9XRgRr57gtj/Qbl2Spja/CicVIiKeS3VHsBeN4TNDmfl4l/cfzMyNJqI+SZPHnrs04CLiwfrpScCz6/tZvzMihqK6d/ml9Y1E3lwv/9z6PtJfo7oYBhHxH/XNfq4ZvuFPRJwEbFi3d2bruuqrZH0iIq6OiKsi4lUtbf8kIr4V1T3Tz6yvxEVEnBQRv6pr+eRk/htJ082MpguQNGFOoKXnXof0HzPzGRExG7g4Ii6ol90P2DOr238CvDEz742IDYFLI+LszDwhIo7PzH3arOvlVDcn2RvYsv7MRfV7TweeRnXd64uBgyLiV1SXQ/2TzMyI2GzCt17SE+y5S+U6nOpa1EuobiW5BdU1twF+2RLsAG+LiKVU9zffvmW5Tp4FnJXVTUruBP4beEZL28uyunnJEmAhcD/wCPCFiHg58HCbNiVNEMNdKlcAf5OZ+9SPnTJzuOf+0BMLVefqDwUOyMy9gSuAOT203cmjLc8fB2Zk5iqqowVnAy8Fzh/TlkgaE8NdKscDwMYtr38A/HV9W0kiYvf6bn0jbQrcl5kPR8SfAPu3vLdy+PMjXAS8qj6vvwA4mC53GavvWb1pZp4HvIPqkL6kPvGcu1SOK4FV9eH104F/ojokfnk9qG05Va95pPOBt0TElVR3DPtFy3unAldGxOWZeWzL/HOAA6juwpXAezPzjvqPg3Y2Br4TEXOoev3vHN8mSuqFX4WTJKkwHpaXJKkwhrskSYUx3CVJKozhLklSYQx3SZIKY7hLklQYw12SpMIY7pIkFeb/A0Sm/g9W07tWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbf78490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solution converge \n",
      "[[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 9.03101384e-02 1.87755980e-01 2.99935563e-01\n",
      "  4.35504586e-01 6.04785174e-01 8.20468562e-01 1.09843968e+00\n",
      "  1.45874442e+00 1.92670598e+00 2.53417257e+00 3.32084902e+00\n",
      "  4.33565280e+00 5.63812564e+00 7.30038579e+00 9.41171716e+00\n",
      "  1.20933523e+01 1.55514520e+01 2.02915174e+01 2.82057597e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.19305424e-01 2.48036687e-01 3.96232143e-01\n",
      "  5.75329916e-01 7.98974476e-01 1.08394890e+00 1.45127089e+00\n",
      "  1.92748145e+00 2.54613627e+00 3.34947732e+00 4.39021830e+00\n",
      "  5.73334741e+00 7.45792876e+00 9.65935592e+00 1.24542114e+01\n",
      "  1.59953099e+01 2.05217989e+01 2.65264530e+01 3.52944061e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.34915608e-01 2.80491949e-01 4.48085447e-01\n",
      "  6.50640795e-01 9.03607953e-01 1.22600459e+00 1.64166917e+00\n",
      "  2.18074054e+00 2.88137889e+00 3.79170496e+00 4.97187708e+00\n",
      "  6.49616588e+00 8.45489875e+00 1.09564636e+01 1.41308145e+01\n",
      "  1.81396972e+01 2.32085938e+01 2.97148989e+01 3.83601321e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.43888596e-01 2.99151321e-01 4.77908602e-01\n",
      "  6.93983062e-01 9.63885059e-01 1.30795845e+00 1.75173827e+00\n",
      "  2.32756345e+00 3.07646681e+00 4.05032270e+00 5.31416240e+00\n",
      "  6.94847041e+00 9.05119064e+00 1.17392423e+01 1.51499412e+01\n",
      "  1.94445683e+01 2.48199523e+01 3.15341725e+01 3.99082992e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.48451813e-01 3.08645547e-01 4.93098639e-01\n",
      "  7.16095789e-01 9.94715871e-01 1.35003105e+00 1.80853626e+00\n",
      "  2.40386058e+00 3.17879308e+00 4.18761262e+00 5.49863016e+00\n",
      "  7.19671769e+00 9.38540688e+00 1.21879377e+01 1.57465598e+01\n",
      "  2.02196211e+01 2.57759541e+01 3.25804532e+01 4.07308175e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.49629830e-01 3.11104227e-01 4.97056003e-01\n",
      "  7.21912775e-01 1.00294441e+00 1.36149261e+00 1.82444743e+00\n",
      "  2.42603210e+00 3.20994238e+00 4.23184616e+00 5.56216297e+00\n",
      "  7.28887789e+00 9.51985673e+00 1.23836709e+01 1.60272700e+01\n",
      "  2.06078796e+01 2.62747061e+01 3.31307810e+01 4.11495814e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.48005144e-01 3.07737601e-01 4.91711735e-01\n",
      "  7.14232637e-01 9.92446417e-01 1.34758069e+00 1.80644669e+00\n",
      "  2.40327568e+00 3.18195301e+00 4.19868010e+00 5.52501122e+00\n",
      "  7.25103157e+00 9.48808468e+00 1.23697912e+01 1.60489491e+01\n",
      "  2.06860575e+01 2.64225239e+01 3.33282772e+01 4.13120780e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.43969373e-01 2.99358901e-01 4.78362176e-01\n",
      "  6.94931353e-01 9.65814290e-01 1.31178782e+00 1.75916517e+00\n",
      "  2.34166230e+00 3.10270704e+00 4.09825419e+00 5.40009983e+00\n",
      "  7.09951319e+00 9.31061472e+00 1.21721256e+01 1.58445910e+01\n",
      "  2.04976012e+01 2.62780295e+01 3.32480462e+01 4.12886805e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.37825061e-01 2.86596080e-01 4.58007465e-01\n",
      "  6.65454319e-01 9.25041445e-01 1.25679190e+00 1.68613923e+00\n",
      "  2.24580042e+00 2.97813657e+00 3.93810467e+00 5.19685950e+00\n",
      "  6.84592050e+00 9.00144731e+00 1.18073159e+01 1.54339153e+01\n",
      "  2.00663325e+01 2.58708502e+01 3.29254865e+01 4.11102909e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.29832497e-01 2.69989261e-01 4.31507467e-01\n",
      "  6.27043335e-01 8.71838679e-01 1.18488745e+00 1.59039145e+00\n",
      "  2.11961387e+00 2.81326162e+00 3.72454549e+00 4.92305949e+00\n",
      "  6.49953460e+00 8.57122705e+00 1.12869107e+01 1.48285822e+01\n",
      "  1.94031299e+01 2.52106559e+01 3.23693421e+01 4.07837872e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.20232933e-01 2.50039087e-01 3.99659270e-01\n",
      "  5.80849681e-01 8.07792311e-01 1.09820228e+00 1.47472464e+00\n",
      "  1.96673920e+00 2.61272709e+00 3.46339518e+00 4.58579349e+00\n",
      "  6.06865694e+00 8.02904578e+00 1.06197603e+01 1.40352836e+01\n",
      "  1.85097733e+01 2.42918205e+01 3.15672958e+01 4.02975362e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.09260403e-01 2.27231370e-01 3.63236788e-01\n",
      "  5.27992144e-01 7.34445818e-01 9.98810242e-01 1.34187661e+00\n",
      "  1.79073852e+00 2.38109715e+00 3.16038832e+00 4.19206187e+00\n",
      "  5.56144250e+00 7.38363895e+00 9.81370746e+00 1.30580133e+01\n",
      "  1.73816488e+01 2.30959293e+01 3.04872239e+01 3.96211284e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 9.71475138e-02 2.02049202e-01 3.23010291e-01\n",
      "  4.69585801e-01 6.53340691e-01 8.88788478e-01 1.19460035e+00\n",
      "  1.59521270e+00 2.12301821e+00 2.82141603e+00 3.74913719e+00\n",
      "  4.98647100e+00 6.64429204e+00 8.87700590e+00 1.19001300e+01\n",
      "  1.60102187e+01 2.15933758e+01 2.90754286e+01 3.87000839e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 8.41277963e-02 1.74977886e-01 2.79754419e-01\n",
      "  4.06753959e-01 5.66033677e-01 7.70242151e-01 1.03569993e+00\n",
      "  1.38385499e+00 1.84330192e+00 2.45266066e+00 3.26479405e+00\n",
      "  4.35315596e+00 5.82158428e+00 7.81965033e+00 1.05665497e+01\n",
      "  1.43859202e+01 1.97451097e+01 2.72524823e+01 3.74439980e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 7.04359730e-02 1.46505513e-01 2.34249160e-01\n",
      "  3.40629732e-01 4.74099326e-01 6.45308716e-01 8.68038032e-01\n",
      "  1.16046513e+00 1.54695319e+00 2.06065437e+00 2.74743569e+00\n",
      "  3.67202332e+00 4.92799903e+00 6.65466703e+00 9.06624381e+00\n",
      "  1.25021047e+01 1.75058906e+01 2.49074293e+01 3.57038311e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 5.63069041e-02 1.17120762e-01 1.87276333e-01\n",
      "  2.72350921e-01 3.79122887e-01 5.16148406e-01 6.94523773e-01\n",
      "  9.28934839e-01 1.23915623e+00 1.65227512e+00 2.20612749e+00\n",
      "  2.95485653e+00 3.97836916e+00 5.39930317e+00 7.41506131e+00\n",
      "  1.03606033e+01 1.48311212e+01 2.18923873e+01 3.32292391e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 4.19736808e-02 8.73089094e-02 1.39613115e-01\n",
      "  2.03050424e-01 2.82686645e-01 3.84925102e-01 5.18086167e-01\n",
      "  6.93211607e-01 9.25226299e-01 1.23468686e+00 1.65053748e+00\n",
      "  2.21468739e+00 2.99008432e+00 4.07594426e+00 5.63865327e+00\n",
      "  7.97928390e+00 1.16904629e+01 1.80236106e+01 2.95873289e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 2.76651777e-02 5.75466831e-02 9.20236688e-02\n",
      "  1.33843735e-01 1.86351856e-01 2.53780516e-01 3.41637227e-01\n",
      "  4.57245591e-01 6.10533067e-01 8.15229563e-01 1.09078389e+00\n",
      "  1.46560974e+00 1.98296632e+00 2.71248159e+00 3.77490079e+00\n",
      "  5.40124962e+00 8.09254950e+00 1.31056528e+01 2.40045079e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 1.36032860e-02 2.82964353e-02 4.52496212e-02\n",
      "  6.58146329e-02 9.16378968e-02 1.24804002e-01 1.68027917e-01\n",
      "  2.24924348e-01 3.00401855e-01 4.01266445e-01 5.37195244e-01\n",
      "  7.22405553e-01 9.78729362e-01 1.34178618e+00 1.87475661e+00\n",
      "  2.70355630e+00 4.12437710e+00 7.02187519e+00 1.50513116e+01\n",
      "  5.00000000e+01]\n",
      " [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n",
      "  0.00000000e+00]]\n",
      "r = 0.0\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.000000\n",
      "0.400000 0.000000\n",
      "0.600000 0.000000\n",
      "0.800000 0.000000\n",
      "1.000000 0.000000\n",
      "1.200000 0.000000\n",
      "1.400000 0.000000\n",
      "1.600000 0.000000\n",
      "1.800000 0.000000\n",
      "2.000000 0.000000\n",
      "2.200000 0.000000\n",
      "2.400000 0.000000\n",
      "2.600000 0.000000\n",
      "2.800000 0.000000\n",
      "3.000000 0.000000\n",
      "3.200000 0.000000\n",
      "3.400000 0.000000\n",
      "3.600000 0.000000\n",
      "3.800000 0.000000\n",
      "4.000000 50.000000\n",
      "r = 0.1\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.090310\n",
      "0.400000 0.187756\n",
      "0.600000 0.299936\n",
      "0.800000 0.435505\n",
      "1.000000 0.604785\n",
      "1.200000 0.820469\n",
      "1.400000 1.098440\n",
      "1.600000 1.458744\n",
      "1.800000 1.926706\n",
      "2.000000 2.534173\n",
      "2.200000 3.320849\n",
      "2.400000 4.335653\n",
      "2.600000 5.638126\n",
      "2.800000 7.300386\n",
      "3.000000 9.411717\n",
      "3.200000 12.093352\n",
      "3.400000 15.551452\n",
      "3.600000 20.291517\n",
      "3.800000 28.205760\n",
      "4.000000 50.000000\n",
      "r = 0.2\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.119305\n",
      "0.400000 0.248037\n",
      "0.600000 0.396232\n",
      "0.800000 0.575330\n",
      "1.000000 0.798974\n",
      "1.200000 1.083949\n",
      "1.400000 1.451271\n",
      "1.600000 1.927481\n",
      "1.800000 2.546136\n",
      "2.000000 3.349477\n",
      "2.200000 4.390218\n",
      "2.400000 5.733347\n",
      "2.600000 7.457929\n",
      "2.800000 9.659356\n",
      "3.000000 12.454211\n",
      "3.200000 15.995310\n",
      "3.400000 20.521799\n",
      "3.600000 26.526453\n",
      "3.800000 35.294406\n",
      "4.000000 50.000000\n",
      "r = 0.30000000000000004\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.134916\n",
      "0.400000 0.280492\n",
      "0.600000 0.448085\n",
      "0.800000 0.650641\n",
      "1.000000 0.903608\n",
      "1.200000 1.226005\n",
      "1.400000 1.641669\n",
      "1.600000 2.180741\n",
      "1.800000 2.881379\n",
      "2.000000 3.791705\n",
      "2.200000 4.971877\n",
      "2.400000 6.496166\n",
      "2.600000 8.454899\n",
      "2.800000 10.956464\n",
      "3.000000 14.130814\n",
      "3.200000 18.139697\n",
      "3.400000 23.208594\n",
      "3.600000 29.714899\n",
      "3.800000 38.360132\n",
      "4.000000 50.000000\n",
      "r = 0.4\n",
      "z\tu_numeric\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.000000 0.000000\n",
      "0.200000 0.143889\n",
      "0.400000 0.299151\n",
      "0.600000 0.477909\n",
      "0.800000 0.693983\n",
      "1.000000 0.963885\n",
      "1.200000 1.307958\n",
      "1.400000 1.751738\n",
      "1.600000 2.327563\n",
      "1.800000 3.076467\n",
      "2.000000 4.050323\n",
      "2.200000 5.314162\n",
      "2.400000 6.948470\n",
      "2.600000 9.051191\n",
      "2.800000 11.739242\n",
      "3.000000 15.149941\n",
      "3.200000 19.444568\n",
      "3.400000 24.819952\n",
      "3.600000 31.534173\n",
      "3.800000 39.908299\n",
      "4.000000 50.000000\n",
      "r = 0.5\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.148452\n",
      "0.400000 0.308646\n",
      "0.600000 0.493099\n",
      "0.800000 0.716096\n",
      "1.000000 0.994716\n",
      "1.200000 1.350031\n",
      "1.400000 1.808536\n",
      "1.600000 2.403861\n",
      "1.800000 3.178793\n",
      "2.000000 4.187613\n",
      "2.200000 5.498630\n",
      "2.400000 7.196718\n",
      "2.600000 9.385407\n",
      "2.800000 12.187938\n",
      "3.000000 15.746560\n",
      "3.200000 20.219621\n",
      "3.400000 25.775954\n",
      "3.600000 32.580453\n",
      "3.800000 40.730817\n",
      "4.000000 50.000000\n",
      "r = 0.6\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.149630\n",
      "0.400000 0.311104\n",
      "0.600000 0.497056\n",
      "0.800000 0.721913\n",
      "1.000000 1.002944\n",
      "1.200000 1.361493\n",
      "1.400000 1.824447\n",
      "1.600000 2.426032\n",
      "1.800000 3.209942\n",
      "2.000000 4.231846\n",
      "2.200000 5.562163\n",
      "2.400000 7.288878\n",
      "2.600000 9.519857\n",
      "2.800000 12.383671\n",
      "3.000000 16.027270\n",
      "3.200000 20.607880\n",
      "3.400000 26.274706\n",
      "3.600000 33.130781\n",
      "3.800000 41.149581\n",
      "4.000000 50.000000\n",
      "r = 0.7\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.148005\n",
      "0.400000 0.307738\n",
      "0.600000 0.491712\n",
      "0.800000 0.714233\n",
      "1.000000 0.992446\n",
      "1.200000 1.347581\n",
      "1.400000 1.806447\n",
      "1.600000 2.403276\n",
      "1.800000 3.181953\n",
      "2.000000 4.198680\n",
      "2.200000 5.525011\n",
      "2.400000 7.251032\n",
      "2.600000 9.488085\n",
      "2.800000 12.369791\n",
      "3.000000 16.048949\n",
      "3.200000 20.686058\n",
      "3.400000 26.422524\n",
      "3.600000 33.328277\n",
      "3.800000 41.312078\n",
      "4.000000 50.000000\n",
      "r = 0.7999999999999999\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.143969\n",
      "0.400000 0.299359\n",
      "0.600000 0.478362\n",
      "0.800000 0.694931\n",
      "1.000000 0.965814\n",
      "1.200000 1.311788\n",
      "1.400000 1.759165\n",
      "1.600000 2.341662\n",
      "1.800000 3.102707\n",
      "2.000000 4.098254\n",
      "2.200000 5.400100\n",
      "2.400000 7.099513\n",
      "2.600000 9.310615\n",
      "2.800000 12.172126\n",
      "3.000000 15.844591\n",
      "3.200000 20.497601\n",
      "3.400000 26.278029\n",
      "3.600000 33.248046\n",
      "3.800000 41.288680\n",
      "4.000000 50.000000\n",
      "r = 0.8999999999999999\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.137825\n",
      "0.400000 0.286596\n",
      "0.600000 0.458007\n",
      "0.800000 0.665454\n",
      "1.000000 0.925041\n",
      "1.200000 1.256792\n",
      "1.400000 1.686139\n",
      "1.600000 2.245800\n",
      "1.800000 2.978137\n",
      "2.000000 3.938105\n",
      "2.200000 5.196859\n",
      "2.400000 6.845920\n",
      "2.600000 9.001447\n",
      "2.800000 11.807316\n",
      "3.000000 15.433915\n",
      "3.200000 20.066332\n",
      "3.400000 25.870850\n",
      "3.600000 32.925487\n",
      "3.800000 41.110291\n",
      "4.000000 50.000000\n",
      "r = 0.9999999999999999\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.129832\n",
      "0.400000 0.269989\n",
      "0.600000 0.431507\n",
      "0.800000 0.627043\n",
      "1.000000 0.871839\n",
      "1.200000 1.184887\n",
      "1.400000 1.590391\n",
      "1.600000 2.119614\n",
      "1.800000 2.813262\n",
      "2.000000 3.724545\n",
      "2.200000 4.923059\n",
      "2.400000 6.499535\n",
      "2.600000 8.571227\n",
      "2.800000 11.286911\n",
      "3.000000 14.828582\n",
      "3.200000 19.403130\n",
      "3.400000 25.210656\n",
      "3.600000 32.369342\n",
      "3.800000 40.783787\n",
      "4.000000 50.000000\n",
      "r = 1.0999999999999999\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.120233\n",
      "0.400000 0.250039\n",
      "0.600000 0.399659\n",
      "0.800000 0.580850\n",
      "1.000000 0.807792\n",
      "1.200000 1.098202\n",
      "1.400000 1.474725\n",
      "1.600000 1.966739\n",
      "1.800000 2.612727\n",
      "2.000000 3.463395\n",
      "2.200000 4.585793\n",
      "2.400000 6.068657\n",
      "2.600000 8.029046\n",
      "2.800000 10.619760\n",
      "3.000000 14.035284\n",
      "3.200000 18.509773\n",
      "3.400000 24.291821\n",
      "3.600000 31.567296\n",
      "3.800000 40.297536\n",
      "4.000000 50.000000\n",
      "r = 1.2\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.109260\n",
      "0.400000 0.227231\n",
      "0.600000 0.363237\n",
      "0.800000 0.527992\n",
      "1.000000 0.734446\n",
      "1.200000 0.998810\n",
      "1.400000 1.341877\n",
      "1.600000 1.790739\n",
      "1.800000 2.381097\n",
      "2.000000 3.160388\n",
      "2.200000 4.192062\n",
      "2.400000 5.561442\n",
      "2.600000 7.383639\n",
      "2.800000 9.813707\n",
      "3.000000 13.058013\n",
      "3.200000 17.381649\n",
      "3.400000 23.095929\n",
      "3.600000 30.487224\n",
      "3.800000 39.621128\n",
      "4.000000 50.000000\n",
      "r = 1.3\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.097148\n",
      "0.400000 0.202049\n",
      "0.600000 0.323010\n",
      "0.800000 0.469586\n",
      "1.000000 0.653341\n",
      "1.200000 0.888788\n",
      "1.400000 1.194600\n",
      "1.600000 1.595213\n",
      "1.800000 2.123018\n",
      "2.000000 2.821416\n",
      "2.200000 3.749137\n",
      "2.400000 4.986471\n",
      "2.600000 6.644292\n",
      "2.800000 8.877006\n",
      "3.000000 11.900130\n",
      "3.200000 16.010219\n",
      "3.400000 21.593376\n",
      "3.600000 29.075429\n",
      "3.800000 38.700084\n",
      "4.000000 50.000000\n",
      "r = 1.4000000000000001\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.084128\n",
      "0.400000 0.174978\n",
      "0.600000 0.279754\n",
      "0.800000 0.406754\n",
      "1.000000 0.566034\n",
      "1.200000 0.770242\n",
      "1.400000 1.035700\n",
      "1.600000 1.383855\n",
      "1.800000 1.843302\n",
      "2.000000 2.452661\n",
      "2.200000 3.264794\n",
      "2.400000 4.353156\n",
      "2.600000 5.821584\n",
      "2.800000 7.819650\n",
      "3.000000 10.566550\n",
      "3.200000 14.385920\n",
      "3.400000 19.745110\n",
      "3.600000 27.252482\n",
      "3.800000 37.443998\n",
      "4.000000 50.000000\n",
      "r = 1.5000000000000002\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.070436\n",
      "0.400000 0.146506\n",
      "0.600000 0.234249\n",
      "0.800000 0.340630\n",
      "1.000000 0.474099\n",
      "1.200000 0.645309\n",
      "1.400000 0.868038\n",
      "1.600000 1.160465\n",
      "1.800000 1.546953\n",
      "2.000000 2.060654\n",
      "2.200000 2.747436\n",
      "2.400000 3.672023\n",
      "2.600000 4.927999\n",
      "2.800000 6.654667\n",
      "3.000000 9.066244\n",
      "3.200000 12.502105\n",
      "3.400000 17.505891\n",
      "3.600000 24.907429\n",
      "3.800000 35.703831\n",
      "4.000000 50.000000\n",
      "r = 1.6000000000000003\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.056307\n",
      "0.400000 0.117121\n",
      "0.600000 0.187276\n",
      "0.800000 0.272351\n",
      "1.000000 0.379123\n",
      "1.200000 0.516148\n",
      "1.400000 0.694524\n",
      "1.600000 0.928935\n",
      "1.800000 1.239156\n",
      "2.000000 1.652275\n",
      "2.200000 2.206127\n",
      "2.400000 2.954857\n",
      "2.600000 3.978369\n",
      "2.800000 5.399303\n",
      "3.000000 7.415061\n",
      "3.200000 10.360603\n",
      "3.400000 14.831121\n",
      "3.600000 21.892387\n",
      "3.800000 33.229239\n",
      "4.000000 50.000000\n",
      "r = 1.7000000000000004\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.041974\n",
      "0.400000 0.087309\n",
      "0.600000 0.139613\n",
      "0.800000 0.203050\n",
      "1.000000 0.282687\n",
      "1.200000 0.384925\n",
      "1.400000 0.518086\n",
      "1.600000 0.693212\n",
      "1.800000 0.925226\n",
      "2.000000 1.234687\n",
      "2.200000 1.650537\n",
      "2.400000 2.214687\n",
      "2.600000 2.990084\n",
      "2.800000 4.075944\n",
      "3.000000 5.638653\n",
      "3.200000 7.979284\n",
      "3.400000 11.690463\n",
      "3.600000 18.023611\n",
      "3.800000 29.587329\n",
      "4.000000 50.000000\n",
      "r = 1.8000000000000005\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.027665\n",
      "0.400000 0.057547\n",
      "0.600000 0.092024\n",
      "0.800000 0.133844\n",
      "1.000000 0.186352\n",
      "1.200000 0.253781\n",
      "1.400000 0.341637\n",
      "1.600000 0.457246\n",
      "1.800000 0.610533\n",
      "2.000000 0.815230\n",
      "2.200000 1.090784\n",
      "2.400000 1.465610\n",
      "2.600000 1.982966\n",
      "2.800000 2.712482\n",
      "3.000000 3.774901\n",
      "3.200000 5.401250\n",
      "3.400000 8.092549\n",
      "3.600000 13.105653\n",
      "3.800000 24.004508\n",
      "4.000000 50.000000\n",
      "r = 1.9000000000000006\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.013603\n",
      "0.400000 0.028296\n",
      "0.600000 0.045250\n",
      "0.800000 0.065815\n",
      "1.000000 0.091638\n",
      "1.200000 0.124804\n",
      "1.400000 0.168028\n",
      "1.600000 0.224924\n",
      "1.800000 0.300402\n",
      "2.000000 0.401266\n",
      "2.200000 0.537195\n",
      "2.400000 0.722406\n",
      "2.600000 0.978729\n",
      "2.800000 1.341786\n",
      "3.000000 1.874757\n",
      "3.200000 2.703556\n",
      "3.400000 4.124377\n",
      "3.600000 7.021875\n",
      "3.800000 15.051312\n",
      "4.000000 50.000000\n",
      "r = 2.0000000000000004\n",
      "z\tu_numeric\n",
      "0.000000 0.000000\n",
      "0.200000 0.000000\n",
      "0.400000 0.000000\n",
      "0.600000 0.000000\n",
      "0.800000 0.000000\n",
      "1.000000 0.000000\n",
      "1.200000 0.000000\n",
      "1.400000 0.000000\n",
      "1.600000 0.000000\n",
      "1.800000 0.000000\n",
      "2.000000 0.000000\n",
      "2.200000 0.000000\n",
      "2.400000 0.000000\n",
      "2.600000 0.000000\n",
      "2.800000 0.000000\n",
      "3.000000 0.000000\n",
      "3.200000 0.000000\n",
      "3.400000 0.000000\n",
      "3.600000 0.000000\n",
      "3.800000 0.000000\n",
      "4.000000 0.000000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFqpJREFUeJzt3X9wVeWdx/HPlySEDGooEiSCbAxjGdtOKm7Gqgw71WzFKUWY0XbYH5XuuMPM/piiXexKnbq0szva0am668502NpZ3LVFR6kQbRfdiFN2dGmD2KgFRDK0Em4alCYgG0JCvvvHPTdNYsI5N8n9kSfv1wxz733uczhfj5xPTp7znHPM3QUAmPymFboAAMDEINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSjN58rmzJnjNTU1+VwlAEx6e/fufd/dq+L65TXQa2pq1NzcnM9VAsCkZ2a/TtKPIRcACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQiETTFs3siKRTks5J6nP3ejObLekpSTWSjkj6krv/LjdlAgDiZHOEfoO7X+Xu9dHneyQ1ufsVkpqizwAw5T23r01LH3hZl9/zgpY+8LKe29eWl/WOZ8hllaQt0fstklaPvxwAmNye29emjdveVFtnt1xSW2e3Nm57My+hnjTQXdKLZrbXzNZFbZe4e0qSote5uSgQACaTB3ceVHfvuSFt3b3n9ODOgzlfd9JL/5e6+zEzmyvpJTM7kHQF0Q+AdZK0cOHCMZQIAJPHsc7urNonUqIjdHc/Fr12SPqxpGsk/dbMqiUpeu0YZdnN7l7v7vVVVbH3lgGASe3SWRVZtU+k2EA3s5lmdmHmvaSbJL0laYektVG3tZK256pIAJgs7l6+WBVlJUPaKspKdPfyxTlfd5Ihl0sk/djMMv1/6O7/ZWa/kPS0md0h6TeSvpi7MgFgcli9ZL6k9Fj6sc5uXTqrQncvXzzQnkvm7jlfSUZ9fb1z+1wAyI6Z7R00ZXxUXCkKAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AApE40M2sxMz2mdnz0efLzWyPmR0ys6fMbHruygQAxMnmCH29pP2DPn9H0sPufoWk30m6YyILAwBkJ1Ggm9kCSSskfT/6bJJulPRM1GWLpNW5KBAAkEzSI/RHJH1dUn/0+WJJne7eF30+Kmn+BNcGAMhCbKCb2Rckdbj73sHNI3T1UZZfZ2bNZtZ8/PjxMZYJAIiT5Ah9qaRbzOyIpK1KD7U8ImmWmZVGfRZIOjbSwu6+2d3r3b2+qqpqAkoGAIwkNtDdfaO7L3D3GklrJL3s7n8maZek26JuayVtz1mVAIBY45mH/veSvmZm7yo9pv74xJQEABiL0vguv+fur0h6JXrfKumaiS8JADAWXCkKAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAEOgAEAgCHQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgYgPdzGaY2c/N7Jdm9raZfStqv9zM9pjZITN7ysym575cAMBokhyh90i60d0/LekqSTeb2bWSviPpYXe/QtLvJN2RuzIBAHFiA93TPow+lkV/XNKNkp6J2rdIWp2TCgEAiSQaQzezEjN7Q1KHpJckHZbU6e59UZejkuaPsuw6M2s2s+bjx49PRM0AgBEkCnR3P+fuV0laIOkaSVeO1G2UZTe7e72711dVVY29UgDAeWU1y8XdOyW9IulaSbPMrDT6aoGkYxNbGgAgG0lmuVSZ2azofYWkP5a0X9IuSbdF3dZK2p6rIgEA8Urju6ha0hYzK1H6B8DT7v68mf1K0lYz+0dJ+yQ9nsM6ASCvnm0/oftbU2rr6dX88jJtrK3WrfNmF7qs84oNdHdvkbRkhPZWpcfTASAoz7af0IaD76m7P31q8GhPrzYcfE+SijrUuVIUAIa5vzU1EOYZ3f2u+1tTBaooGQIdAIZp6+nNqr1YEOgAMMz88rKs2osFgQ4Aw2ysrVbFNBvSVjHNtLG2ukAVJZNklgsATCmZE5/BzXIBgKno1nmziz7Ah2PIBQACQaADQCAIdAAIBIEOAIEg0AEgEAQ6AASCQAeAQBDoABAIAh0AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQBDoABIJAB4BAlBa6AADIhVT7drUefkhnelKaUV6t2kUbVD1vVaHLyqnYI3Qzu8zMdpnZfjN728zWR+2zzewlMzsUvX4s9+UCQLxU+3YdOHCvzvQck+Q603NMBw7cq1T79kKXllNJhlz6JP2du18p6VpJf2Nmn5B0j6Qmd79CUlP0GQAKrvXwQ+rv7x7S1t/frdbDDxWoovyIDXR3T7n769H7U5L2S5ovaZWkLVG3LZJW56pIAMjGmZ5UVu2hyOqkqJnVSFoiaY+kS9w9JaVDX9LcUZZZZ2bNZtZ8/Pjx8VULAAnMKK/Oqj0UiQPdzC6Q9KykO939ZNLl3H2zu9e7e31VVdVYagSArNQu2qBp0yqGtE2bVqHaRRsKVFF+JJrlYmZlSof5k+6+LWr+rZlVu3vKzKoldeSqSADIRmY2y1Sb5RIb6GZmkh6XtN/dvzvoqx2S1kp6IHoN+/QxgEmlet6q4AN8uCRH6EslfVnSm2b2RtT2DaWD/Gkzu0PSbyR9MTclAgCSiA10d/8fSTbK1w0TWw4AYKy49B8AAkGgA0AgCHQACASBDgCBINABIBAEOgAEgkAHgEAQ6AAQCAIdAALBI+gAFK2WlhY1NTWpq6tLlZWVamhoUF1dXaHLKloEOoCi1NLSosbGRvX29kqSurq61NjYKEmE+igYcgFQlJqamgbCPKO3t1dNTU0Fqqj4EegAilJXV1dW7SDQARSpysrKrNpBoAMoUg0NDSorKxvSVlZWpoYG7to9Gk6KAihKmROfzHJJjkAHULTq6uoI8Cww5AIAgSDQASAQBDoABIJAB4BAEOgAEAhmuQDIqdP7OnRy5xGd6+xRyaxyXbS8RjOXzC10WUEi0AHkzOl9Hercdkje2y9JOtfZo85thySJUM8BhlwA5MzJnUcGwjzDe/t1cueRwhQUOAIdQM6c6+zJqh3jQ6ADyJmSWeVZtWN8CHQAOXPR8hpZ2dCYsbJpumh5TWEKChwnRQHkTObEJ7Nc8oNAB5BTM5fMJcDzJHbIxcx+YGYdZvbWoLbZZvaSmR2KXj+W2zIBAHGSjKH/u6Sbh7XdI6nJ3a+Q1BR9BgAUUGygu/vPJJ0Y1rxK0pbo/RZJqye4LgBAlsY6hn6Ju6ckyd1TZsYAGRCw/bt3affWJ3Tqg/d14cVztGzN7bpy2Q2FLgvD5PykqJmtk7ROkhYuXJjr1QGYYPt379KLmx9T39n0xUCn3j+uFzc/JkmEepEZ6zz035pZtSRFrx2jdXT3ze5e7+71VVVVY1wdgELZvfWJgTDP6Dvbo91bnyhQRRjNWAN9h6S10fu1krZPTDkAis2pD97Pqh2Fk2Ta4o8kvSZpsZkdNbM7JD0g6XNmdkjS56LPAAJ04cVzsmpH4cSOobv7n4zyVcME1wKgCC1bc/uQMXRJKp1ermVrbi9gVRgJV4oCOK/MiU9muRQ/Ah1ArCuX3UCATwIEOjBFvLOnXa9tP6wPT/Togtnlum7VIn38M/MKXRYmEIEOTAHv7GnXricPqO9s+ulBH57o0a4nD0gSoR4Q7ocOTAGvbT88EOYZfWf79dr2wwWqCLlAoANTwIcnRn7k22jtmJwIdGAKuGD2yI98G60dkxOBDkwB161apNLpQ3f30unTdN2qRQWqCLnASVFgCsic+GSWS9gIdGAS6WpsVMfDj6gvlVJpdbXm3nWnKleuTLTsxz8zjwAPHIEOTBJdjY1KffM++ZkzkqS+Y8eU+uZ9kpQ41BE2xtCBSaLj4UcGwjzDz5xRx8OPFKgiFBsCHZgk+lKprNox9RDowCRRWl2dVTumHgIdyLMXWl/QTc/cpLotdbrpmZv0QusLiZabe9edshkzhrTZjBmae9eduSgTkxAnRYE8eqH1BW16dZPOnEuPhadOp7Tp1U2SpBW1K867bObE51hnuSB85u55W1l9fb03NzfnbX1AsbnpmZuUOv3RMe/qmdV68bYXC1ARJgMz2+vu9XH9GHIB8qj9dHtW7UA2CHRgLFqelh7+lLRpVvq15elEi82bOfKFPaO1A9kg0IFstTwtNX5V6npPkqdfG7+aKNTXX71eM0qGnticUTJD669en6NiMZUQ6EC2mr4t9XYPbevtTrfHWFG7Qpuu36TqmdUymapnVmvT9ZtiT4gCSTDLBVPWc/va9ODOgzrW2a1LZ1Xo7uWLtXrJ/PgFu45m1z7MitoVBDhygiN0TEnP7WvTxm1vqq2zWy6prbNbG7e9qef2tcUvXLkgu3YgTwh0TGrPtp9Q/atvq3rXG6p/9W09234i0XIP7jyo7t5zQ9q6e8/pwZ0H4xduuE8qqxjaVlaRbgcKiCEXTFrPtp/QhoPvqbs/fS3F0Z5ebTj4niTp1nmzz7vssc7urNqHqPtS+rXp2+lhlsoF6TDPtAMFQqCj4FLt29V6+CGd6UlpRnm1ahdtUPW8VbHL3d+aGgjzjO5+1/2tqdhAv3RWhdpGCO9LZ1WM0HsEdV8iwFF0GHJBQaXat+vAgXt1pueYJNeZnmM6cOBepdq3xy7b1tObVftgdy9frIqykiFtFWUlunv54kR1A8WII3RMiJaWFjU1Namrq0uVlZVqaGhQXV1d7HKthx9Sf//QI+X+/m61Hn4o9ih9fnmZjo4Q3vPLy2LXm5nNMqZZLkCRItADs3/3Lu3e+oROffC+Lrx4jpatuV1XLrsh0bKn93Xo5M4jOtfZo5JZ5bpoeY1mLpkbu1xLS4saGxvV25sO166uLjU2NkpSbKif6Rn5Xt6jtQ+2sbZ6yBi6JFVMM22sTXY72dVL5hPgCAqBXoTe2dM+pof57t+9Sy9ufkx9Z3skSafeP64XNz8mSbGhfnpfhzq3HZL39kuSznX2qHPbIUmKDfWmpqaBMM/o7e1VU1NTbKDPKK+Ohls+2h4nM05+f2tKbT29ml9epo211bHj50CoxjWGbmY3m9lBM3vXzO6ZqKKKxVjvWy2ln/946MYG7b/yEzp0Y4O6oiPWOO/sadeuJw/owxPpUP7wRI92PXlA7+yJv3nT7q1PDIR5Rt/ZHu3e+kTssid3HhkI8wzv7dfJnUdil+3q6sqqfbDaRRs0bdrQE5HTplWodtGG2GWldKg3X/9JpW64Ss3Xf5Iwx5Q25iN0MyuR9K+SPifpqKRfmNkOd//VRBWXMeYr+qT0/TXGML1sPPetHs/DfF/bflh9Z4cGa9/Zfr22/XDsUfqpD97Pqn2wc509WbUPVllZOWJ4V1ZWxi6bGScfyywXAEON5wj9Gknvunuru5+VtFXShO+F47qibxw3UXr09UcHwjzjzLkzevT1R2OXHc/DfDNH5knbB7vw4jlZtQ9WMqs8q/bBGhoaVFY29ERkWVmZGhoaYpeV0qG+dOluNdz4rpYu3U2YA2M0nkCfL+m9QZ+PRm0TalxX9I3jJkrjuW/1eB7me8HskQN0tPbBlq25XaXTh/YrnV6uZWtuj132ouU1srKh/xysbJouWl4Tu2xdXZ1Wrlw5cEReWVmplStXJprlAmDijOekqI3Q9pHHH5nZOknrJGnhwoVZr2RcV/SN4yZK82bOG/HJMknuW11aXa2+Yx890ZfkYb7XrVqkXU8eGDLsUjp9mq5btSh22cyJz7HMcsmc+BzLLBcpHeoEOFBY4wn0o5IuG/R5gaSPpJi7b5a0WUo/gi7blYzrir7KBdFwywjtMdZfvX7IGLqU/L7Vc++6c8gYupT8Yb6ZcfKxzHKR0qGedJricDOXzE0c4ACKz3gC/ReSrjCzyyW1SVoj6U8npKpB7l6+WBu3vTlk2CXxFX0N96XHzAcPuyS8iVLmxOejrz+q9tPtmjdzntZfvT7RbU/H+zDfj39mXuIAB4CMcT0k2sw+L+kRSSWSfuDu/3S+/mN9SHQhZrkAQLFI+pDocQV6tsYa6AAwlSUNdG7OBQCBINABIBAEOgAEgkAHgEAQ6AAQiLzOcjGz45J+PY6/Yo6k+DtN5Vcx1iRRVzaKsSapOOsqxpqk8Ov6A3eviuuU10AfLzNrTjJ1J5+KsSaJurJRjDVJxVlXMdYkUVcGQy4AEAgCHQACMdkCfXOhCxhBMdYkUVc2irEmqTjrKsaaJOqSNMnG0AEAo5tsR+gAgFEUXaDHPXjazMrN7Kno+z1mVlMkdX3FzI6b2RvRn7/MQ00/MLMOM3trlO/NzP45qrnFzK7OdU0J6/qsmXUN2lbx9zMef02XmdkuM9tvZm+b2UdubJ/v7ZWwpkJsqxlm9nMz+2VU17dG6JP3/TBhXXnfD6P1lpjZPjN7foTv8ret3L1o/ih9G97DkmolTZf0S0mfGNbnryV9L3q/RtJTRVLXVyQ9luft9UeSrpb01ijff17ST5V+utS1kvYUSV2flfR8nrdVtaSro/cXSnpnhP+Hed1eCWsqxLYySRdE78sk7ZF07bA+hdgPk9SV9/0wWu/XJP1wpP9X+dxWxXaEnuTB06skbYnePyOpwcxGehxevuvKO3f/maQT5+myStITnva/kmaZWfxz8HJfV965e8rdX4/en5K0Xx99Bm5et1fCmvIu+u//MPpYFv0ZfrIt7/thwrryzswWSFoh6fujdMnbtiq2QE/y4OmBPu7eJ6lL0sVFUJck3Rr9qv6MmV02wvf5lpcHeY/RddGvzj81s0/mc8XRr7xLlD7CG6xg2+s8NUkF2FbREMIbkjokveTuo26rPO6HSeqS8r8fPiLp65L6R/k+b9uq2AI9yYOnEz2ceoIlWWejpBp3r5P03/r9T+RCKsS2SuJ1pS9l/rSkf5H0XL5WbGYXSHpW0p3ufnL41yMskvPtFVNTQbaVu59z96uUflbwNWb2qWFdCrKtEtSV1/3QzL4gqcPd956v2whtOdlWxRboSR48PdDHzEolVSr3v97H1uXuH7h7T/Tx3yT9YY5rSiLRg7zzzd1PZn51dvefSCozszm5Xq+ZlSkdnE+6+7YRuuR9e8XVVKhtNWj9nZJekXTzsK8KsR/G1lWA/XCppFvM7IjSQ7E3mtl/DuuTt21VbIE+8OBpM5uu9AmEHcP67JC0Nnp/m6SXPTrbUMi6ho213qL0eGih7ZB0ezR741pJXe6eKnRRZjYvM4ZoZtco/e/wgxyv0yQ9Lmm/u393lG553V5JairQtqoys1nR+wpJfyzpwLBued8Pk9SV7/3Q3Te6+wJ3r1E6F1529z8f1i1v26o0F3/pWLl7n5n9raSd+v2Dp982s29Lanb3HUrvAP9hZu8q/VNuTZHU9VUzu0VSX1TXV3Jdl5n9SOlZEHPM7Kikf1D6RJHc/XuSfqL0zI13Jf2fpL/IdU0J67pN0l+ZWZ+kbklr8vBDeamkL0t6MxqDlaRvSFo4qK58b68kNRViW1VL2mJmJUr/AHna3Z8v9H6YsK6874cjKdS24kpRAAhEsQ25AADGiEAHgEAQ6AAQCAIdAAJBoANAIAh0AAgEgQ4AgSDQASAQ/w+8YdZSkCXkgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd83370>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def main():\n",
    "    import numpy as np\n",
    "    from numpy import sin , pi, sinh\n",
    "    global aE,aW,aN,aS,aP,U,U1,n,m,p,p1,h,k,r\n",
    "    import seaborn as sns\n",
    "    from numpy import sin , pi, sinh\n",
    "    import matplotlib.pyplot as plt\n",
    "    from mpl_toolkits import mplot3d\n",
    "    n=20\n",
    "    m=20\n",
    "    p=3\n",
    "    r=2.0\n",
    "    z=4\n",
    "    h=r/n\n",
    "    print(\"h=\",h)#delta r\n",
    "    k=z/m\n",
    "    print(\"k=\",k)#delta z\n",
    "    r0=0\n",
    "    z0=0\n",
    "    p1=1000\n",
    "    eps=0.001\n",
    "    R=np.zeros(p1)\n",
    "    R1=np.zeros(p1)\n",
    "    aP=np.zeros((n,m))\n",
    "    aE=np.zeros((n,m))\n",
    "    aW=np.zeros((n,m))\n",
    "    aN=np.zeros((n,m))\n",
    "    aS=np.zeros((n,m))\n",
    "    U=np.zeros((n+1,m+1))\n",
    "    r=np.zeros(n+1) \n",
    "    z=np.zeros(m+1)\n",
    "    U1=np.zeros((n+1,m+1))\n",
    "    R=np.zeros(p1)\n",
    "    r[0]=r0\n",
    "    for i in range(1,n+1):\n",
    "        r[i]=r[i-1]+h\n",
    "    print(\"r=\",r)\n",
    "    z[0]=z0\n",
    "    for j in range(1,m+1):\n",
    "        z[j]=z[j-1]+k\n",
    "    print(\"z=\",z)    \n",
    "    \n",
    "    for i in range(0,n+1):\n",
    "        U[i,0]=0\n",
    "        U[i,m]=50\n",
    "    for j in range(0,m+1):\n",
    "        U[n,j]=0\n",
    "    #print(\"u=\",U)\n",
    "    #for j in range(0,m+1):\n",
    "        #U[n,j]=50\n",
    "    coeff()\n",
    "    #outer iteration loop\n",
    "    for k2 in range(0,p1):\n",
    "        for i in range(0,n+1):\n",
    "            for j in range(0,m+1):\n",
    "                U1[i,j]=U[i,j] #storing the values to check convergence\n",
    "        gauss()\n",
    "        max1=0\n",
    "        flag=0\n",
    "        for i in range(1,n):\n",
    "            for j in range(1,m):\n",
    "                max2=abs(U[i,j]-U1[i,j])\n",
    "                if(max2>max1):\n",
    "                    max1=max2\n",
    "        R[k2]=max1 #storing the residues in an array\n",
    "        R1[k2]=k2+1 #storing the iterations in an array\n",
    "        if(max1<eps):\n",
    "            print(\"max1=\",max1)\n",
    "            flag=1\n",
    "            break\n",
    "    print(\"iteration=\",k2)\n",
    "    #plot of residue vs iterations\n",
    "    plt.figure(figsize=(8,8))\n",
    "    plt.plot(R1[0:k2],R[0:k2],color=\"red\")\n",
    "    plt.xlabel(\"iterations\")\n",
    "    plt.ylabel(\"residue\")\n",
    "    plt.ylim([-0.06,15])\n",
    "    plt.xlim([0,200])\n",
    "    plt.title(\"control volume 20x20\")\n",
    "    plt.show()\n",
    "    if flag==1:\n",
    "        print(\"solution converge \")\n",
    "    else:\n",
    "        print(\"diverge\")\n",
    "    print(U)\n",
    "   \n",
    "    for i in range(0,n+1):\n",
    "        r1=r[i]\n",
    "        print ('r =',r1)\n",
    "        print(\"z\\tu_numeric\")\n",
    "        for j in range(0,m+1):\n",
    "            print('{0:1f} {1:2f}'.format(z[j], U[i,j]))\n",
    "    #numerical plot\n",
    "    for i in range(0,n+1):\n",
    "        for j in range(0,m+1):\n",
    "            if(i==1.0):\n",
    "                plt.scatter(z[j],U[i,j])\n",
    "    plt.show()\n",
    "main()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}