{
"cells": [
{
"cell_type": "markdown",
"id": "c905588b",
"metadata": {},
"source": [
"## Simple example of relaxation method used to solve Laplace's equation\n",
"### From Computational Physics by Mark Newman, Section 9.1\n",
"#### Program `laplace.py` downloaded directly from Computational Physics - Programs and data - University of Michigan\n",
"\n",
"This is a 2-D problem with $V = 1$ on one side, and $V=0$ on other three sides. Griffiths, Section 3.3.1 solves a very similar problem analytically using separation of variables and Fourier series; in Griffiths, the direction \n",
"perpendicular to the $V=!$ side extends to infinity. A simple extension of Newman's code that enlarges the \n",
"grid in one dimension leads to results that agree with Griffiths.\n",
"\n",
"#### M.L.'s notebook `laplace.ipynb` does the same calculation as this notebook, but without Newman's explicit looping.\n",
"Instead of looping, the array `phi` is shifted up one grid step, shifted down one grid step, shifted \n",
"left one grid step, and shifted down one grid step, and the results are summed and then divided by 4.\n",
"\n",
"`for i in range(20000):\n",
" phi_new = (np.roll(phi_old, 1, axis=0) + np.roll(phi_old, -1, axis=0) \\\n",
" + np.roll(phi_old, 1, axis=1) + np.roll(phi_old, -1, axis=1))/4\n",
" phi_new[:, 0] = phi_new[:,m] = 0 # Reset boundaries\n",
" phi_new[n,:] = 0 # Reset boundary\n",
" phi_new[0,:] = 1 # Reset boundary`\n",
" \n",
" This is MUCH faster than Newman's looping procedure!"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "d9fb6e81",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAGgCAYAAAAtsfn1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA11UlEQVR4nO3df5BW1X3H8c8K+vAjyzbosMsGMOvMOhoxiQHLFJ2ANdBpNB3HTkxEI2n6BxYxbLYNQkwrcXBX+YM61UqqkxGnlsHpRFvbSTNsfhS1zFQDkhicwaYhQKzbrQ1lIeIS4PYPx2fuc4Bz99x77o/nPO/XzM54n3vufe5z2X2O9/s953vaoiiKBABAAM4r+wIAAPCFTg0AEAw6NQBAMOjUAADBoFMDAASDTg0AEAw6NQBAMOjUAADBoFMDAASDTg0AEIxSO7XHHntMPT09mjRpkubNm6cXX3yxzMsBADS5iWW98TPPPKO+vj499thjuuaaa/Q3f/M3+v3f/329/vrrmjNnjvXY06dP67/+67/U3t6utra2gq4YAOBLFEU6evSouru7dd55/p6v2soqaLxgwQJ94hOf0ObNm+uvXX755brppps0ODhoPfaXv/ylZs+enfclAgBydujQIc2aNcvb+Up5Ujtx4oR27dqltWvXNry+dOlS7dy584z2Y2NjGhsbq2+/3w8fOnRI06ZNy/diAQDejY6Oavbs2Wpvb/d63lI6tbffflunTp1SZ2dnw+udnZ0aHh4+o/3g4KC+8Y1vnPH6tGnT6NQAoIn5TiGVllOTzvwwURSd9QOuW7dO/f399e33e/iOjo76a/fdd1/9v+fNm9dwfFdXV8N2/P8MJk2a1LBv4sTGW3L++efX/3vChAkN+8w4cHzb/BxmW5d/SFu82XYel/fI8ouVV16TfGlryyszkuW8Lsfa2tr2nT59OvV7mMfG95v7zO1Tp07V//s3v/lNw76TJ082bL/77rv1/z569GjDPvPBZNeuXfX/PtvDiW+ldGoXXXSRJkyYcMaHHxkZOePpTZJqtZpqtVpRlwcAaFKlDOm/4IILNG/ePA0NDTW8PjQ0pIULF5ZxSQCAAJQWfuzv79cXvvAFzZ8/X7/zO7+jxx9/XAcPHtSdd96Z6ny7d+8+574Pf/jDDdsXXXRR/b/NnNzkyZMbtuPhyXgoUjozVBnftoUmpcZQZlJoMr5t2+fadrz7sh6btm1efA4fhlu4LC9FhAVdj41v2/a5to3f73jI0NxnbpshRHM7HnKMhxcl6fjx4w3bo6Oj9f9+++23G/b94he/aNi2fTfnobRO7XOf+5z+93//V/fff7/eeustzZ07V9/5znd08cUXl3VJAIAmV+pAkZUrV2rlypVlXgIAICDEYQAAwSj1Sc2nt956q/7fBw4csLaNx4c/+MEPNuz7wAc+0LA9ZcqU+n9fcMEFDfvMEZnxnJuZbzOnA8S3XaYKuEwNyCtXZ3Jpm3Rs2vP4PDaP8zQbX8Pp8xo+nyWHljYvZm675L6S3ife1mXovZlTM7fjeTNzmH68oIX0XlGM973zzjsN+44dO9awffjw4fp///d//3fDPvP7N/7dXASe1AAAwaBTAwAEg04NABCMYHJq8Zjv//zP/zTsM0thxZlxZnN+xtSpU+v/bc5hM3Nq8ZybOafN3Lbl1GzbSfPf4jmgpJyaraxXljyZ7bymMua/tWpZryqUnSpjPlmWslO2bVtJKnN/Ur7NNvfMtm3uM7/P4tvxnJl0Zk4tPtbg17/+dcO++Lw0SfrVr35V/28zp2Z+/5r5uLzxpAYACAadGgAgGMGEH21Vo+PDT6UzQ4FxZumY+CO7+bhuhjXj4Uhz+L8tHOky/N9Wbsvc7zP86NLWts9X6a68Qn1eV+AtORzpM9yYVymstCFGlzChy7FJIcW04UfbsHxzf5Zh+ua27fvLTLXEh/GbIUNb+NH8fjW/f833yRtPagCAYNCpAQCCQacGAAhGMDk1lxIv8VyYmZOyxcVtSzWY22ZOzZZjsy1hY15j0vB/W04tS/mttMP/8yq3lVdprixtXfg6b7MP08+SNxtv27yG6SeVs3LJqdmG6duWjLHl0MztpOVk4t+bZl7MzKnFt83vV/P717ymvPGkBgAIBp0aACAYdGoAgGAEk1OLx6Ft8y+kxhiwOX/MlutIiovHY91mCS0z9h1/X/MabDk2l5Ja5mcxz1vGnDYbl6VyTM2cY3M5Zyvn0NKWuypq7pmZ+4of61L6ypZDkxq/S5JyavG5aUk5tXhpLDOnZubN4tvm96v5PuZnzRtPagCAYNCpAQCCEWT40Xwktw1dNYfaJ61CHWcLRyaFEOLva+5zCT+6DP833yev4f+2If0uK3e7nDfpWBuX0liEH7NVvbftzxKqtIUNXcKaWYbpu6xQbf4tpg0/JlXeTxt+NEOKZtX++H7zPOb3L+FHAABSolMDAASDTg0AEIwgc2pJ5axspWPM4fW24elZlpqIX6OZ1zPzZLaSWuZnTVtSKymX6GtFbduxPpepKWJJm6JKaFUhj5bXEjEubW15MltOzecK1fG/4yylr1yWk7Eth2UriyU1fr+5lMky95nHxt8n6XrJqQEAkBKdGgAgGHRqAIBgBJNTi8fCk3JqtrkbZk4tnmtymc/kMt/FjDmb1xDfb+a+bNfrM6fmq6SWyzy1LPk3276q5dTKmqeWV07N1xIxtvlweeXUXOaeZcmp2eZy2Za0MrdtOTRzv5kns+XNzPOY2/Hv0KScWl454XPhSQ0AEAw6NQBAMIIJP9rCeS5Vrs0yM7Zwni2slfTInXZVXHNIv9k2fo22Elrmthnqs33WpPuQtqRWlvJbJpew5niPy9LW57FxWUI7vsKPLmFC236XclZZSl+ZbeN/Q7Z95rlchulnWc3aZUi/S5kss218v7nPFgJNGsLvUlbNB57UAADBoFMDAASDTg0AEIxgcmq2Iby2+LUZk7atQm3mqFyGirustmteg0scP36NtnybuZ00pD++3zxvXiW1XHJqvvJvzbRi9tkUsYK1zzxZ2rY+S1/Z/r5clpPJspq1Ladm2zZzX7Yh/rZ8m3lsUvmt+PWanyXp3yZvPKkBAIJBpwYACAadGgAgGC2RU7PNEXFZpsbMqdlKSbnkUHyW8YlvJ5XJSptTs31uczvL0jO+lqnxmVOrYh7tXHzm18rIqdn+LrKUyapiTs0278vMZ9nKZGXJk7ksJ2Obg0dODQAAT+jUAADBIPyYEH6MhxzNx/Us4UeX4clpK4ObIUVb2SyX8KP52WzHlrVKtu08JsKPfsOPLmWziljNOulvJn5slvCjrVxUUeFH23ZSW1v5QMKPAACUgE4NABAMOjUAQDBaIqdmi4tnKV+TlIfyJe3yHGYOzbwP8etNWqbGlquznbesnFoRq2T7kmVF9bzktZp1FXJqWVazjn8fuLRNWvnapUyWy9IztmH6tmtI+l603TNyagAAeEKnBgAIBp0aACAYQebUsiw1YYslJ81ps+WSknJLcVlyFPHcmPm5bXkz8zzmZ7UtaWP7rC4ltUxm7i6vnJpLjm28+1zamvczC19z09Lm0MztLG1t98Wl9JVLW/P33nZNec1Tc1l6Jmk+Wdq5Z7bPItnvr0s+Ng88qQEAgkGnBgAIRjDhx7ikx1+XSttph/QnhRvTDhV3Ca0mDem3lckyt+Pvm6VMlktpMfN6XcpZ+VpR2+QrHFmGtOFGc3+WFapdzuur8n6WMlm2bZfvjrxWvnZpm5Q+sU1BsG2XHW408aQGAAgGnRoAIBh0agCAYASTU3OJzactk2XmqMy2ReTUXMoDudyHpPxbfNv8bC5lsopaJTt+TVnuvW3KQRWXofE1pN+lHFfVV7POsvJ1ljJZ8ZxVUpms+H6fOTWXYfq+ymSVnWPjSQ0AEAw6NQBAMOjUAADBCCanFuczjm+LM2eZnxXnkm9xmU9kss1Tc5n/lvTZ4ud1yS3aymKZ5yqrTFarzFMro0xWUm7GZW6ULU+WJadWxjw1s218ORmX87q0dcmTlb3UjIknNQBAMLx3aoODg7r66qvV3t6uGTNm6KabbtK+ffsa2kRRpPXr16u7u1uTJ0/W4sWLtXfvXt+XAgBoMd7Djzt27NBdd92lq6++WidPntS9996rpUuX6vXXX9fUqVMlSRs3btSmTZu0ZcsWXXrppdqwYYOWLFmiffv2qb29PdX7pl0dOik0ET82KfwY37YN95caQ1VmuMHkq6RWlqHXtvCjbSVslyH9tqkBJpdVEMoKN5YdjswSqs5r5euk4fW2fXmVyXL5PsirSr/LaiAu53Vpa7u/tvvi8j1TBO+d2ne/+92G7SeffFIzZszQrl279MlPflJRFOnhhx/Wvffeq5tvvlmS9NRTT6mzs1Nbt27VihUrfF8SAKBF5J5TO3LkiCRp+vTpkqT9+/dreHhYS5curbep1WpatGiRdu7cedZzjI2NaXR0tOEHAABTrp1aFEXq7+/Xtddeq7lz50qShoeHJUmdnZ0NbTs7O+v7TIODg+ro6Kj/zJ49O8/LBgA0qVyH9K9atUo/+clP9NJLL52xz8w3RFF0zhzEunXr1N/fX98eHR21dmxZckm2FZ/N3JGt9I2Z4zHPG+eSx/E5nNqWd7DlyWz7pPRlspKG/9vyby730CX/ZipiJWwXRax0be5PypPZzptlaHja0ldJbW0ln8oY0u+ySnZS6StbqS6X8oFZvmeKllundvfdd+v555/XCy+8oFmzZtVf7+rqkvTeE9vMmTPrr4+MjJzx9Pa+Wq2mWq2W16UCAALhPfwYRZFWrVqlZ599Vj/4wQ/U09PTsL+np0ddXV0aGhqqv3bixAnt2LFDCxcu9H05AIAW4v1J7a677tLWrVv1j//4j2pvb6/nyTo6OjR58mS1tbWpr69PAwMD6u3tVW9vrwYGBjRlyhQtW7bM9+UAAFqI905t8+bNkqTFixc3vP7kk0/qi1/8oiRpzZo1On78uFauXKnDhw9rwYIF2r59e+o5apK/eWq2JVhsJbQk+zw1W34oqW1aLvchKU8Sv16zrUuezNfSMy7lt8x9tnxclntvm1dn47NMmo3LcjK2980rr1vW0jO28lAubfOap2Zux983S1vbZ7WVKEtqW3aOzXunNp4P0NbWpvXr12v9+vW+3x4A0MKo/QgACEZLVOl3CY/YhukmlcmKt00Ku8XP5WtVbFPSU7N5/XG2UljmebOsVmAL/dmq9mdZJTuvsllZwnu29/AVvqlClX6XFTPMv8W05beavUyW7X1cVhVImr5k+2zNNKSfJzUAQDDo1AAAwaBTAwAEI8icmsklHmyKx5bN4f5m3Nk25NxlKLvL8HSbpM+WV5ks28rXLsP0zeu3fXaX/FteObWkY4uWpYRWXjk182/Gdh6XIecuq1n7KpOVZUi/S/ktl9yXS/4ty7+Ny3do0XhSAwAEg04NABAMOjUAQDCCyamlLeOTVBYn3taMOZt5HJfl0F1yanmxzVOzzUVLmqcW35+UU7Odx2XumUv+LanElm2fr6VnbJptnprL0jMu89SS/oZsbcuYp+aSz8oyT812DS6lxWw5NpfvxaqVyeJJDQAQDDo1AEAwggk/2uRVxsf2OJ9UxslWJsuUV9ks2xBplyH9tqr9LkP6k0KIaSv6u1TpN7nc+5Cq9JdRJstl5WufZbJspaRsYUOX8yaVyXIZ/p92pYCk++vyb8GQfgAACkCnBgAIBp0aACAYQebU8lqZ11YWS7KvDm07Num8tlydGUNPyzxP0urW423rUi7Mttq2ud/cl2XKhEsZMl9LzxS1vNB42yadx9eUGVMVVr52KZNl+z6w5b6SVp32dQ22e5iUq0u7mnXZQ/hNPKkBAIJBpwYACAadGgAgGEHm1JKkXRbeZSmXpDyZbckV27FJuRhbjs0W+06ae2Yrk2XLk7nk1JLuQ/wakuaaxfe7lLrKUkLLpqjSZ3nl2FxKYaXNv9nKYJnX4DOnVkSZLJc5bT6vIW0JM5d5alXDkxoAIBh0agCAYAQTfkwb8sgypN9WSipLyMOUduXrLMP9bffFDNHZQpdJYSyXUGX8XEn3IX5NLuHHvEpoVVFeVftdQoppUwEuJbVchsi7DP9PqrzvMvw/bUgxSxi2qCH9VOkHACAlOjUAQDDo1AAAwQgmp2aTVzzY1zBil+HpSTkKG5fh/ra8mUtJpaTSVy5tbdMgbKtFJ+W+bCW1TGlzbFUY0u+rTFZS6au0OWufZbJcVnx2GdJvKzvla+h9XtMVfP67sfQMAAAFoFMDAASDTg0AEIwgc2pFLaPhK9adNJ8sr3xM/H3NuWY2SffMNvfMdqwtL2buT5ozaJvb56ukVhLbcj1FcFkKx+RrnlrSEky28xS19IzLXC6XJWJccmpp57/5vA++xhqUnWPjSQ0AEAw6NQBAMIIMP5p8DenPK+RhC4cltfUZLrNxKWdkG/5vO9a20rW53yVUabJNB3AZ0p/U1vbvlpe0VfqzhOzTVvA3t11KX7mU1HJZHdrnyte+woQu11tG+LHscKOJJzUAQDDo1AAAwaBTAwAEoyVyajZ5LbHgEvM3mfmXLKWxxitLTsU2dD3pnsU/a1JezDb837adJf9msn3WtHkzl+Oy5C9cSr65nCevnLWvnFpey7P4XKHaZapA2hWqff67VRlPagCAYNCpAQCCQacGAAhGMDm1vGLHLm3jse8spZlc8m1lzH+yzUsz2yblPuI5qqRcl21Om8k2p83lHrr8uyUdWzSXPGmWfS55MpMtT+YyPyvL3DOXuahpy2RlOW+WzxbfLut7seh8HE9qAIBg0KkBAIIRTPjRRdqVeV3KOLmc12QLj2WpvO7Cdv1JoT+X4f+20I+tnFWW4f82SStqj3ef6/vmwWf4Me2Qf5e/A5cq/UUN6Xcpv5Vl+L/t78B2rEuoL8uUibJDii54UgMABINODQAQDDo1AEAwWiKn5mt4cpa2ZlzcVsbJJd9mU0a5Lclt+H/8GpNWyXbJqbmsfD3e98yiqFWwfeVcXXIoWf5m0g7pT/r7yiv/Zit9ZbumvIb0J11DEWWyqpZf40kNABAMOjUAQDDo1AAAwQgyp+Zzfk4RZbJMWeapxduePHmyYd/EiY3/3LYcm0vpIzNPZitRZXIpk+XS9lzHJbXNUkLLpqi8g8v7ZMk12/bZfkezLD3jUs7KZd5X2jJZ5vWZf295Xa/L3D5fZbJMVc6x8aQGAAgGnRoAIBhBhh9d+Bq67NrWpfxWnEvIywwL2sKRLsP9zeu1hXOShum7lLOKfx6XyvtZhulnWem6mctk+VoJ2yXklfQ76BKqLKKivy3caB6bZUi/KW2ZrLzaVg1PagCAYNCpAQCCQacGAAhGS+TUfA1dzTIk1oyTpx2m79LWzAfYcmzmPpPLcjJxLqtkJw29t+VcbMdmWfna9vtRds4sq7TDsn3+Hbic16VMlu09fZXJcinVlWVIv0tbX/lNU5bv0KLxpAYACEbundrg4KDa2trU19dXfy2KIq1fv17d3d2aPHmyFi9erL179+Z9KQCAwOXaqb3yyit6/PHH9dGPfrTh9Y0bN2rTpk169NFH9corr6irq0tLlizR0aNH87wcAEDgcsupHTt2TLfddpueeOIJbdiwof56FEV6+OGHde+99+rmm2+WJD311FPq7OzU1q1btWLFirwuaVx8xaRtOTSfbW3M/JUtx+ZzWRpbnsx8H5fSV7Y8mQuX/JvJV9msMmQpbZQ2T2ay5YCy/H255N9ccmq2PJnPpWdc7q+vts2UJ3OR25PaXXfdpRtuuEGf+tSnGl7fv3+/hoeHtXTp0vprtVpNixYt0s6dO896rrGxMY2Ojjb8AABgyuVJbdu2bdq9e7deeeWVM/YNDw9Lkjo7Oxte7+zs1IEDB856vsHBQX3jG9/wf6EAgKB479QOHTqk1atXa/v27Zo0adI5250txHau8M26devU399f3x4dHdXs2bPPOP5s/5207bNt/NHfFkJMamtbJTuJy5D++P6kUKWNbdh+0vB/W5kslzChLXTpskJCkixls9K2tfFVlT/LPcnr78tX2C3pd9llJem8hvS7XIOv1cKL+l4sOpTpvVPbtWuXRkZGNG/evPprp06d0gsvvKBHH31U+/btk/TeE9vMmTPrbUZGRs54entfrVZTrVbzfakAgMB4z6ldf/31eu2117Rnz576z/z583Xbbbdpz549uuSSS9TV1aWhoaH6MSdOnNCOHTu0cOFC35cDAGgh3p/U2tvbNXfu3IbXpk6dqgsvvLD+el9fnwYGBtTb26ve3l4NDAxoypQpWrZsme/LAQC0kFLKZK1Zs0bHjx/XypUrdfjwYS1YsEDbt29Xe3t7GZdjVUSZrKSh7C5L0bjkyWxD+m3HJuXJ4m2Tlp6xrZKddF9s1+CrTJbt/lZ9CH+SZiuTVcSQfpc8WZb8W1KOzXbeMspkNZO2qAk/yejoqDo6Ohpei3+pmV/Itu2kti7ntbU1v3Rt12C2jW/bziM1ftHa9pn7bfvM/bbrS2prGwySpfMpqlMb775mQKeW3DavTs08Nm39Sdt5zP1Jc/DStrXtSzqvJB05ckTTpk074/W0qP0IAAgGnRoAIBgtsfSMySXObNtXVEktlzltafNkLvk3l+VvkthKarmUybJdk21uWdJ5bW1d9lVBXmWyTGlDjC5zrpJ+B13Cj3mFFF3Om/Z6fX4njfc8SceWjSc1AEAw6NQAAMFoifBjllJCZZTUcgkpTpzY+E9oKztluw8uocqkcJ6NLWyYtEq2y4jGOJepAaaqlcVKUkbZLF9hraRwnst7pl0J2yWk6BKyNeUVqszr+8umaqFIntQAAMGgUwMABINODQAQjJbIqdlkiYvnFb+25cKKWkbF15B+20rXZtuk3OJ4z5N03rRlsZKObWZl/R34Gv6f1zI1LlNxXFbUNpXxPZPlGqqMJzUAQDDo1AAAwaBTAwAEI8icWl6x47zKzLjkxXzm1OKS5p6lnaeWpa0tF5b0uW3V/00uZbJs1+e6P29Z/g5sbV3O65JLypJ/c1nKxWWemkupLpecWpbz5lX6arz7sh6bN57UAADBoFMDAAQjyPCjiyzVyfMaamur2p9UzspWfitLqDLtkH6fbV0WFHX5rGlLYWUJL/oKTWYJ9aQNIfksoeUrVJnXkH5T2pCiS/mtvFazzlJ5v+yQogue1AAAwaBTAwAEg04NABCMlsip+Vp6xuXYospkZVlOxmW1aNsq2VlWnba19bW6dVL+Lf7ZsuS6sizJU4QsOVVfQ8XTlsVKOtZlSH+zLVOTV57M9j4sPQMAQAXQqQEAgtES4UeTr8fsoobaugzTtw17L6PyflFD+osYwp+0v2oV/X2GkIoY4l+FIf0+h/+7DNMvYzXrotIyReNJDQAQDDo1AEAw6NQAAMFoyZyaTRklaZLapi2TZcoyTN+Wh8rS1nZcXkP681otPOl9i+aS9/A53L+MIf0uubosbfMa/j/e6zO3i2rbTHhSAwAEg04NABAMOjUAQDDIqXlUxDw1M6eTZf5Q/FxZcnVZ8nrNNk/N10rjVRPaPDXbe7rks9Lm0MxzJd3fIuaptQqe1AAAwaBTAwAEI5jwY9WG3tuuz/W8tn0uZbJsK2q7lNRKaps2TBjakH7bNeQlbbiplYb0247NUvrKpW0ZQ++rMFWgCDypAQCCQacGAAgGnRoAIBjB5NSKkCXObFuxOktOzZS2TJa5z6WtLWeVlINiSH81NMOQ/vGexzxXllydy4raLmWy8spR5bX0TDPhSQ0AEAw6NQBAMOjUAADBIKdmsMWdk/ItvvJkWXJqaXNUWdqmXWrGPNa8vz7nuLlckw3z1Nze0yV/lXRsESW1fOXQbNeTtF3WHDGWngEAoGLo1AAAwWj58KOvUlfmsT7bjnefZB9O71ImyxZKKav0VdrrTVrZwCU02KpD+l3apg0x+iqLleXYLKWvTHkN6R/vebK0TTq2ynhSAwAEg04NABAMOjUAQDBaMqeWdpi+rdRVlra293TZl9TWtmRMUpksX6WvfC0RU9YQftv9L2rYflouvztVH+JfVP7N17D9oob0264pSym0LG2LxpMaACAYdGoAgGDQqQEAgtESObUsc9Fs88lclpNJeh8f+1zbupTUSpsnM2XJk6UtAca8NHe+fs/yWk7G5RqSzpNXrq6MnFoVcvdl40kNABAMOjUAQDBaIvxo4xJSTBqmnzZUmXTe8e4z+RxOn7YSf1ErX7tc03jfQypmCL/LeXyGeooY4u8Shi8rVJlX9f8ihun72pf12CrhSQ0AEAw6NQBAMOjUAADBaPmcmqnqebIssW2XY205tqLyZL6Wl0nKWblcr4vxTnXIK1+R13SELLmZola+Hu95ks5V9WH6VZg6VDU8qQEAgkGnBgAIRi6d2ptvvqnbb79dF154oaZMmaKPf/zj2rVrV31/FEVav369uru7NXnyZC1evFh79+7N41IAAC3Ee07t8OHDuuaaa3TdddfpX/7lXzRjxgz953/+p37rt36r3mbjxo3atGmTtmzZoksvvVQbNmzQkiVLtG/fPrW3t2e+hiy5L1/v41KiymdbmyzlrVzKZNnur0uuznZ/XfJtpiz5N5ss11SEvPKxvpapSbq+LCWrfJ3HltfLkidzuQZfbbOocj6uLfJ8BWvXrtW//du/6cUXXzzr/iiK1N3drb6+Pt1zzz2SpLGxMXV2duqhhx7SihUrzjhmbGxMY2Nj9e3R0VHNnj27oY3tC89WUzCp3qDLQIYqtLXdB1/nNblM6k57HvPYLJ2Py/W6CGk9NZdj6dSS2+Z13rLaxvdnmYQuSUeOHNG0adPOeD0t7+HH559/XvPnz9dnP/tZzZgxQ1dddZWeeOKJ+v79+/dreHhYS5curb9Wq9W0aNEi7dy586znHBwcVEdHR/3H7NAAAJByCD/+/Oc/1+bNm9Xf36+vfe1revnll/XlL39ZtVpNd9xxh4aHhyVJnZ2dDcd1dnbqwIEDZz3nunXr1N/fX98+25NaWj5DlWnbmlymA2QR/z+qLE+AWcpkjfc85rmS/t1s+3yuhB2X9t+prDJZcT5DpWlDUy4lqpL2NfMwfZfrTeKrbdkhRRfeO7XTp09r/vz5GhgYkCRdddVV2rt3rzZv3qw77rij3u5sX9zn+uOu1Wqq1Wq+LxUAEBjv4ceZM2fqIx/5SMNrl19+uQ4ePChJ6urqkqT6E9v7RkZGznh6AwDAhfdO7ZprrtG+ffsaXnvjjTd08cUXS5J6enrU1dWloaGh+v4TJ05ox44dWrhwoe/LAQC0EO/hx6985StauHChBgYGdMstt+jll1/W448/rscff1zSe2HHvr4+DQwMqLe3V729vRoYGNCUKVO0bNky35fjVV5D74vKoRUlbZmsJGmH7bvk30xVH6afF1+5mCzlrGz7fZbJqlpOzeU+FJWbaybeO7Wrr75azz33nNatW6f7779fPT09evjhh3XbbbfV26xZs0bHjx/XypUrdfjwYS1YsEDbt2/3MkcNANC6vM9TK8Lo6Kg6OjoaXks7Ty1LW1/zvpLmyrlcQ9XaJnFZ+DPtebOcp+pzz/LCk1p129rmiOU1r85l7llw89QAAChLSyw945JTcZk/5uvYJnxYTi3LkjZpz5vEZe5ZqE9uLr+DviqIuF6Dy9wz274qPn2lzZOZfB2b11y4IvCkBgAIBp0aACAYLRF+NLmUszrXcWdTRoixCiWWbO+T16rTLhX9be/p+r4usoRTfShq+kFeIUWX9/Q1tN3n0HubMq7B55D+qoUc43hSAwAEg04NABAMOjUAQDBaMqcW53PIdpbSWGW8ZxXi4rb8m+0eZSm/5bNslk0zldTKawi3z+H/vkpA+cpnZTmvL1muN8v7VBlPagCAYNCpAQCCQacGAAhGy+fUTFnKZPl6n2aTpQyZr7a+5p4lzS0r6vcjD1l+53zlB7PMafO1PItL2yrMf8srT2YqIgdYBJ7UAADBoFMDAASD8KODZn4kP5u0K3U3uyxTBdIqqmRWGdMIfE0HaIbwo21fFSv6+6r+30x4UgMABINODQAQDDo1AEAwyKk1mVDj4FLyEH5bSa28pgqY0g7br3rJrLzKZLnks3ztM/c3Q04N/vCkBgAIBp0aACAYdGoAgGCQU0tQtdh3Uk4n7fVWcVmaLGz5N19tTSGXySoj51aFnJpL2yx5syrMU2v2v/n38aQGAAgGnRoAIBjBhB/zenQuI6SU5bNUrfRVluHzzXbeKoQmq7DScZawYdrzVCFUmdfw/zI0c1iTJzUAQDDo1AAAwaBTAwAEI5icWtX4yotlOW+zD9P3lTfLK/+W5RqqrmpD+JPaViGn5tLW12ctO39VRTypAQCCQacGAAgGnRoAIBjk1BxkyWfldd68ji0qH+dSoqoMVZh7lpcy5rT5yn2Vdd60bYvKvxV13irjSQ0AEAw6NQBAMAg/JvAVHssSxkp7bFJ4oYiSWmWEMc/2vmnb+nrPKqhauNG1ra+q9y7n9RWazOtYn39PoYQjeVIDAASDTg0AEAw6NQBAMMipGcpYeqQKua8smr0cly9Vy7EVlb8sqm3Vlq3xufRM2rY+j/VxXBXwpAYACAadGgAgGHRqAIBgtHxOzWcOLW3ezGWeWpa2WcpiVT2vl0VepbqKKLFV1L9FEbmZsnJqaa+hqOvNcg1p38PnsUXjSQ0AEAw6NQBAMFoy/Jg23JRXKSafQ/ptbas49D7ksKZNs33WKoQfXdpWoaK/bV9RQ/p9tfV5bN54UgMABINODQAQDDo1AEAwWjKn5sLXcPq82maJbTdb/q0IVSt1VZaQl6mpWrktn9dUVNsq40kNABAMOjUAQDDo1AAAwWiJnFpR88vSts1r7lleebFWLqnVKpo9p+bSNqQ5bWmvp6i2ReBJDQAQDDo1AEAwWiL86IvP0GQZ1fTLCFWiORF+dG9bher/zRw29IUnNQBAMOjUAADB8N6pnTx5Ul//+tfV09OjyZMn65JLLtH999+v06dP19tEUaT169eru7tbkydP1uLFi7V3717flwIAaDHec2oPPfSQvvnNb+qpp57SFVdcoR/96Ef6oz/6I3V0dGj16tWSpI0bN2rTpk3asmWLLr30Um3YsEFLlizRvn371N7e7vuSrHyWRXLJUbnkvsZ7nrPt98XXMH2f+cMitErZrLzubVHLmxSVZ6ra8Hqf+cIizlOEtsjz1d54443q7OzUt771rfprf/iHf6gpU6bob//2bxVFkbq7u9XX16d77rlHkjQ2NqbOzk499NBDWrFixRnnHBsb09jYWH17dHRUs2fPPuc1JH3x2Pab+1zaVu0asrxnWe+Ttq3PY8s4b9no1MbXNqSaklW4Bkk6cuSIpk2bZm3jwnv48dprr9X3v/99vfHGG5KkH//4x3rppZf06U9/WpK0f/9+DQ8Pa+nSpfVjarWaFi1apJ07d571nIODg+ro6Kj/2Do0AEDr8h5+vOeee3TkyBFddtllmjBhgk6dOqUHHnhAt956qyRpeHhYktTZ2dlwXGdnpw4cOHDWc65bt079/f317aQnNQBAa/LeqT3zzDN6+umntXXrVl1xxRXas2eP+vr61N3dreXLl9fbnS1fcq7QTq1WU61W832p41JEnszlGlzkNafN9X1QfSGHH/NqW1So0kUV8m9l896pffWrX9XatWv1+c9/XpJ05ZVX6sCBAxocHNTy5cvV1dUl6b0ntpkzZ9aPGxkZOePpDQAAF95zau+8847OO6/xtBMmTKgP6e/p6VFXV5eGhobq+0+cOKEdO3Zo4cKFvi8HANBCvD+pfeYzn9EDDzygOXPm6IorrtCrr76qTZs26Utf+pKk90JVfX19GhgYUG9vr3p7ezUwMKApU6Zo2bJlvi+nstKGNW3nMdvmNfyf0GR4mj38mOXYMkY0pj2P67lakfdO7ZFHHtGf//mfa+XKlRoZGVF3d7dWrFihv/iLv6i3WbNmjY4fP66VK1fq8OHDWrBggbZv3174HDUAQFi8z1MrwujoqDo6Os653+ccsaqf1+f7FH2eqh5bxnnLxpNac5wny/tU8bxSE8xTAwCgLEEuPeNzyZW8lnIpYjh9UUva+JIl/1bWKt+havYnt7Tv4zOfVYUnwKo91RWBJzUAQDDo1AAAwaBTAwAEI8icWhXklX9zeZ+k2HZeebIqXAOyqWJOLe15Qmtb5XxWFfCkBgAIBp0aACAYhB8T+AoTjvc9xvM+422bdD2+yluVMeydIfzFqlo4supTA5KOLeoaXM4Tyu86T2oAgGDQqQEAgkGnBgAIRsvn1PJaoTqvYfoubV2G0yed13ZcUeWsGP5fniLudzMUP057XBWH/7topr83ntQAAMGgUwMABINODQAQjJbPqbnwmRcrY45YGfmsKi4nY3ufZsodlKWoe1T1PFmWY5utbTPhSQ0AEAw6NQBAMAg/FqSMIf1Jbcd7XJZjqx5uzPN9WkUzhCN9naeZh/+3Cp7UAADBoFMDAASDTg0AEIyWyKlVIZ9lO7asJWLGe1yW93Q5Fzm0MFS9pFZe56rCkja+jmvmvB5PagCAYNCpAQCCQacGAAhGS+TUskib+0qbX8t6rO08WfJmeZzH9Vxp36dqMf9Wxpy26l5DKH8nPKkBAIJBpwYACAbhRwe+SlJlOdZXaDLLNSRhmD7Gq9mG/+d13qqvKtBMeFIDAASDTg0AEAw6NQBAMFoyp5Z26L3tPK7n8jX83yav4fQM088uy+9dXLPfB5tmG/5f9fNWIQdYBJ7UAADBoFMDAASDTg0AEIyWzKnFZcmL2c5VVp6saufN8r5VjttLfvOxZajCNaRV1rVXLU9W1Hmb6XeFJzUAQDDo1AAAwWj58KPJVzgyS4mqLOdN+x5FhRTzfF8fmvna89TMnzW0UGUR79HM/948qQEAgkGnBgAIBp0aACAY5NQS+CqpZTtvkiLybz6voQrx+CrcMzTf/azC9YacAywCT2oAgGDQqQEAgkGnBgAIBjk1B3nNPcv6vudSVg7Ql2YvQ4Vsmv3fsGrXX7XryQtPagCAYNCpAQCCQfjRoyKG6bsoK9zAqs6okpB/j0L+bGnxpAYACAadGgAgGHRqAIBgkFMrSRGx8CLydmdDnB/g76AsPKkBAIJBpwYACAadGgAgGOTUAkZMH0Cr4UkNABAM507thRde0Gc+8xl1d3erra1N//AP/9CwP4oirV+/Xt3d3Zo8ebIWL16svXv3NrQZGxvT3XffrYsuukhTp07VH/zBH+iXv/xlpg8CAIBzp/brX/9aH/vYx/Too4+edf/GjRu1adMmPfroo3rllVfU1dWlJUuW6OjRo/U2fX19eu6557Rt2za99NJLOnbsmG688UadOnUq/ScBACDKQFL03HPP1bdPnz4ddXV1RQ8++GD9tXfffTfq6OiIvvnNb0ZRFEX/93//F51//vnRtm3b6m3efPPN6Lzzzou++93vnvV93n333ejIkSP1n0OHDkWS+OGHH374afKfI0eOZOmGzuA1p7Z//34NDw9r6dKl9ddqtZoWLVqknTt3SpJ27dql3/zmNw1turu7NXfu3Hob0+DgoDo6Ouo/s2fP9nnZAIBAeO3UhoeHJUmdnZ0Nr3d2dtb3DQ8P64ILLtAHP/jBc7YxrVu3TkeOHKn/HDx40OdlAwBKEnkepZ3LkH6zPFMURYklm2xtarWaarVafXt0dDT7RQIASnf06FF1dHR4O5/XTq2rq0vSe09jM2fOrL8+MjJSf3rr6urSiRMndPjw4YantZGRES1cuHBc79Pd3a1Dhw4piiLNmTNHhw4d0rRp0zx+knCMjo5q9uzZ3KME3Kdk3KNk3KNk79+jgwcPqq2tTd3d3V7P77VT6+npUVdXl4aGhnTVVVdJkk6cOKEdO3booYcekiTNmzdP559/voaGhnTLLbdIkt566y399Kc/1caNG8f1Puedd55mzZpVf2KbNm0av0AJuEfjw31Kxj1Kxj1K1tHRkcs9cu7Ujh07pp/97Gf17f3792vPnj2aPn265syZo76+Pg0MDKi3t1e9vb0aGBjQlClTtGzZMknvfZA//uM/1p/+6Z/qwgsv1PTp0/Vnf/ZnuvLKK/WpT33K3ycDALQc507tRz/6ka677rr6dn9/vyRp+fLl2rJli9asWaPjx49r5cqVOnz4sBYsWKDt27ervb29fsxf/uVfauLEibrlllt0/PhxXX/99dqyZYsmTJjg4SMBAFqVc6e2ePFi62iVtrY2rV+/XuvXrz9nm0mTJumRRx7RI4884vr2DWq1mu67776GQSRoxD0aH+5TMu5RMu5RsrzvUVvkezwlAAAloaAxACAYdGoAgGDQqQEAgkGnBgAIBp0aACAYTd2pPfbYY+rp6dGkSZM0b948vfjii2VfUikGBwd19dVXq729XTNmzNBNN92kffv2NbSJxrF4aysZHBxUW1ub+vr66q9xj97z5ptv6vbbb9eFF16oKVOm6OMf/7h27dpV39/q9+nkyZP6+te/rp6eHk2ePFmXXHKJ7r//fp0+fbreptXuUaUWj/a6kE2Btm3bFp1//vnRE088Eb3++uvR6tWro6lTp0YHDhwo+9IK93u/93vRk08+Gf30pz+N9uzZE91www3RnDlzomPHjtXbPPjgg1F7e3v07W9/O3rttdeiz33uc9HMmTOj0dHREq+8HC+//HL04Q9/OProRz8arV69uv469yiKfvWrX0UXX3xx9MUvfjH693//92j//v3R9773vehnP/tZvU2r36cNGzZEF154YfTP//zP0f79+6O///u/jz7wgQ9EDz/8cL1Nq92j73znO9G9994bffvb346kxnU2o2h89+POO++MPvShD0VDQ0PR7t27o+uuuy762Mc+Fp08edLpWpq2U/vt3/7t6M4772x47bLLLovWrl1b0hVVx8jISCQp2rFjRxRF41u8tVUcPXo06u3tjYaGhqJFixbVOzXu0Xvuueee6Nprrz3nfu5TFN1www3Rl770pYbXbr755uj222+Pooh7ZHZqeS0efS5NGX48ceKEdu3a1bDQqCQtXbr0nAuNtpIjR45IkqZPny5pfIu3toq77rpLN9xwwxl1RrlH73n++ec1f/58ffazn9WMGTN01VVX6Yknnqjv5z5J1157rb7//e/rjTfekCT9+Mc/1ksvvaRPf/rTkrhHprwWjz6XXNZTy9vbb7+tU6dOWRcjbVVRFKm/v1/XXnut5s6dK8m+eOuBAwcKv8aybNu2Tbt379Yrr7xyxj7u0Xt+/vOfa/Pmzerv79fXvvY1vfzyy/ryl7+sWq2mO+64g/sk6Z577tGRI0d02WWXacKECTp16pQeeOAB3XrrrZL4XTKN536kWTz6XJqyU3tfmsVIQ7dq1Sr95Cc/0UsvvXTGvla+X4cOHdLq1au1fft2TZo06ZztWvkeSdLp06c1f/58DQwMSJKuuuoq7d27V5s3b9Ydd9xRb9fK9+mZZ57R008/ra1bt+qKK67Qnj171NfXp+7ubi1fvrzerpXv0dn4Xjz6XJoy/HjRRRdpwoQJZ/Tg8cVIW9Hdd9+t559/Xj/84Q81a9as+uvxxVvjWul+7dq1SyMjI5o3b54mTpyoiRMnaseOHfqrv/orTZw4sX4fWvkeSdLMmTP1kY98pOG1yy+/XAcPHpTE75IkffWrX9XatWv1+c9/XldeeaW+8IUv6Ctf+YoGBwclcY9M47kf8cWjz9VmvJqyU7vgggs0b948DQ0NNbw+NDQ07tWzQxJFkVatWqVnn31WP/jBD9TT09OwP7546/veX7y1Ve7X9ddfr9dee0179uyp/8yfP1+33Xab9uzZo0suuaTl75EkXXPNNWdMB3njjTd08cUXS+J3SZLeeecdnXde41fnhAkT6kP6uUeNxnM/4otHv+/9xaOd71mq4S0V8P6Q/m9961vR66+/HvX19UVTp06NfvGLX5R9aYX7kz/5k6ijoyP613/91+itt96q/7zzzjv1Ng8++GDU0dERPfvss9Frr70W3XrrrUEPMR6P+OjHKOIeRdF70x0mTpwYPfDAA9F//Md/RH/3d38XTZkyJXr66afrbVr9Pi1fvjz60Ic+VB/S/+yzz0YXXXRRtGbNmnqbVrtHR48ejV599dXo1VdfjSRFmzZtil599dX6FKvx3I8777wzmjVrVvS9730v2r17d/S7v/u7rTWkP4qi6K//+q+jiy++OLrggguiT3ziE/Uh7K1G0ll/nnzyyXqb06dPR/fdd1/U1dUV1Wq16JOf/GT02muvlXfRFWB2atyj9/zTP/1TNHfu3KhWq0WXXXZZ9Pjjjzfsb/X7NDo6Gq1evTqaM2dONGnSpOiSSy6J7r333mhsbKzeptXu0Q9/+MOzfgctX748iqLx3Y/jx49Hq1atiqZPnx5Nnjw5uvHGG6ODBw86XwvrqQEAgtGUOTUAAM6GTg0AEAw6NQBAMOjUAADBoFMDAASDTg0AEAw6NQBAMOjUAADBoFMDAASDTg0AEAw6NQBAMP4fsrn3GG19ylQAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from numpy import empty,zeros,max\n",
"from pylab import imshow,gray,show\n",
"\n",
"# Constants\n",
"M = 100 # Grid squares on a side\n",
"V = 1.0 # Voltage at top wall\n",
"target = 1e-6 # Target accuracy\n",
"\n",
"# Create arrays to hold potential values\n",
"phi = zeros([M+1,M+1],float)\n",
"phi[0,:] = V\n",
"phiprime = empty([M+1,M+1],float)\n",
"\n",
"# Main loop\n",
"delta = 1.0\n",
"count = 0\n",
"while delta>target:\n",
" count += 1\n",
" # Calculate new values of the potential\n",
" for i in range(M+1):\n",
" for j in range(M+1):\n",
" if i==0 or i==M or j==0 or j==M:\n",
" phiprime[i,j] = phi[i,j]\n",
" else:\n",
" phiprime[i,j] = (phi[i+1,j] + phi[i-1,j] \\\n",
" + phi[i,j+1] + phi[i,j-1])/4\n",
"\n",
" # Calculate maximum difference from old values\n",
" delta = max(abs(phi-phiprime))\n",
"\n",
" # Swap the two arrays around\n",
" phi,phiprime = phiprime,phi\n",
"\n",
"# Make a plot\n",
"imshow(phi)\n",
"gray()\n",
"show()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9dcd97bc",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}