{ "cells": [ { "cell_type": "markdown", "id": "KypQvVV84qNl", "metadata": { "id": "KypQvVV84qNl" }, "source": [ "# Simulating Molecular Motors walking down a microtubule" ] }, { "cell_type": "markdown", "id": "0Xaj-CEi4qNn", "metadata": { "id": "0Xaj-CEi4qNn" }, "source": [ "In lecture, we discussed the evolution of the positional distribution function of a molecular motor that processively walks down a biopolymer (e.g. Kinesin V walking down a microtubule). In the simplest models, this distribution function satisfies a Smoluchowski equation of the following form. The probability density to observe the walker at position $x$ at time $t$ is given by" ] }, { "cell_type": "markdown", "id": "YUz7uHkn4qNn", "metadata": { "id": "YUz7uHkn4qNn" }, "source": [ "$$ \\partial_t p(x,t)= D \\partial_x^2 p - v \\partial_x p $$" ] }, { "cell_type": "markdown", "id": "r9GLs43z4qNn", "metadata": { "id": "r9GLs43z4qNn" }, "source": [ "where $D$ is the diffusivity of the walker and $v$ is its velocity." ] }, { "cell_type": "markdown", "id": "c1lXmJB44qNo", "metadata": { "id": "c1lXmJB44qNo" }, "source": [ "Here, we'd like to directly time integrate the above equation to see how the pdf evolves over time. " ] }, { "cell_type": "markdown", "id": "DoIQFR8v4qNo", "metadata": { "id": "DoIQFR8v4qNo" }, "source": [ "Initializing the simulation" ] }, { "cell_type": "code", "execution_count": null, "id": "ajVRcc634qNo", "metadata": { "id": "ajVRcc634qNo" }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# Diffusivity of the random walkers\n", "D = .2\n", "# mean speed of the random walkers\n", "v = 1.\n", "# total length of microtubule (a.u.)\n", "w = 40.\n", "# simulation step size \n", "dx = 0.01\n", "nx = int(w/dx)\n", "# time step of the simulation\n", "dx2 = dx*dx\n", "dt = dx2 / (2 * D )\n", "\n", "\n", "\n", "\n", "uini = np.zeros(nx)\n", "u = uini.copy()\n", "\n", "x = np.linspace(0,w,nx)\n", "\n", "# Initial conditions - ensemble of random walkers sharply peaked at (cx\n", "r, cx = .2, 10\n", "r2 = r**2\n", "for i in range(nx):\n", " p2 = (i*dx-cx)**2 \n", " if p2 < r2:\n", " uini[i] = 1/(2*r)\n" ] }, { "cell_type": "markdown", "id": "a8dbbUx44qNp", "metadata": { "id": "a8dbbUx44qNp" }, "source": [ "The core simulation algorithm." ] }, { "cell_type": "code", "execution_count": null, "id": "aE-WBpJN4qNp", "metadata": { "id": "aE-WBpJN4qNp" }, "outputs": [], "source": [ "def do_timestep(u0, u):\n", " # Propagate with forward-difference in time, central-difference in space\n", " u[1:-1] = u0[1:-1] + D * dt * (\n", " (u0[2:] - 2*u0[1:-1] + u0[:-2])/dx2 ) - v * dt * (\n", " (u0[2:] - u0[:-2] )/(2*dx) )\n", "\n", " u0 = u.copy()\n", " return u0, u\n" ] }, { "cell_type": "markdown", "id": "IMqBXxPd4qNq", "metadata": { "id": "IMqBXxPd4qNq" }, "source": [ "Running the simulations and plotting the output at regular time intervals" ] }, { "cell_type": "code", "execution_count": null, "id": "C83h97bL4qNq", "metadata": { "id": "C83h97bL4qNq", "outputId": "4c8b90fa-de78-4349-bf2b-12fd8524eda1" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1f3/8ddn1ux7QhIS9rAquyyiVMUF3K0L1n2pS1tta1vb2m9b236/ra39tbWtVuuudccNVCwiIILKTljDEgJIQiA72ZeZOb8/7hBiCBAkZODO5/l4xJk5987MJzfh7c2Zc88RYwxKKaVOfo5QF6CUUqpraKArpZRNaKArpZRNaKArpZRNaKArpZRNuEL1xikpKaZPnz6henullDoprVy5sswYk9rRtpAFep8+fVixYkWo3l4ppU5KIrLzUNu0y0UppWxCA10ppWxCA10ppWxCA10ppWxCA10ppWziiIEuItkiskBE8kRkg4j8oIN9zhKRfSKSG/z69fEpVyml1KF0ZtiiD/ixMWaViMQCK0VkrjFmY7v9FhljLu76EpVSSnXGEQPdGFMMFAfv14hIHtATaB/o6gSTu6uK+Xl7rQciXDoigwFpsaEtSil13BzVhUUi0gcYBSztYPNEEVkD7AZ+YozZ0MHz7wTuBOjVq9fR1qqO0j/nbWXephJEwBgorWnioW+eGuqylFLHSac/FBWRGOAt4IfGmOp2m1cBvY0xI4B/Au929BrGmCeNMWONMWNTUzu8clV1Ib8xjMiKZ/tDF5ERH0EgoIuZKGVnnQp0EXFjhfnLxpi32283xlQbY2qD92cDbhFJ6dJKlVJKHVZnRrkI8AyQZ4z56yH2SQ/uh4iMC75ueVcWqo5e+9UFDXqGrpSddaYPfRJwI7BORHKDbb8AegEYY54ArgK+IyI+oAG41uhipUop1a06M8plMSBH2OdR4NGuKkp1IesPp8P/AJVStqBXioYR/ZtJKXvTQLcxzW+lwosGus3t72oJfmatlLIxDXSllLIJDXQbaz/QSLtglLI3DXSllLIJDXSb065zpcKHBnoY0WGLStmbBrpSStmEBrrNHRi2GNIylFLdQAM9jOjkXErZmwa6jWmfuVLhRQPd5vZfIapdLkrZnwa6UkrZhAa6jR3UZ65dMErZmga6UkrZhAa6zbUOW9QlLpSyPQ30MKI9LkrZmwa6jemwRaXCiwa6ze0frqjDFpWyPw30MNJ+fnSllL1ooNuY5rdS4UUDXSmlbEID3eb2D1fULnSl7E8D3cbaXymqPTBK2ZsGulJK2YQGut21DlvUThel7E4DPYzoqBel7E0D3cY0wJUKLxroNiftbpVS9qWBHkb0hF0peztioItItogsEJE8EdkgIj/oYB8RkX+ISL6IrBWR0cenXHU0NMCVCi+uTuzjA35sjFklIrHAShGZa4zZ2GafaUBO8Gs88HjwVimlVDc54hm6MabYGLMqeL8GyAN6ttvtMuBFY1kCJIhIRpdXq46aaCe6UmHjqPrQRaQPMApY2m5TT2BXm8eFHBz6iMidIrJCRFaUlpYeXaXq6LVfUlSHvShla50OdBGJAd4CfmiMqW6/uYOnHJQexpgnjTFjjTFjU1NTj65SpZRSh9WpQBcRN1aYv2yMebuDXQqB7DaPs4Ddx16eOlY6OZdS4aMzo1wEeAbIM8b89RC7zQJuCo52mQDsM8YUd2Gdqgtoh4tS9taZUS6TgBuBdSKSG2z7BdALwBjzBDAbuBDIB+qBW7u+VHW0rNkW9dxcqXBxxEA3xizmCKlgrE/bvtdVRSmllDp6eqWozUnb2Ra1z0UpW9NAtzEdpahUeNFAV0opm9BAt7nWLpfQlqGU6gYa6GGk/RqjSil70UC3MY1vpcKLBrrNtV4pqn0uStmeBnoY0VEvStmbBrqN6eyKSoUXDXSllLIJDXSbOzBsUbTLRSmb00C3Mc1vpcKLBrpSStmEBnqY0GGLStmfBnoY0StFlbI3DXQb0w9BlQovGuhKKWUTGug2J206z/WMXSl700C3Mc1vpcKLBrpSStmEBrrN7e9wERE9Y1fK5jTQlVLKJjTQ7Uw/BVUqrGig25yuKapU+NBADyN6wq6UvWmg25jmt1LhRQNdKaVsQgPd5g4MWwQ9Z1fK3jTQbUz7zJUKLxroSillE0cMdBF5VkRKRGT9IbafJSL7RCQ3+PXrri9TfV37J+fSBS6Usj9XJ/Z5HngUePEw+ywyxlzcJRWp40a7YJSytyOeoRtjPgUquqEW1cV0hSKlwktX9aFPFJE1IvKhiAw71E4icqeIrBCRFaWlpV301koppaBrAn0V0NsYMwL4J/DuoXY0xjxpjBlrjBmbmpraBW+tjqR12CI626JSdnfMgW6MqTbG1AbvzwbcIpJyzJWpY6Z95kqFl2MOdBFJl+BQChEZF3zN8mN9XaWUUkfniKNcRORV4CwgRUQKgQcBN4Ax5gngKuA7IuIDGoBrjdFzwxNF62yLAvpjUcrejhjoxphvHWH7o1jDGpVSSoWQXilqY3pCrlR40UC3PWnzX6WUnWmghxE9YVfK3jTQbUwDXKnwooGulFI2oYFuc9JmhQv9kFQpe9NAtzEdd65UeNFAV0opm9BAtzlpc6vn60rZmwa6UkrZhAa6UkrZhAa6UkrZhAa6zelsi0qFDw10G9P8Viq8aKArpZRNaKDbnOhsi0qFDQ10pZSyCQ10GzN6KZFSYUUD3eYOjHLRybmUsjsNdKWUsgkNdBvTM3KlwosGulJK2YQGus219qGjH5IqZXca6Dam8a1UeNFAV0opm9BAt7nWK0VFPyRVyu400JVSyiY00G1Mp8tVKrxooIcRzXel7E0D3e5ahy3qfItK2Z0Guo3pCblS4eWIgS4iz4pIiYisP8R2EZF/iEi+iKwVkdFdX6ZSSqkj6cwZ+vPA1MNsnwbkBL/uBB4/9rJUV5E2d/RKUaXs7YiBboz5FKg4zC6XAS8ayxIgQUQyuqpApZRSndMVfeg9gV1tHhcG2w4iIneKyAoRWVFaWtoFb60OS0/IlQorXRHoHQ2f6DBKjDFPGmPGGmPGpqamdsFbqyMRObCmqA5bVMreuiLQC4HsNo+zgN1d8LpKKaWOQlcE+izgpuBolwnAPmNMcRe8rjpGekKuVHhxHWkHEXkVOAtIEZFC4EHADWCMeQKYDVwI5AP1wK3Hq1h1bDTglbK3Iwa6MeZbR9hugO91WUWqS+3/gENnW1TK/vRKURvTybmUCi8a6EopZRMa6DYnbSfn0hN2pWxNA10ppWxCA93G9IRcqfCigR5GdHIupexNA93m2g5bVErZmwa6jemoRaXCiwZ6GNGAV8reNNBtrnW2Re1yUcr2NNCVUsomNNBtTEe1KBVeNNBtrnWUC6LxrpTNaaCHK79PPyVVymY00G3skHld8Ak8lAWv36ChrpSNaKCHkdbpdD/6JfgaYNP7sHl2aItSSnUZDXS72z/b4v7O9IZK2LMOBl1obcx9JVSVKaW62BFXLFInrw57U3Z8Zt0Onw7Vu2HLnG6tSSl1/OgZehgxAAULrAf9zoKs0yDQAnXloStKKdVlNNBtTmh3ieiupRCVApEJ0GeS1Va0svsLU0p1OQ30cFNfAQm9rPvZE6zbHYtCV49SqstooIcRr2mE6iLIOc9qiE23bhsqQleUUqrLaKDbXOuaoiL0bPnSehCdemBjz7GQPz80xSmlupQGehgZ6Nts3ckae6DRE2WNSVdKnfQ00G3MtBu3GBeotu4k5xxo7DPZGpveWN2NlSmljgcN9DAyuCXPGuHijTnQGJVo3ZZvDU1RSqkuo4Fuc9LmNtrUgrT7kWeOtm5Lt3RnWUqp40ADPYykBEqh/zlfbUzsY91Wbu/2epRSXUsD3cba9qB7Ao0kBSohIv6rO0UlgTva6kdXSp3UNNBtbv+wxcRA8PL+pH4H7xSXYU2pq5Q6qXUq0EVkqohsFpF8Efl5B9tvEZFSEckNfn2760tVxyKrZYd1JyXn4I2eaAyGOTvmMGPLDFr8Ld1am1KqaxxxtkURcQKPAecBhcByEZlljNnYbtfXjTH3HIca1dfUdtRirL/KuhOfffCOWeP4v4I3eWPhTwD4eOfHPDblMVwOnYxTqZNJZ87QxwH5xpgCY0wz8Bpw2fEtS3W1JH+pdWf/h6BtrPbX8kZcLGemj+cnY3/C57s/54k1T3RvgUqpY9aZQO8J7GrzuDDY1t6VIrJWRN4UkQ5OA0FE7hSRFSKyorS09GuUq47W/tkWc5o20IQXXJ6D9nnCvxeAXw28npuG3oRLXLyU9xIBE+jWWpVSx6YzgS4dtLVfOuE9oI8xZjjwMfBCRy9kjHnSGDPWGDM2NTX16CpVR820+TGJMVQ4kg7ap6a5hs+r8zm9voGMqiJEhHtG3UNdSx2fFX3WneUqpY5RZwK9EGh7xp0F7G67gzGm3BjTFHz4FDCma8pTXWVA00a2u/oe1P5+wfsAXFVT23r5/1UDrwJgxpYZ3VegUuqYdSbQlwM5ItJXRDzAtcCstjuISEabh5cCeV1XojoWIkDAj5sWAh38sTW7wFokekpDE5RZV4vGe+MZkjSEBbsWHDQfjFLqxHXEQDfG+IB7gDlYQf2GMWaDiPxORC4N7vZ9EdkgImuA7wO3HK+C1ddQb813vt516leaW/wt5JbmMiFjAg6nByoKWrdN6TUFgHVl67qvTqXUMenUuDRjzGxgdru2X7e5/wDwQNeWpo5V68l1mTVtbrUj9ivbl+5ZCsDZ2WfDzgJoqW/dNrXvVN75+J9U334vW/cJEUOGkP6bB3Gnp3dL7Uqpo6dXitqcCFBbAkCJI+0r2+Z9OQ+A83qfB0l9oXB567b08gB/eMFP9M5SosaOpXbhQnZcMx1fuS4ordSJSgM9HAT7xnc7M7/SvGT3EuI8caRGpYIzOJyxpRFjDLvuupvIZvjttQ7S/vwQWY//C19JCUX3/ai7q1dKdZIGuo21fpxZZS09VysH5kFv9jdTWFvI+IzxVkN28LaulH1vv03Ll19SPu00tmUKK/asIPass4g971zqly2jfuXK7vsmlFKdpoEeDhr3sdeViR9na9OKvSsA+EbWN6yGyAQAzJ51lD76GAA5v/49AMv2LAOgxy9/BcCW/3uMNfN2sXnpHlqa/d3yLSiljkwn67A9gZ2fUeHq/ZXWhbsWAnB65ulWQ+oQAGo/X46vuJikm2+mR2I2SRFJLCxcyA9G/wBnairbzvoxO+kHM6wVjiLfdDP1zlPIzEnsvm9JKdUhPUMPBwE/fnF+pWnZnmVEOCOs/nOAOKt/vfyNOQAk330XAIMSB7G1civ+gJ8F/8ljJ/3IKP6MC2IWcPl9oxCHMPORXEp26pqkSoWaBrqNGQOuQBM0VbPVe2rrMEZjDPlV+YxIG3FgZ28MviahYUclUaedhivROuMelzEOgLmzVrHpiz0MnpjOiMByWt5/k8z+sVx5/xgcTuGN3y8mf+VKijbnUVNe1t3fqlIKDXTb69G0A4BmR2Rr25ZKa9RLa3dLUFWJNVd6ynfubm2bnDWZqOY4tv23huh4D2fdMJjEq62pAareeYfYJC/ZA/Np2vckMx9+kNd+fT9Pfe82/vuvv9FUX49SqvtoH7rNxfisq0QLvEOg0Wr7tPBTACb3nPyVfStyGxCnIWrixNa2nIQcztl2AwDn3joUp9NBwvTp7P3DQ1S+9TafbM+jYOUyohIG0+IbypnX5FC1ZxOrPpzF7q2bmf7gH4n0RCORLkQ6mudNKdVVNNBtzdCrbj0AVc7k1mGMS4utK0T7J/Rv3bOpoAB/nZ/4vg2IMa1r15V+WUNW1SCK4/LJGmwtMO3weomePJnPt+dRtLKG0y69khHnX8tLv/qCvC/cXPebO8hKG0zZ7M2UPbQap7iQSBeRg5OIO7cXruQDfy0opbqOBrrNeQINAJS7erD/FH1j+UYGJAz4yhlzxQsvApA0qBb27YJEa1TM4jes0SwL+r1KWcM3SYlMAWDnoL4U7SsiJ6svk6+/FYABY9LIX1lC4aubiFwLPeMHsq18FREZCQwacjoNa0upW1uCmZJCS28PbrebtLQ0IiIiuuVYKGV3Gug217duDST1QxxWeJc1lFHTUsNlGV9ddKr6v//FER1BRIIPKrdDYm/KCmso3raPuIFCdWQZy/csZ1rfaZR+uYNlyxYT09jM0K07W1/jzGtziNlYDmtKiRiSRNLVA9n4/EpWLXyT6EkZNF2QyOIFC9n3aV3rc5xOJ4MHD+acc84hOTm5ew6KUjalgW5jxoA70AgOb2vbosJFAHwj+xutbQ3rNxDYt4+kay8H/gXFa6DfWXz2Zj4A500fycMfw4IvFzC1z1Rm/vl/ATg9ow8tnyzCV1mJKzGR5oVF9PY62NkUoP/kLBxRbs659S42LPmcd+fMJRAVQ0ZGBiMbBhK7V/BMSKPIU8KuXe/w8bwnSE31ExUVQUREFklJk0jvcRkez8GLciilOqaBbnNpTTthxF1QZg1X/Gy3tQrRqLRRrftUvWktZJF40+3w4r+gtoTq8gYKN1WSNTiR9J6JpEens2LvCla+/w77SvYy4rwL6TNgKIWfLKLy5VeImnAltYuLcA1IIHdFKaXvFnD5faOorK6hedBIAs3NZHmE2+64AwnA3idyKd85h+jh79B/QDk+XxxlZXHEx3kxcTtZXL6RtdvWURV9HvFxQ4h0uhkYHcGUpDh6eN0hOZZKneg00G0szQRnRnRHti5tsWrvKnpE9cDrtM7ajTFUv/8BrswMPP0GWItIF69h6UxrbvTxl/YDrP8BfLzlvyx861lcbg9n33IHDqf161Pz8VKaiobjSoukxy3DyNyXS9HmSrZt/pLX33oJRMgONFC1ZiPFWzeT1q8nxac/SkXVQiKq+jFk7MMkZZzJ+++/z7vrd7DplLHkeaLBQExNDc7aQholhiZjjbO9JC2Bn/ZNp3+U9r0r1ZaOQ7ex3hRZd/qcCUCABkobSg9MyAU0bthIoLaW+IsuthqcXporStmybC/JPaNJ7xcPWOPRJ2ywuj8u+M4PcLrciAix50/F3f8GcAgpt56CuBycfuUA/I5GXn39JYwx3HbbbVx1rzVL48fP/Z7lyy+jomohvZPvodeSX+J/PYaGgGHRoFG8M/obbMfFFS3VLJ0whM9OqeUJuZPnuZmZg33clZ3K3PJqzlq2mT8VFOML6IpKSu2ngW5jg812605CLwCa3JsAOCf7nNZ9Kp591tpl+jVWQ58zWFk0FoAx0/q07jfcmUP/3TEEkiMZPOlA/7s752rEHYkrdjuuROuMOSHDS03qenz+Zr41/XrS09OJT0tn+NTx9DhjGQ2NhQw/9QkGjLiPuCl92FrfxDmfbuCl4gpu75nCb6t30ePz+exZm0uPtGmMHfs2TqeXuk3Xc0/Sdr4YP4QLUuL42869XLJqK3ubWo7L8VPqZKOBbmMDCI5ASR6ACDS41wAwIXMCEOxumT0bd3Y2nqwsAAIRyayqu5LIGCcDxhxYEGPRk/8G4ItR+1rb6teV0lLkx1+xldq5L7S+5uuvv06L1BNXNYTiNT4AGhp24e0/E6cnQP57vUhKPBuA1aMSuXFiFCWBAM9kpPP7gVnc8M0ryMjI4MMPP6SgoIDYmMGMO+1dPJ4U1qy5DXfdcp4+pS9/HJjFmpp6zl2xmY21DcfxSCp1ctBAt7EMU0ajIxoc1o+5yZ1HZnQm0e5oAGrnWSsWxV92YAjjuspJAIwcUds6Tn177kr25G+hqV8sWz3F1DTXEKhvoeIV64zfnVFMS1ERzTt2MH/+fAoKCpg8eTIxpJM790saGvayYuU1+Py1xPm/S+3uSBa/9iLv7q3k2nUFJHrdvPRFHaNm7MD4DU6nkxtvvJGoqChee+01qqqqiIjIZOyYN3G7k8hdczv79uVyS88U3ho1gAZ/gItWbmFJVW23HVulTkQa6DY2ho1si7W6T5rNPgKOWiZmHrisv3KGNbol6aYbAevsesUXggMfI3ptam378NG/AHDq9VcDMP/L+ZS/thkMJE0fRNJ1Vvvyf/+bRYsWMWDAAM4++2wmXN4fh7ueJZ9fQ3NzCSNGPM24qfcRFZ/As1t3cvfGneREeflo/CBOnZhNoK6Fqve2ARAVFcUNN9xAc3MzL7/8Mj6fj4iIDMaMfh2nM4LVuTdRV1fAxIQYZo8ZSLTTydW521hYUdMNR1apE5MGul01VAJQ67I+yCzDutz/vN7nARCoq6Nu4adEjRuHMy4OgM1L9tBY72dE/Fyc2xcAsHzWWzTUVDPusqs4d/A0AIqWbKFpSyXeAQlEjUojYtBA6vv0YZ7LRXx8PFdffTUiwtAzU8me/AgBKWTY0L+TnHQGIkLjt3/CnLOuIKuuig/GDCTV4ybu3F44kyOoW1JMU3Aq3szMTC677DJKS0uZOXMmAFFRvRk96mWM8bFq9bdoaiphYHQEH44dSIrHxbfWbOPjcp3KV4UnDXS72mWtMrQtzpr+ttxYKxTtnw63/JlnAEi87jrAOhNfEhyqeFrGp1Czh4baGha98jxubwSTpt9IpCuSkTHDuXT9BBBIvtFaFKO5uZlPJ52OEeGy3r3xer0YY9iw8V4ik7ezd/W1lG0eCcBzRWX8qQ6ya8q58tV/UJ1v/SUgIqTedgoAZc+tx/gCAIwaNYrRo0ezbt06li2zvqfY2KGMGP40zc1lrFp9HT5fDdkRHj4YnUOG181NawuYW3agr1+pcKGBbldfLgFgT2QOzf5matiKp2UgbocbY0xrd0vsBecDsGXpHuqqmhhxbjbu/hOhfCtzHrO6WqZ+94c4nE6MMdxXeD0AVZd6cXhdGGOYMWMG+/x+Jn7xBeappwHYtPl/KCubR3bWHVTmT2HxjHye2lXCA1sKGRkbxXvjh+HxNfP+3/6ICU7U7kqOJP7ifphGP1UfFLR+KxdddBEZGRnMnj2bHTt2AJCUdDqnDPs79fXbWZ17C4FAE5kRHj4YM5DMCDc3rdvORxrqKsxooNtVmTXneZUnnbk75wIQ2XQaALULFuAvLSP5jm8jIgQChgUvbwZg/CX9IHMUO2oT2LZqBZkDhzBwwhnW8z7fTVpxLLMSP+GNwPsAzJs3j61btzJp0iSGDBpE8/btbF7+ALt3v056+uUMHPhzxl/Sl0X9Pfwqfzdj4qJ4e9QA0nukM/aSb1JXVcnnb7zUWnbsGT2JGJRI3RfFNORZF0Y5nU6uv/56IiMjeeWVV6istLqTevS4mIE5v6K6Opc1a+/CGD/pXjezRw8kO8LDTeu2815JVTccbKVODBrodrXpAxYzEhHhjc1vABDZPAaA0kcfBSD59tsBWPHBdvwtAcZf1g+314mv5wRm7x4MwMX3/QyA5l017HuvAHdmNP/NWcY7+e+wcuVKFi9ezIABAzj33HNJ+/GPqZnqp7DmDVKSz2HokD8D8PEgL/NGRpFd2sIrg3oT5bR+7c687mYi4+JZ8vbrlBftai096bohOGI9VLyyCV+FNUNkTEwMN954Iz6fjxdeeIGGBmuYYnb2LfTtcy8VFYtYt+57GOMnzevmgzE59Iv0cseGHbxaXH5cD7VSJwoNdDuqKAAM+fSiKVDNqpJVJDAMwUPdkqU0bcwj7uKLcSYkUFfVxPIPdhAR7Wb0BdaUufNnfkSD382UIX5ik1LwVzdT9tx6xOsk5ZZhXNj/IjLqMnjvvfdIS0tj+vTpiAhFvE/NpX48m4RBCQ8QQPjRpi955MsSJrq83PhJDV+8vKW1TIfDyTd//hsAPvj7wwQCfqvd6yTl1mEYX4CyZ9YRaLbaMzMzueaaa6iqquKFF16gubkZgH79fkiv7NspLZvL2nXfJRDwkeqxQn1YTAT3bdrFP3bu7b7jr1SIaKDb0YZ3AVjIGLY2fQBApuM8TMBQ/D//A0Da/fcD8PHzGwE4//ZhOBxC3qIFrJs/h8FZLkbyOYGGZkqfWkug3kfKrcNwxnkZ7xrPxJKJtHhbuPnmm3G5nGze8ju2b3+EBO8okh9zkf/bPzA9dxuvFFdwbXoSb54xmMzecWxbXUpx/oFukPT+OZx53S2U7tzOguefam33ZMaQNH0QvvJGyp7bgPFb/eyDBw/mkksuYc+ePbz44outoZ6T8wt69fo2ZWUfk7vmFvz+BhLdLmaNzuHMxBj+UFDMvXk7aQ4EjvPBVyp0NNDtKG8WAKvpz5amWSRFJJEkIzlzwye0FBWReN11uHuksWFREYWbKhkwNo3soUkUbc5j9mN/JTGjJ+dfcQEBE0Hpo5/jK20g6dpBePvEs3btWmbOmInxGj5K+4gv67eydu1dFBa+QFrahYya+BrFF13B9RdOZ3FVLb/ol8EjQ3rhFGHa3acC8OG/19HU4Gstd9xlVzFw/CRy57zPmrmzW9ujRqYRd35vmrfvo/w/G1tDfcyYMUybNo3CwkKeffZZ6oNrl+YMeIB+/X5EZeUXLF9xBY2Nu4l2OnltRH9uykxmxp5KLl65ld2Nzd31k1CqW2mg201DFexeDQOn4U+yPgy9e8TdRNdVcemydxCvl7Sf/ZTd+VV88vJmkntGc+6tQ9mzbStv/t8viYiJ5apf/i+OgVdS2vxnWsqdJFw+gIjhKcybN4+3336b1NRUpt80nYTIejauuZ6y8vn07v0dBg/5O48XVnDbhdOpi4zk4Rf+xT2psa2lRcd7ueTeETTUtPDBo2sI+A+cLU+79yek9e3Px0//i42fzm9tjzunF7HfyKJxUwVlz61v7X4ZP35865n6E088wZ49ewDo2+d7DB36F+rrt7N02cWUly/EKcLDg7L5y6Bs8uoaOXv5Zt7cU9E6ukYpu9BAt5uFDwNQMvxKiF9AhCRyTd9vcukrD+P1NZP9739TUtTIe3/PJSLazcX3jGRH7gpee/CnON0upj/4EN66CPY+VUCL6U2i+6/4+gV46aWXWLRoETk5Odx2241ENs/j/vRm3KYRR8b32JN0F1NXbeF/t+3mjMQ4Pmwq57Qliyi89/uYNt0cvYYlM/GK/hRv28fsx9fhD4a6y+3m6l/9nuSsXnz42F9ZNvPN1ufET+tL3AW9acqvouTR1Zy5KNcAAAvTSURBVLSUWmfkY8aM4YYbbqCxsZGnn36aZcuWEQgEyEi/nNPGvo3bnUDumtvYmPczmpvLuT4zmfmnDaJvpJd78r7kmjXb2KBzwCgbkVCdpYwdO9asWLEiJO9tW8318IdMmjwxXH/KBDZXbmaK635++MEcGlat4tUzr+eGO77NvOfz8ES5uPT7w8hbPIvls94iLiWVK3/8GxzrW6j9fDeOGA9REypZtegRPpfx4HQzZcoZZGUXsnPn4zQ17SEu6Vzu2FpLZewlNHkHk+l186v+mVyeloCIUPLII5Q/8W9izzuXnn/5C+LxtJa6ZOY2Vn64k/R+8Uy98xSiE6z52Zvq63n34d9RmLeenPGnc86tdxOTaF3t2rC+jIoZWzB+Q9zZ2cSc2ROHx0lVVRXvvvsuO3bsICMjgylTptC/f38CgUYKCv7Gl7uew+mMIjv7ZrKzbsbhSuK5ojL+34497PP5uSg1nruyUjktPvor66wqdSISkZXGmLEdbutMoIvIVODvgBN42hjzx3bbvcCLwBigHJhujNlxuNfUQO9ixsCMmynb9B4/GDqetbWFDFs7iR8t3UJs+R4WTrmLoojhpNdBQg8HvYeVsGH++9RWVjBq3IWM6HcOjbkV+JtbqBrsYEdsOevyNuD3N3BK/Gr6j+tBVctq6n2NFEVdyNboa5lbE8XuphYc/kpSGz7lxYnTGZEy+Ctl7f3Tw1Q89xwRp5xC+oO/JvLUU1u3rV9YyKIZW3G5nYyZ2pthk3vijXQRCPhZ8tZrLH1nBk6XixHnX8jI8y8kPi0d374mqmZuo3FjOY5oN9Hj04ka3QNnkpe1a9cyf/58qqur6dGjB6NGjWLIkCE4nXsp2P43Sks/wuHwkJJyLmlpFyIx43l2TxPPFJZS4w8wIMrLJakJnJscx/DYKNwODXd14jmmQBcRJ7AFOA8oBJYD3zLGbGyzz3eB4caYu0XkWuAKY8z0w72uBnrXaGmoYte2/7JjyQvk78xjnS+ZtPIEJpb1wVcVoDKhB409B9BU34RHmkiMbaSlvhaHM4K4xEw8kRE0U0WDt4qq+GpKnfU0ugyBaHDFu2hwNlNNNHtNOrsdg9kuPWkxgkeEyUmxXJ2eSKo/n/s/+SF1LXVM6zuN83ufz9DkoaRGpeIQB/tmzWLvHx7CX1VF5NgxxJ4zhcjhp+Lu1Yvqlgg+e6uAXRsrcLkdZA9NoufARJJ7RhPwV5A7Zwb5y77AmACpffqRPfRU0vr0I0FSceb58O2oAwOulEg8feJw9oggr3oHq7evZ29ZCQDJyclkZWWRktKC2/0ZjU2L8Pv3AUJMzGBc0SP4LDCWOfWZrKr3YBAiHcKpMZEMiYmkb6SXXpEe0jxukt0u4lxOYl0OPA7tsVTd71gDfSLwG2PMBcHHDwAYYx5qs8+c4D5fiIgL2AOkmsO8+NcN9N8++Qfe7T/yoHaDdHj/UPsc/vltdaa9zX051D5Heq5pV9/Rve/Xr9va6sdFs3g5lITmak6p28rw2i2MqsrjlOptePw+CHZRVDng1XgHH8U4aAie2TqMISpgfVAT2WT4Rq7h9PWGzDbX+QQEml1QGZ9NSY+JVMcNo8Wb8tXvy1dBoGkjfv8uAv4SwN+6LdIZS1b0INIj+5LkTSfCGdW6rVLq2OUoo9hRQZnU0ODYvxBGgLi4UhIS9xAXW0p0TCUej3UBUw2xbGQYmxjKdvpTRDb1Et3hMXEYP26aceHHiQ+n8eMgcJjfsiM5+u7Pr/te8jXe61icDO/XnTVOKNnE49/60dd67uECvTNrivYEdrV5XAiMP9Q+xhifiOwDkoGydoXcCdwJ0KtXr04V315UADKaSg+8ZtvXP8QPRIL/XzEc/A/g4EiUjl/LHPwe0sH2w9ZxqFo7eO7BtZrWPVqfa9ruI195jhHrvkMcOBxCAIPPWCFuRAiIEBAXAePBh5eA8eIyLlx+Nw6/F6ffgycAHgOeALgCUfRqimRAgxAdSKY4KhKX+epKQZc2w7SKAAWuJopcTVQ5fDRIoPWbLxsFs0ZBVF2AlBIfsdUBIusDuJsNLl8JXt9MelS+i1/i8DnT8TsT8TtiMI5IAuIBRzZGsvDTSMA0EaCZgPFRWFvErtpdgMHriCTaHUuEMxqvMwK3w0M/cZEjKfgchganjyZHgOaaRHy7c6iVANViEFcTrsganN46+nkayXHnIc514PLR4HJR6Y5jnyuaOlckDRJBg8NLs3hoETc+cRLAgT94+3UcLkpEvk7UH+69Ovu/AdPuUfd2QR3x/bqwnK/3vbU/Ceu82Mbjs8pWZwK9o4rb/351Zh+MMU8CT4J1ht6J9z7I/Xf/gvu/zhNVFxgLXHfEvSYd/0KUUh3ozOlEIZDd5nEWsPtQ+wS7XOKBiq4oUCmlVOd0JtCXAzki0ldEPMC1wKx2+8wCbg7evwqYf7j+c6WUUl3viF0uwT7xe4A5WMMWnzXGbBCR3wErjDGzgGeA/4hIPtaZ+bXHs2illFIH60wfOsaY2cDsdm2/bnO/Ebi6a0tTSil1NHQgrVJK2YQGulJK2YQGulJK2YQGulJK2UTIZlsUkVJg59d8egrtrkI9QWhdR+9ErU3rOjpa19E5lrp6G2NSO9oQskA/FiKy4lBzGYSS1nX0TtTatK6jo3UdneNVl3a5KKWUTWigK6WUTZysgf5kqAs4BK3r6J2otWldR0frOjrHpa6Tsg9dKaXUwU7WM3SllFLtaKArpZRNnHSBLiJTRWSziOSLyM9DXc9+IrJDRNaJSK6IhGyxVBF5VkRKRGR9m7YkEZkrIluDt4knSF2/EZGi4DHLFZELQ1BXtogsEJE8EdkgIj8Itof0mB2mrpAeMxGJEJFlIrImWNdvg+19RWRp8Hi9Hpxq+0So63kR2d7meB28fmX31OcUkdUi8n7w8fE5XsaYk+YLa/rebUA/wAOsAYaGuq5gbTuAlBOgjsnAaGB9m7aHgZ8H7/8c+NMJUtdvgJ+E+HhlAKOD92OxFkQfGupjdpi6QnrMsFYniwnedwNLgQnAG8C1wfYngO+cIHU9D1wVyt+xYE0/Al4B3g8+Pi7H62Q7Qx8H5BtjCowxzcBrwGUhrumEYoz5lINXi7oMeCF4/wXg8m4tikPWFXLGmGJjzKrg/RogD2uN3JAes8PUFVLGUht86A5+GeAc4M1geyiO16HqCjkRyQIuAp4OPhaO0/E62QK9owWrQ/5LHmSAj0RkZXAx7BNJD2NMMVhBAaSFuJ627hGRtcEumW7vCmpLRPoAo7DO7k6YY9auLgjxMQt2H+QCJcBcrL+aq4wxvuAuIfl32b4uY8z+4/X74PH6m4h4u7su4BHgp0BwxXSSOU7H62QL9E4tRh0ik4wxo4FpwPdEZHKoCzoJPA70B0YCxcBfQlWIiMQAbwE/NMZUh6qO9jqoK+THzBjjN8aMxFpfeBwwpKPdureqg+sSkVOAB4DBwGlAEvCz7qxJRC4GSowxK9s2d7Brlxyvky3QO7NgdUgYY3YHb0uAd7B+0U8Ue0UkAyB4WxLiegAwxuwN/iMMAE8RomMmIm6s0HzZGPN2sDnkx6yjuk6UYxaspQr4BKuvOiG4QDyE+N9lm7qmBruujDGmCXiO7j9ek4BLRWQHVhfxOVhn7MfleJ1sgd6ZBau7nYhEi0js/vvA+cD6wz+rW7VdxPtmYGYIa2m1PzCDriAExyzYn/kMkGeM+WubTSE9ZoeqK9THTERSRSQheD8SOBerf38B1gLxEJrj1VFdm9r8T1mw+qm79XgZYx4wxmQZY/pg5dV8Y8z1HK/jFepPf7/Gp8UXYn3ivw34n1DXE6ypH9aImzXAhlDWBbyK9ad4C9ZfNLdj9dnNA7YGb5NOkLr+A6wD1mIFaEYI6joD68/dtUBu8OvCUB+zw9QV0mMGDAdWB99/PfDrYHs/YBmQD8wAvCdIXfODx2s98BLBkTCh+ALO4sAol+NyvPTSf6WUsomTrctFKaXUIWigK6WUTWigK6WUTWigK6WUTWigK6WUTWigK6WUTWigK6WUTfx/vPZWscRUO7wAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "u0 = uini.copy()\n", "\n", "# total simulation time\n", "TotTime = 7\n", "# Times of outputs\n", "nsteps = 10\n", "\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "\n", "for m in range(int(TotTime/dt)):\n", " u0, u = do_timestep(u0, u)\n", " if m in range(0, int(TotTime/dt), int(TotTime/(dt * nsteps))):\n", " ax.plot(x[1:-1], u.copy()[1:-1])\n", " \n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "SQIboYWF4qNq", "metadata": { "id": "SQIboYWF4qNq" }, "source": [ "### In the co-moving frame" ] }, { "cell_type": "code", "execution_count": null, "id": "oXvfytkM4qNq", "metadata": { "id": "oXvfytkM4qNq", "outputId": "4dcc2db1-d07a-4f69-ba4f-00cbea2cb07c" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAfDUlEQVR4nO3deXycZbn/8c+VyZ62SZe0TZd0o9BSllJKAQsIIgiIgAjKooAHX/2pcABxOYCCeFxYFBTFH5wqHBaRTRaLVlkKWpA1XehK6UqbbkmbNF2yZ67zx0zbGNImaSeZ5p7v+/WazrPNzDV3J988ued+nsfcHRER6f7Skl2AiIgkhgJdRCQQCnQRkUAo0EVEAqFAFxEJRHqyXrhfv34+fPjwZL28iEi3NGvWrE3uXtjauqQF+vDhwykpKUnWy4uIdEtm9tGe1qnLRUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEG0GupkNNbPXzGyxmS00s2tb2eZkM6sys7nx2y2dU66IiOxJe4YtNgLfdvfZZtYTmGVmL7v7ohbbve7uZye+RBERaY82A93d1wPr49PbzGwxMBhoGehygJm7ZguvLt4YmzHjnCOLOKh/z+QWJSKdpkMHFpnZcOAo4J1WVh9vZu8D64DvuPvCVh4/BZgCUFxc3NFapYN+M2MpMz4owwzcoXxbHbedf3iyyxKRTtLuL0XNrAfwDHCdu29tsXo2MMzdjwR+Azzf2nO4+1R3n+juEwsLWz1yVRKoyZ0jh+Sz8rbPUpSfTTSqi5mIhKxdgW5mGcTC/DF3f7blenff6u7b49PTgQwz65fQSkVEZK/aM8rFgAeAxe5+9x62GRjfDjObFH/ezYksVDqu5dUFHe2hi4SsPX3ok4GvAPPNbG582U1AMYC73w9cAHzDzBqBGuAi18VKRUS6VHtGubwBWBvb3Avcm6iiJIFifzjt/T9QRIKgI0VTiP5mEgmbAj1gym+R1KJAD9zOrpb4d9YiEjAFuohIIBToAWs50EhdMCJhU6CLiARCgR44dZ2LpA4FegrRsEWRsCnQRUQCoUAP3O5hi0ktQ0S6gAI9hejkXCJhU6AHTH3mIqlFgR64nUeIqstFJHwKdBGRQCjQA/axPnN1wYgETYEuIhIIBXrgdg1b1CUuRIKnQE8h6nERCZsCPWAatiiSWhTogds5XFHDFkXCp0BPIS3Pjy4iYVGgB0z5LZJaFOgiIoFQoAdu53BFdaGLhE+BHrCWR4qqB0YkbAp0EZFAKNBDt2vYojpdREKnQE8hGvUiEjYFesAU4CKpRYEeOGtxLyLhUqCnEO2wi4StzUA3s6Fm9pqZLTazhWZ2bSvbmJn92syWmdk8M5vQOeVKRyjARVJLeju2aQS+7e6zzawnMMvMXnb3Rc22ORMYHb8dC9wXvxcRkS7S5h66u69399nx6W3AYmBwi83OBR7xmLeBAjMrSni10mGmTnSRlNGhPnQzGw4cBbzTYtVgYE2z+VI+HvqY2RQzKzGzkvLy8o5VKh3X8pKiGvYiErR2B7qZ9QCeAa5z960tV7fykI+lh7tPdfeJ7j6xsLCwY5WKiMhetSvQzSyDWJg/5u7PtrJJKTC02fwQYN3+lyf7SyfnEkkd7RnlYsADwGJ3v3sPm00DLouPdjkOqHL39QmsUxJAHS4iYWvPKJfJwFeA+WY2N77sJqAYwN3vB6YDZwHLgGrgq4kvVToqdrZF7ZuLpIo2A93d36CNVPDYt21XJaooERHpOB0pGjhrfrZF9bmIBE2BHjCNUhRJLQp0EZFAKNADt6vLJblliEgXUKCnkJbXGBWRsCjQA6b4FkktCvTA7TpSVH0uIsFToKcQjXoRCZsCPWA6u6JIalGgi4gEQoEeuN3DFk1dLiKBU6AHTPktkloU6CIigVCgpwgNWxQJnwI9hehIUZGwKdADpi9BRVKLAl1EJBAK9MBZs85z7bGLhE2BHjDlt0hqUaCLiARCgR64nR0uZqY9dpHAKdBFRAKhQA+ZvgUVSSkK9MDpmqIiqUOBnkK0wy4SNgV6wJTfIqlFgS4iEggFeuB2D1sE7bOLhE2BHjD1mYukFgW6iEgg2gx0M3vQzMrMbMEe1p9sZlVmNjd+uyXxZcq+2nlyLl3gQiR86e3Y5iHgXuCRvWzzurufnZCKpNOoC0YkbG3uobv7TKCiC2qRBNMVikRSS6L60I83s/fN7G9mNm5PG5nZFDMrMbOS8vLyBL20iIhAYgJ9NjDM3Y8EfgM8v6cN3X2qu09094mFhYUJeGlpy65hi+hsiyKh2+9Ad/et7r49Pj0dyDCzfvtdmew39ZmLpJb9DnQzG2jxoRRmNin+nJv393lFRKRj2hzlYmaPAycD/cysFPghkAHg7vcDFwDfMLNGoAa4yF37hgeKXWdbNNB/i0jY2gx0d7+4jfX3EhvWKCIiSaQjRQOmHXKR1KJAD541+1dEQqZATyHaYRcJmwI9YApwkdSiQBcRCYQCPXDW7AoX+pJUJGwK9IBp3LlIalGgi4gEQoEeOGt2r/11kbAp0EVEAqFAFxEJhAJdRCQQCvTA6WyLIqlDgR4w5bdIalGgi4gEQoEeONPZFkVShgJdRCQQCvSAuQ4lEkkpCvTA7R7lopNziYROgS4iEggFesC0Ry6SWhToIiKBUKAHblcfOvqSVCR0CvSAKb5FUosCXUQkEAr0wO06UtT0JalI6BToIiKBUKAHTKfLFUktCvQUonwXCZsCPXS7hi3qfIsioVOgB0w75CKppc1AN7MHzazMzBbsYb2Z2a/NbJmZzTOzCYkvU0RE2tKePfSHgDP2sv5MYHT8NgW4b//LkkSxZhM6UlQkbG0GurvPBCr2ssm5wCMe8zZQYGZFiSpQRETaJxF96IOBNc3mS+PLPsbMpphZiZmVlJeXJ+ClZa+0Qy6SUhIR6K0Nn2g1Stx9qrtPdPeJhYWFCXhpaYvZ7muKatiiSNgSEeilwNBm80OAdQl4XhER6YBEBPo04LL4aJfjgCp3X5+A55X9pB1ykdSS3tYGZvY4cDLQz8xKgR8CGQDufj8wHTgLWAZUA1/trGJl/yjgRcLWZqC7+8VtrHfgqoRVJAm18wsOnW1RJHw6UjRgOjmXSGpRoIuIBEKBHjhrfnIu7bCLBE2BLiISCAV6wLRDLpJaFOgpRCfnEgmbAj1wzYctikjYFOgB06hFkdSiQE8hCniRsCnQA7frbIvqchEJngJdRCQQCvSAaVSLSGpRoAdu1ygXTPEuEjgFuohIIBToAdOoFpHUokBPITqdrkjYFOih23m2RQ1bFAmeAj1g2iEXSS0K9BSifBcJmwI9cIb6WkRShQJdRCQQCnQRkUAo0AO365qiZvqSVCRwCnQRkUAo0AOmA4lEUosCPYUo3kXCpkAPnLW4F5FwKdBFRAKhQA/Yx7pY1KcuEjQFeuBMJ+cSSRntCnQzO8PMlpjZMjO7oZX1V5hZuZnNjd++lvhSRURkb9Lb2sDMIsBvgdOAUuA9M5vm7otabPqku1/dCTXKPnKP/1Nbpe4WkRTQnj30ScAyd1/h7vXAE8C5nVuWJEKhb+La5V+D24v5Xtn3yItuS3ZJItKJ2hPog4E1zeZL48ta+oKZzTOzP5nZ0NaeyMymmFmJmZWUl5fvQ7nSEXc3/pTBtUth/KUcWjuX72y7M9kliUgnak+gt/Z1Wsu/318Ahrv7EcArwMOtPZG7T3X3ie4+sbCwsGOVSscsfYXR/hElBWfAef+f+dkTmdAwGzbMT3ZlItJJ2hPopUDzPe4hwLrmG7j7Znevi8/+Djg6MeXJPnv9LgD+MvA/AXi0z3/+23IRCU97Av09YLSZjTCzTOAiYFrzDcysqNnsOcDixJUoHVZTCavf5F82gdqMXgBsyhzM0shBsPA5aGpIcoEi0hnaDHR3bwSuBl4kFtRPuftCM/tvMzsnvtk1ZrbQzN4HrgGu6KyCpR1mPQTA85FP/9viV7NOjU0s+nMXFyQiXaFd49Ddfbq7H+zuo9z9p/Flt7j7tPj0je4+zt2PdPdT3P2Dzixa2rAo9gfUP5n0b4tfyj4tNjH/6a6uSES6gI4UDU1DLaybDQd9Gsx2HykK1Fk2DDwClr2icekiAVKgh2bFa7H7see0vn7MZyHaqNEuIgFSoIdm4fOx+0PPbf385+POj2/3bFdVJCJdRIEemmWvQH4x5BR8bJU7UHgwRLLgwxe7vjYR6VQK9JDs2AzVm2DESQA4jWxoKOGpJU9RZfNxorHtRpwIZYs0fFEkMG2enEu6kSXTY/djzmJp5VLqB9xJSW05JW/HFmf0GMza7QcxePRnYnvyK2fCQacmr14RSSjtoYfkozcBWNFvOF+e/mU8so2js69jxoUzGMlXaUjbwCV/vYSy4cfHtl/5zyQWKyKJpkAPyao3qB4wjitfvZqGaAMZZVczKHMS/XP7U8hJFG6/horaCq588yYa8wp3/QIQkTAo0ENRUwlVq7mpIJtNNZu446Q7SGv495NiZjeN5sZJN7Jq6yruKBoCpe+pH10kIAr0UHz4IrOysphRX86E/hM4bdhpmDfRq2wtdcuW0WNrBe7OxWMuZlivYTzRtJkVGemw5p1kVy4iCaIvRQPhS1/he/37AnD7J35GxSOPcvdz99GvupIVv4T/B5Tn96dq6LXcdfIvuOAvF3JjYV+e/PDvMPyE5BYvIgmhQA/EtPL3KMtN54r+51D71WvYsmgRm/qPYumZX+KCyaN59G9zOWjuTNZ///vkHnssZ114ItM3v87bpa9zXLKLF5GEUKAHwGuq+GVmAwMq0vjc72ZSW1FB0W238eV5uZw5ZgBMOpj5Gwr566jJ/KH3KjbedjtXrhvMv85z7szcwLPRKKSp902ku9NPcQD+9M6dRGvS+PnjEZqqqii45wHm1R3KlzZGKHptM4/+4C0Of3cbJy9vYnGPyRT8/F6a1qzl7scirI1m8K95rV5gSkS6GQV6N+fu/M/yv3Lz402k70hj3ZX38tyz1Sx6Yx3lmbD1kDw+ddlY1g3LYnum8f6MNfz5xQgbp/yG3C1Rbnyqibvm/y7Zb0NEEkBdLt3c9JXTuXRaPb3qh1Fy2nfY8UEDxWOqychcSek7s0nfUscby3PJq26ilnRGTZpA9fZhLFiSTf7pP2fszNs58W9llJxWwsSBE5P9dkRkP5gn6bzYEydO9JKSkqS8dkh++K1jOKXkYBaOvRxLW0562iy2bd5Az5y+1OWMZkBBf4bn57KkrIJttZWkbf+QHXVV5PcfSl3dBIwhHDF/Kq+ctY1bb34l2W9HRNpgZrPcvdW9L+2hd2NvvPYoJ8w7mvmHnIzXPkl2tIFxg05gaNEYbAfUUE8tDViVMTGzPxmZ6WT1OoWmPOejqgUsqniN+rRezDnsy5ww4zkWfv41xh12SrLflojsIwV6NxWtqaHsrndYM6SYnnUzOKLfKWRm9+WjSDmL0uezqUcFRLaSnl6HpUWJNkVobMwi4r3pZ/kUZedzYvFlbK9ey4LN01k4fDxbv/8Yhz5zIpauj4VId6Sf3G7qhW/9grI+mRzVsy91+SN4J2c59HmJgoIN5PeupykzhyrLZwe5pNNIBg1k0UAv30F+XQ0VldmsqBxIesUYDs47nbSqchZFe/PKTfdw2p3fTvbbE5F9oEDvhl697WHwCIOKD2XxwDfxQZtZnT+URXY8y+0QKuiz5wcbkAOFOeUcVLSEsb4Qq3wT1g1k2Ibj2Fa+lnemTuPYKXu4hJ2IHLAU6N3Mq3c9wdYd5ZROXsIHRVW8Fzmbj2wkAIMz4ZMFBRzeM5eDcrO4/uFZnDq6kKtOHc0tLyyktKqWG88fx4rqOuZt68Xblf14q/EE6Auj+nzIpDHvMLq0Aluax9sPTee4K85K8rsVkY5QoHcT0aYoL/z8HuaMXMPb48cxz64hahGOrPmQKw4r4vR++YzKzf63x6TvaCTfjYPzsilogMrqJj5buPvSdO7Okupa/l5WwXOLl/N4zleIDG/g6GElTNr4ClvuXsNnvjUFM+vqtysi+0CB3g0sWzyHexf9lX9OOor1dgo9mrYyOW0Bd7x7ByOPPA+Kv9iu5/EWl402M8bk5TBmxGCuK3mBxQvmcv3E65jPEbxbdDzFA1fx6tM/4dpjL2bwsIM6462JSALpSNEDWGnVBr4x/X7OXL+NJ/qcTUa9c/6/fk+fjTfweOVbjKxZC8d/c6/PsXPfus2d7OOvZmz1Sp5pXE5e6TVc+NbjNDTm8Ejh5/j08o+4/u9T2VxdmZD3JSKdQ4F+AFq5rYKv/fNpPjFrFc/lHEdx3Tq+/sY/+Ondd/LmoH9w2yduIX3e4zDgMMgfkpgXHXw05BWS+9a93DTp27zd76/cc9sd/Md7MylsqOSPWZM45u1FXPPGs6yv2Z6Y1xSRhFKXywFkduUm7lz0LjPrBgAjOaZuNp+aU0avBbkcteZ5vnlJBccMOZ5T130Ye8AxV+71+VoeBLzXg4LNYMJl8PpdXBLN4ekBo7n1vKXc9fizjJh/PlVHfciMcYN4miN47q1FnJZbwXcPPZ6xvfL36z2LSOIo0JOsIepM27CG+5YvZUFjX7K9NyfXzuS0hTU0Lqynx46eHLHmab53SR3R/B7cfdIv4I4RsQdPuCKxxZzwLXj9LvjbDdz3jX9y5jNnctdF9fzXo39kTvQyLlywjPOOXMvLh/TiZSYxfdZKJmZW8M1RY/nMgIFE9OWpSFIp0JNk2Y5aHl69hD9t3EGl51LojZy343lO/CCDrKX1rGscyKiNKxhV9RI/uyKHtbnVPHHGY/R8+z7wJjj1h+06h/nOESrtytqsnjDxSih5gIFLXub+0+7nay99jQeuHMaUBx5iSfHZrJmbx6WLK/jcIS/y+sHGGz6Z/1i8kf5LVnBxUQFfLh7N0OzM/WwdEdkXCvQutKa2nmnr1/KndetYXJ+HeRNH+mLOr1zAhCUDyS7NZsWOKnLyjuHo+Y9SWFDDT6bkM8dK+dUnf8WY9F6xPeic3jD52jZfr+WJ19p1GrbTfwIlD8LLN3PstfO4+bib+fHbP8avGs1VD8+gT/kglhx+AY1z3uYLywby6eI3mTO6kvcKxnPP2sO5Z+0iJmTv4PzBQ/ncgCIGZGXsW2OJSIcp0DtRdVOU97ZsZ0bZal7dXMWyhjwAhvlGzq+fzdjSTYxYPZ5IeRGrqtfTlDmR4Y3pDJ/5Y+yUo7nu+BWsbdrMHSfcwalFx8HvPgVN9fDFZyEt0jlFZ+bChf8LT18Bf/wSX7z8BRqiDdz+7u2snzKUH81Ip+CN21kx4assrc0la/FcTl1dzOTCKlYMv4/FRQN4zyfxg+V5/GB5BYdm7eDUfn04pXAIR+fnkaUrI4l0mnYFupmdAdwDRIDfu/vtLdZnAY8ARwObgS+5+6rElnpg29rYxLIdtcyvKmfOljLmbq1laUMeTUSIeCOjWcX59YsYVb6RIaUj6Lt5JNu2GCvrl0HmEQwin+JZD9CjD7x21QTuy3uH/hn9eeDUB5iYVQgPngHlH8C5v4URJ3a4PqMD/dvjPg8VK2HGj+APn+fSLzzIoLxB3PzmzVzyyXVcd/jRHPvsHxlUm8eqwy/ig1onbfl8Rq4tYlzvYZzct4Q1Q//M8r4DWeiHcW9pDr9Zu5wMGjkks5qjeuUyvvcADs/vy6jcLPIinfTLSSTFtBnoZhYBfgucBpQC75nZNHdf1GyzK4FKdz/IzC4C7gC+1BkFdyZ3pzoapbopyvbGBrbWV7O1fgdVdTvYUruDqrpqqupq2VZfz+aGKJua0iiLZlJmeWy1HrueJ9sbGcFqTm9cxZBtGxlW3kifstHkVQ2jblsWG+oq2WBZ5KSNZUT5+wxaP5WG/s6Mz+bwxIgNRDKquGzImXy99xH0fPcRmPdU7InPuw/GX9I1jXHi9RDJhFd+CL8+ilOOupTnj/gO95a/ya/sZezyBr6yPIcT3/ofRm7LZd2gyazvfzC1G1eTt3Eb41YM48ieAzm9YDWb+r/KR4VZlPYYyEofyZPlI3h002Ziv/uhwLcx0Krpnx5lYEYahdlZ9M/OoTAnjz5ZefTOyqNHRjY90rPJS4+QlZZGuqEjWEVaaPMCF2Z2PHCru38mPn8jgLvf1mybF+PbvGVm6cAGoND38uT7eoGLH039Gc+PGo9ju25As3lg57Ttno9i8eUt53ff6sjCrX1dAlleQx8q6O2b6dtUSe/6LRTW7KBoaz39K7PJqOxPZHsa0doatjbUs6UxQmZTNgU1jfTeuobs6iWs77WOxcVGyWhjXT9nXH0TJ+2o5wvbK+kbjQJQQxYzM07g0ayL2Jg2oENttax8O5cfP5xbzxnHVY/N5qVFGxjeN69DzzGsaTWX1j3J5Ma3yKQRgNJIJs/27MXruVl8mGGM2GAcvSzK2NXQf3sxO3oeTGWvwVRlR2iM1NA74vTMzMRycmjq0UB9QTkb+tSxoWcmm3J6UJFZwOZIXyrpQwV9aLCsNusyj5JOI0aUNJw0ovFbE2m+e952fXOw838+/vhm3yjsaRsRAOuEiwAdW76E+y6+fp8eu78XuBgMrGk2Xwocu6dt3L3RzKqAvsCmFoVMAaYAFBcXt6v4lnKjUFRXTtrOHz5vHuu75yEe106zKHfYOe+7fxWkeWx9pjeQHa0js6mBrKYmMhqjZNVHyaqHrAbIrTMya5yMGidSA9F6wxvAm6I0RaM0ueHRRraxiTRfS5pXE82opS67lqaCWiryYG1uhLrsCFlZmeT7QQyNZnK4Z5JflUGENEiDf/UuYENGMRszh1KeMZjGtEx6Ab062FYHD+jJ544cBMCFE4d87ND/9jmUp/kR06LVDK9dzMD61eQ3buZwr+Wo7fVEvZaNOTsoH9/IhxMamcN2Mra/R862t8jaHiWjNo3q+t7U1/SC7XlEyzKJeoQelsGhEUiPNGLplVh6FZ69FHLqqctrZHtOhNqcCNWZEWozItSlR2iIRKhPS6cxLY1Gi9CUlkbUjKjFPgFRi8V71OJRbs1/ycc+EbH/8WbT8Q+OK8aDkPDo7aSPRc/ahk553vYEemtvqWW7tWcb3H0qMBVie+jteO2P+e7Xb+K7+/LAFHfyIf05+ZD++/ksHe+7F5Gu057+hVJgaLP5IcC6PW0T73LJByoSUaCIiLRPewL9PWC0mY0ws0zgImBai22mAZfHpy8AXt1b/7mIiCRem10u8T7xq4EXiQ1bfNDdF5rZfwMl7j4NeAB41MyWEdszv6gzixYRkY9r1zh0d58OTG+x7JZm07XAhYktTUREOkKH7YmIBEKBLiISCAW6iEggFOgiIoFo89D/Tnths3Lgo318eD9aHIV6gFBdHXeg1qa6OkZ1dcz+1DXM3QtbW5G0QN8fZlayp3MZJJPq6rgDtTbV1TGqq2M6qy51uYiIBEKBLiISiO4a6FOTXcAeqK6OO1BrU10do7o6plPq6pZ96CIi8nHddQ9dRERaUKCLiASi2wW6mZ1hZkvMbJmZ3ZDsenYys1VmNt/M5ppZx6+tl7g6HjSzMjNb0GxZHzN72cyWxu97HyB13Wpma+NtNtfMzkpCXUPN7DUzW2xmC83s2vjypLbZXupKapuZWbaZvWtm78fr+lF8+QgzeyfeXk/GT7V9INT1kJmtbNZe47uyrmb1Rcxsjpn9JT7fOe3l7t3mRuz0vcuBkUAm8D5waLLrite2Cuh3ANRxEjABWNBs2Z3ADfHpG4A7DpC6bgW+k+T2KgImxKd7Ah8Chya7zfZSV1LbjNjVyXrEpzOAd4DjgKeAi+LL7we+cYDU9RBwQTI/Y/Gargf+CPwlPt8p7dXd9tAnAcvcfYW71wNPAOcmuaYDirvP5ONXizoXeDg+/TBwXpcWxR7rSjp3X+/us+PT24DFxK6Rm9Q220tdSeUx2+OzGfGbA58C/hRfnoz22lNdSWdmQ4DPAr+Pzxud1F7dLdBbu2B10j/kcQ68ZGaz4hfDPpAMcPf1EAsKYH8vLppIV5vZvHiXTJd3BTVnZsOBo4jt3R0wbdaiLkhym8W7D+YCZcDLxP5q3uLujfFNkvJz2bIud9/ZXj+Nt9cvzSyrq+sCfgV8D4jG5/vSSe3V3QK9XRejTpLJ7j4BOBO4ysxOSnZB3cB9wChgPLAeuCtZhZhZD+AZ4Dp335qsOlpqpa6kt5m7N7n7eGLXF54EjG1ts66t6uN1mdlhwI3AGOAYoA/wX11Zk5mdDZS5+6zmi1vZNCHt1d0CvT0XrE4Kd18Xvy8DniP2QT9QbDSzIoD4fVmS6wHA3TfGfwijwO9IUpuZWQax0HzM3Z+NL056m7VW14HSZvFatgD/INZXXRC/QDwk+eeyWV1nxLuu3N3rgP+l69trMnCOma0i1kX8KWJ77J3SXt0t0NtzweouZ2Z5ZtZz5zRwOrBg74/qUs0v4n058Ock1rLLzsCM+zxJaLN4f+YDwGJ3v7vZqqS22Z7qSnabmVmhmRXEp3OATxPr33+N2AXiITnt1VpdHzT7pWzE+qm7tL3c/UZ3H+Luw4nl1avufimd1V7J/vZ3H74tPovYN/7Lge8nu554TSOJjbh5H1iYzLqAx4n9Kd5A7C+aK4n12c0Alsbv+xwgdT0KzAfmEQvQoiTUdQKxP3fnAXPjt7OS3WZ7qSupbQYcAcyJv/4C4Jb48pHAu8Ay4Gkg6wCp69V4ey0A/kB8JEwybsDJ7B7l0intpUP/RUQC0d26XEREZA8U6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gE4v8AwLOdLLei2ZcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "u0 = uini.copy()\n", "\n", "# total simulation time\n", "TotTime = 15\n", "\n", "# Times of outputs\n", "nsteps = 10\n", "\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "\n", "for m in range(int(TotTime/dt)):\n", " u0, u = do_timestep(u0, u)\n", " if m in range(0, int(TotTime/dt), int(TotTime/(dt * nsteps))):\n", " indexshift=int(v*m*dt/dx)\n", " ushifted=np.pad(u.copy()[indexshift:], [(0, indexshift)], mode='constant')\n", " ax.plot(x[1:-1], ushifted[1:-1])\n", " \n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "id": "zCMOSJSx4qNr", "metadata": { "id": "zCMOSJSx4qNr" }, "source": [ "### Diffusion equation without bias" ] }, { "cell_type": "markdown", "id": "QUYpWB2T4qNr", "metadata": { "id": "QUYpWB2T4qNr" }, "source": [ "Only works far away from the boundaries." ] }, { "cell_type": "code", "execution_count": null, "id": "Ay1_obuO4qNr", "metadata": { "id": "Ay1_obuO4qNr" }, "outputs": [], "source": [ "def do_timestep(u0, u):\n", " # Propagate with forward-difference in time, central-difference in space\n", " u[1:-1] = u0[1:-1] + D * dt * (\n", " (u0[2:] - 2*u0[1:-1] + u0[:-2])/dx2 ) \n", "\n", " u0 = u.copy()\n", " return u0, u\n" ] }, { "cell_type": "code", "execution_count": null, "id": "tyqOg9Le4qNr", "metadata": { "id": "tyqOg9Le4qNr", "outputId": "aad543ed-62c7-401c-87cd-f561e1ac2eb7" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAe40lEQVR4nO3deXxU9b3/8dcnCyFA2EIIu6yiIAKaKly8ilbrWu2iv+qtre21l1urta1drl1+Lr2/Vm9va621V+Vqa7V1qbZWbLUVxbqLBgg7SNgkbAlLFkhCkpnP748ZMIZAEplhku+8n49HzJxzvnPmM9+Qtyff+Z5zzN0REZGuLyPVBYiISGIo0EVEAqFAFxEJhAJdRCQQCnQRkUBkpeqFBwwY4CNHjkzVy4uIdEkLFizY4e4FrW1LWaCPHDmS4uLiVL28iEiXZGYbD7VNQy4iIoFQoIuIBEKBLiISCAW6iEggFOgiIoFoM9DNbLiZvWRmK81suZl9rZU2M82sysxK4l83JadcERE5lPZMW2wCvunuC80sD1hgZnPdfUWLdq+6+0WJL1FERNqjzUB3963A1vjjGjNbCQwFWga6dDIlmyqZt3J7bMGMiycPZuzAvNQWJSJJ06ETi8xsJDAVmN/K5ulmthjYAnzL3Ze38vxZwCyAESNGdLRW6aBfvriGF1eVYwbuUFGzj9s+NSnVZYlIkrT7Q1Ez6wX8Efi6u1e32LwQOMbdJwO/BP7c2j7cfba7F7l7UUFBq2euSgJF3Jk8rA/rb7uQwX26E43qZiYiIWtXoJtZNrEw/727/6nldnevdvc98cfPAtlmNiChlYqIyGG1Z5aLAQ8AK939jkO0GRRvh5mdEt/vzkQWKh3X8u6Cjo7QRULWnjH0GcDngKVmVhJf9z1gBIC73wtcClxjZk1AHXC562alIiJHVXtmubwGWBtt7gbuTlRRkkCxP5wO/wMUkSDoTNE0or+ZRMKmQA+Y8lskvSjQA7d/qCX+mbWIBEyBLiISCAV6wFpONNIQjEjYFOgiIoFQoAdOQ+ci6UOBnkY0bVEkbAp0EZFAKNAD9/60xZSWISJHgQI9jejiXCJhU6AHTGPmIulFgR64/WeIashFJHwKdBGRQCjQA3bQmLmGYESCpkAXEQmEAj1wB6Yt6hYXIsFToKcRjbiIhE2BHjBNWxRJLwr0wO2frqhpiyLhU6CnkZbXRxeRsCjQA6b8FkkvCnQRkUAo0AO3f7qihtBFwqdAD1jLM0U1AiMSNgW6iEggFOihOzBtUYMuIqFToKcRzXoRCZsCPWAKcJH0okAPnLX4LiLhUqCnER2wi4StzUA3s+Fm9pKZrTSz5Wb2tVbamJndZWalZrbEzE5KTrnSEQpwkfSS1Y42TcA33X2hmeUBC8xsrruvaNbmfGBc/OtU4J74dxEROUraPEJ3963uvjD+uAZYCQxt0ewS4CGPeQvoa2aDE16tdJhpEF0kbXRoDN3MRgJTgfktNg0FNjVbLuPg0MfMZplZsZkVV1RUdKxS6biWtxTVtBeRoLU70M2sF/BH4OvuXt1ycytPOSg93H22uxe5e1FBQUHHKhURkcNqV6CbWTaxMP+9u/+plSZlwPBmy8OALUdenhwpXZxLJH20Z5aLAQ8AK939jkM0mwN8Pj7bZRpQ5e5bE1inJIAGXETC1p5ZLjOAzwFLzawkvu57wAgAd78XeBa4ACgFaoEvJr5U6ajY1RZ1bC6SLtoMdHd/jTZSwWOftl2bqKJERKTjdKZo4Kz51RY15iISNAV6wDRLUSS9KNBFRAKhQA/cgSGX1JYhIkeBAj2NtLzHqIiERYEeMMW3SHpRoAfuwJmiGnMRCZ4CPY1o1otI2BToAdPVFUXSiwJdRCQQCvTAvT9t0TTkIhI4BXrAlN8i6UWBLiISCAV6mtC0RZHwKdDTiM4UFQmbAj1g+hBUJL0o0EVEAqFAD5w1GzzXEbtI2BToAVN+i6QXBbqISCAU6IHbP+BiZjpiFwmcAl1EJBAK9JDpU1CRtKJAD5zuKSqSPhToaUQH7CJhU6AHTPktkl4U6CIigVCgB+79aYugY3aRsCnQA6Yxc5H0okAXEQlEm4FuZr82s3IzW3aI7TPNrMrMSuJfNyW+TPmw9l+cSze4EAlfVjvaPAjcDTx0mDavuvtFCalIkkZDMCJha/MI3d1fAXYdhVokwXSHIpH0kqgx9OlmttjMnjOziYdqZGazzKzYzIorKioS9NIiIgKJCfSFwDHuPhn4JfDnQzV099nuXuTuRQUFBQl4aWnLgWmL6GqLIqE74kB392p33xN//CyQbWYDjrgyOWIaMxdJL0cc6GY2yOJTKczslPg+dx7pfkVEpGPanOViZo8CM4EBZlYG3AxkA7j7vcClwDVm1gTUAZe769iwszhwtUUD/VhEwtZmoLv7FW1sv5vYtEYREUkhnSkaMB2Qi6QXBXrwrNl/RSRkCvQ0ogN2kbAp0AOmABdJLwp0EZFAKNADZ83ucKEPSUXCpkAPmOadi6QXBbqISCAU6IGzZt91vC4SNgW6iEggFOgiIoFQoIuIBEKBHjhdbVEkfSjQA6b8FkkvCnQRkUAo0ANnutqiSNpQoIuIBEKBHjDXqUQiaUWBHrj3Z7no4lwioVOgi4gEQoEeMB2Ri6QXBbqISCAU6IE7MIaOPiQVCZ0CPWCKb5H0okAXEQmEAj1wB84UNX1IKhI6BbqISCAU6AHT5XJF0osCPY0o30XCpkAP3YFpi7reokjoFOgB0wG5SHppM9DN7NdmVm5myw6x3czsLjMrNbMlZnZS4ssUEZG2tOcI/UHgvMNsPx8YF/+aBdxz5GVJolizBzpTVCRsbQa6u78C7DpMk0uAhzzmLaCvmQ1OVIEiItI+iRhDHwpsarZcFl93EDObZWbFZlZcUVGRgJeWw9IBuUhaSUSgtzZ9otUocffZ7l7k7kUFBQUJeGlpi9n79xTVtEWRsCUi0MuA4c2WhwFbErBfERHpgEQE+hzg8/HZLtOAKnffmoD9yhHSAblIeslqq4GZPQrMBAaYWRlwM5AN4O73As8CFwClQC3wxWQVK0dGAS8StjYD3d2vaGO7A9cmrCJJqP0fcOhqiyLh05miAdPFuUTSiwJdRCQQCvTAWfOLc+mAXSRoCnQRkUAo0AOmA3KR9KJATyO6OJdI2BTogWs+bVFEwqZAD5hmLYqkFwV6GlHAi4RNgR64A1db1JCLSPAU6CIigVCgB0yzWkTSiwI9cAdmuWCKd5HAKdDTQTSa6gpE5ChQoAfMHT6+/R74USHTa+eluhwRSTIFesD6eiVn7noMIg1cvftOMqMNqS5JRJJIgR6w85pejj2Y+jm6ez1TGhamtiARSSoFesAmRlfHHpzzQwCm7Xs9hdWISLIp0AM2JbqcDbkToEd/9mTkMa5xdapLEpEkUqCHqqGWflRT1n08ABuyxzIwsl3n/4sETIEeqvIVAGzLGQ3A4txTyKEBqjalsioRSSIFeqg2vgHA+h6TANiRWRhbv+PdVFUkIkmmQA9V7Q4AduQMB2B9TmzoZX/Qi0h4FOihWvcym62QSEY2ALszB8TW11WmsCgRSSYFeqia6mkk68CiWwbrM0fBhtdSWJSIJJMCPVQVqyi2Ez+wKmKZ0FiXooJEJNkU6CGq3gJAveU0u9oiLM6eAlXvaeqiSKAU6CGq2QrAYjv+A6v3WU7swe71R7siETkKFOghKisGYJf1/cDqVVnxgK/eerQrEpGjQIEeooY9AKy1EQfuJWoGuzP6xRbK3klRYSKSTO0KdDM7z8xWm1mpmd3YyvYvmFmFmZXEv76U+FKl3Ta+Cdk92Wc51EQ38fCKh2lgJ5szh8W2N9amtj4RSYqsthqYWSbwK+AcoAx4x8zmuPuKFk0fd/frklCjdFSkATKyaOq2mtfq7uW1d5ws682gzG9Aj3zYsijVFYpIErTnCP0UoNTd17l7A/AYcElyy5IjUrWJ6MjTqO9/Pxlkc+/Z9xKlkW25d8e21+1ObX0ikhTtCfShQPMrOpXF17X0aTNbYmZPmtnw1nZkZrPMrNjMiisqKj5EudKmaBR2reO5jHrIaGR09gXMGDqDfkylKWMnC4dPhdpdqa5SRJKgPYFuraxrOZH5GWCku58IvAD8trUduftsdy9y96KCgoKOVSrtEz/6fiQaC+0x2bE/pob5ZQD8jirYtRaikdTUJyJJ055ALwOaH3EPA7Y0b+DuO919X3zxf4GTE1OedNiudewxY0nDTjLqJpBp3QDIpjfdIiOY27CdCMCe7SktU0QSrz2B/g4wzsxGmVk34HJgTvMGZja42eLFwMrElSgdUl3G3J49AMiqK2o2bdHo2VgEwNvdcw6cfCQi4Wgz0N29CbgO+DuxoP6Duy83sx+a2cXxZteb2XIzWwxcD3whWQVLG/bu4LlesUDPrPvgtVx6NZ4KwLO9esK2pUe9NBFJrjanLQK4+7PAsy3W3dTs8XeB7ya2NPkwfMNrvJmby6T8E9i7vIJz75/Nmrt384meg3hsxuU0Fg7klUgE9u1JdakikmA6UzQwK6N7ATgnaxI3Pfcz8jevo8fUKYwtLeH6J29jZq+T2ZWZya73Xk9xpSKSaAr0wMyrehfcmfazF+nRWM/cq3/A0Dvu4KlLv0GfvZVc+OCqWLumnSmuVEQSTYEemBIaOGOpw6YtzBs3gx3DxwGwbuwUlo6eSs7iNYzf5CxsUKCLhEaBHhBvamBxdgZXvhab2vLIyZ+ie0UDi1/cRK/qJp467QoArv27UxIfmhGRcLTrQ1HpGtaXvcH4DUafqgg5V1zNuTtyyN9WzWuLqjkWyO/Ri25nncegeX/DK7Oo3rOd3r0KU122iCSIjtAD8vLGuXzqjSiNWbm8Wncq/Rph18Q8vviT0yg7Jof8Wuft/E8TycjikreivFH6TKpLFpEEUqAHZMV7izmuDEqm3cieykbm9Y/Q0L2cXVvWUJVXyeJBGeyqaGBZ0Vc5Y5mzoOzNVJcsIgmkIZeADHx9E5uHnE5N1gCOmVDHie88Su+1lTz+AowBqnsWMui4S9m2biwV+ZPIeWUlXJTqqkUkUXSEHojK+kqmLXTePfYzRJtWsPr1e8jwKDWnXMal3/lPtk64kO4NNbxXch+RxvWsPP7zjFlQSWO0MdWli0iCKNADMf/tP7Gz8FNEGjfQUPM3jhl2IhNGXM1VO0eT+cQevrzvBMaN/ncK+o2gce/TNFBD98yzWbrurVSXLiIJokAPRN0jf2XToEk07nmKyYM+yrRu53NCNJMF/erYfUo2C3rWMrkxi5l9LuPYfifTUPM460acTdnjv0t16SKSIBpDD0TTtok05v6RovyPMSx3IgsKNvFO9Vqy9kZYtiTWZmn3bKbnjuGE6BlkR7uxouYZmhbmpbZwEUkYBXoAdpeUsLlPDpNyp5LTdzhP5s4ny8sYOmYbhflVDMhzNu/OZndlf1Zs387KnsM4yybREK1nY3YfGnfuJDs/P9VvQ0SOkAI9AP+4fy5Ds6ohfwqv9PkH48cvoj5vHwuiH2Ee59I3qye1ueWM6rmSKcPewHb34fl3y5lh06nZ8g9e/fljnPX/rk312xCRI6RA7+LcnZ01Oxg29ETWH/Nn+o2p5L6ML7KIKZAJ3Rx6RTKpzGziRTuTDJyP9J/PJ4seZ8macib4BazZWMJZqX4jInLE9KFoF/fqPX9ixOAJbD3hUV4dO4WbMm5nfdZUvjmykMLinXyuOoMVp03iwrII41fv4cvDC1mWMY3vZ/2Uxccdw9rxjzN28ASW/0WX0xXp6hToXZi7s2fNdlYX/YVfDPoKc+18rhycz1vTJvLtUYPJrI9i8Xt8G5BbH+WmsUN4c9pELizoxxz7NA+M+CyrpvyVDc8vSu2bEZEjpkDvwl791VOU/PNCftb/BnYygPurnua/jxtB3+z3R9IOBLq9/7yBOdnMPmE0Pyn7DRsYzZ2DvsKKf1rM4j++fLTfgogkkAK9i6rbXcPzQ97gl32voVvTPh4q+TIXZe7q0D4+3zvCbStuoMG7cdfAL/HXhjlEGiNJqlhEkk2B3kXd+sLPub/flRQ07Wbg5v/kzOq1MOajH2jj7i2WW+wkbxCf27GIUSv+Lz2i9dxd+C/c/sSPk1y5iCSLAr0LuunZ+3hwwMUMbdzK1+79b07t2ye2IX9sx3Y06gwygWGDs/nB3ffR16u4e9CF3PXCQwmvWUSST4Hexfxk/hxm557KsMgm/v33i5h/bBkfq60Hy4SC8Qe13z92vn8s/QMKJwJwdvYAVg5ezDVPLSTfd/LjzBN5cOm8ZL4NEUkCBXoXcueyV7ijdgQjIhu4Zk4JU5e+wJvHGTPLN0HeYMjI/ED7liMsB+nRH3L7c+GOLfz9ZGPqgjf4xtz5FPh2btzRn0fXLkzaexGRxFOgdxF3rl7A7RW9GRldx7deXE637Xm8eOI2Ti48mZyq92DIlDb34a1FfOFEBmwqZlD/Y1g/8F28LJ/vvr6YQt/KN97L4LGNK5PwbkQkGRToXcCPVi3m9i2ZjI2+y7dfX8e2rTlMXP0ET5+awcX9YsMmHHdhq8/dP9BirYy4AAf+R3DRkNP47UyYtPrP7Cyt5dvFJQzxMr6+bh8PbHg3oe9HRJJDgd6JRd35+tISfrnVmRhdwvVvr6d8zQaG76pl9fhG9vQwLqzaHWs8/NQP9yLHngfAZdFctvc3qgZWMqC2L3VLN/G1kmJG+Vq+v76W297VkbpIZ6dA76RqmiJcWryAx3bA9MirXP12GXWrthDtfi7jSp/gnn+u5fRhp5NT+hJk5kD/0Qfto+U0xYOmLQIMnwbAgLWvcHz/47nzzL1MXPk79nX/GDlL9zJr0UJOjC7iF5v38aWSEvZFo0l4tyKSCAr0TmjV3jpmvrmAN/Zk8YnGJ7ns7RqyS6uozi5i8vIHWTGjkN15xg0nzIKKlTDunMOMqbQhMyt2dL/2Ra6f8lU2FRgVE/tw/Jo5lNs48ld144riUs6O/I2/7Iaz31rAe3X7EvuGRSQhFOidSMSdezaW8bG3l7OrsZGv1N7NWW8OoMe6OrZGhzFm23Ly96zh9mnbmJA/gTHr34g9cfLlh9yn2cGn/h/kuNidomdUbmdg7kBuPmM7Q7a9yYia3axv7M6gdb04d36EL+67n/X1cOb8JTyypfygE5dEJLUU6J3Em5V7OGf+Am5dt4PjfCnf3HYXH3nzdDLf28bGpj4M9lxGlj7Nw/96DI1Zxg//6Yfw8k/AMuDY81vd50Fnih7qxU/+AgA2/15unXErVT2NeZ+fyLhlD1HQbSSr6xvpu7GJU944ke/s/ClDomu5YfUWPl5cwuKa2sR1gogcEQV6CkXceXFHFZ98ZwGfXFTK1toqrt53H58uWcTJSy6jYtNbbI4OZVBWIRPm38maf/soz/TbwCfGfoLx5WuhbhdMvTI2bHIkuveGiZ+Esnc4zbszY8gM7h26ml2fmcmk12+jIHci6xpyqd+0gpMWXspnl/2DzzY+yKqaPZxb/C5XLirhjd17dMQukmLtSgIzOw/4BZAJ3O/ut7fYngM8BJwM7AQ+4+4bEltqGOojUYqr9vLcto08s2MP5ZHu9PVq/k/kOY7fuIWhG84kb3cdC2ueJ9rtDEbWbmfMottZe+V0fjDgZYoKi7h50rVwz/TYDs9t37VXWj1TtLmP3gzLn4In/5Wf/9s8PvviLL7CK9x98SlMfvpWVp9yHZszBlK1+Xkm7v0IBTujjBn1PywZNoa5u8/lhcpShmfV8fGBfTm3cDgn9e5JdsaHHNcXkQ+lzUA3s0zgV8A5QBnwjpnNcfcVzZpdDex297FmdjnwX8BnklFwZ7YvGmV3Y4Sd+2rZUV9JeW0V2/dWUl5by4b6COsac9hIHxotmyxvZCIruKh2ISPLahm0eTq9K0extmoB6zNG0zv7LI5d/jg9M7fw8FVD+Ovg+Zwx6FR+WngWWb8+J3Z0/i9PQE6CbvLcfxR8/C545npyf3sxvznzP/jquj9w3XFv85keQ7j4mfsYkDOe0uM+zcKq1fSp2sLEyrMYtK6JKcPvZ93QvhR3L+KezRP4ny2ldPNGxmRWMTbHGdurJyN79WZYz74U5vahT3Z3emdn0s3swBi/iBw5a+vPZDObDtzi7ufGl78L4O63NWvz93ibN80sC9gGFPhhdl5UVOTFxcUdLviW2T/m6TGxk2Ecw5sdebZ8HHtxiy/HHnvzZbNm48rW7Dn7l5s///D7i5BJg+Ucsu7eXslQNjG8sYwRe7YzpqKB3ttH06Myl7qaHWxtbAQGUlBdy+DyBVRmL+flScZLk40R0QhXVdVyyd5KDNhmA7kj9zoWZR3+7NDSij1cNX0kt1w8kWt/v5DnV2xjZH7Pwz7n9MbXuLZ+Nv29kgjwSF4+j/TJZWfUOHuRc/pS6G5T2DZwCjt6ZZPpFQzJySW7V3/q+tewq3Aj6wd0Y2OvIWzKGk4Zw6m1Xq2+lnmUbjSQRROZNJHpETJa/BQ6i8TWlJh9JfJ/hene50f7sGJa+SruueKGD/VcM1vg7kWtbWvPkMtQYFOz5TKg5VksB9q4e5OZVQH5wI4WhcwCZgGMGDGiXcW31DMKg/dVxPbHB3/o5u/HbGyrf+AHZd4ysgFvGd2xdvv38MH9c/A+cXDIwOkZqaNHpJ7cxiZy6yP0rHN61Rk99kK3Gidam0FkX5SmSC6Nns3uaCnV7CWLarr1qKEmr54VY7N4u3cWed1GMTbSnRl7u5Pn3bBs48kBg1ibewJlOWOIWhbj2uirYwvz+PjkIQBcVjSs9VP/W9jKedwSPYvja4sZtq+UYyK7uLmqlm2Z+9gwsZ7XJzURqdtA//I15FU5OfX9qawdgNf2orE8i6x3uzMpqyeTs/fh3ddB3jLq8xrZ3Tub6tws9uRkU5eZTUNWFg0ZWTRmZNJkmUQtg4hlEv3AEfuR/EIn7lc0sbHSRl3tLrtzvr/E7bPzvb82f3YdkFffmLB9NdeeQG/tXbTso/a0wd1nA7MhdoTejtc+yLe//D2+/WGemOZmjh/IzPEDO/CM6UmrRUSSoz2zXMqA4c2WhwFbDtUmPuTSB+jY7XNEROSItCfQ3wHGmdkoM+sGXA7MadFmDnBV/PGlwLzDjZ+LiEjitTnkEh8Tvw74O7Fpi7929+Vm9kOg2N3nAA8AD5tZKbEj80OfuigiIknRrnno7v4s8GyLdTc1e1wPXJbY0kREpCN0pqiISCAU6CIigVCgi4gEQoEuIhKINk/9T9oLm1UAGz/k0wfQ4izUTkJ1dVxnrU11dYzq6pgjqesYdy9obUPKAv1ImFnxoa5lkEqqq+M6a22qq2NUV8ckqy4NuYiIBEKBLiISiK4a6LNTXcAhqK6O66y1qa6OUV0dk5S6uuQYuoiIHKyrHqGLiEgLCnQRkUB0uUA3s/PMbLWZlZrZjamuZz8z22BmS82sxMw6fm+9xNXxazMrN7Nlzdb1N7O5ZrYm/r1fJ6nrFjPbHO+zEjO7IAV1DTezl8xspZktN7OvxdentM8OU1dK+8zMupvZ22a2OF7XrfH1o8xsfry/Ho9farsz1PWgma1v1l+Hv29j8urLNLNFZvaX+HJy+svdu8wXscv3rgVGA92AxcCEVNcVr20DMKAT1HE6cBKwrNm6nwA3xh/fCPxXJ6nrFuBbKe6vwcBJ8cd5wLvAhFT32WHqSmmfEbs7Wa/442xgPjAN+ANweXz9vcA1naSuB4FLU/lvLF7TDcAjwF/iy0npr652hH4KUOru69y9AXgMuCTFNXUq7v4KB98t6hLgt/HHvwU+cVSL4pB1pZy7b3X3hfHHNcBKYvfITWmfHaaulPKYPfHF7PiXA2cBT8bXp6K/DlVXypnZMOBC4P74spGk/upqgd7aDatT/o88zoHnzWxB/GbYnUmhu2+FWFAAHbm5aLJdZ2ZL4kMyR30oqDkzGwlMJXZ012n6rEVdkOI+iw8flADlwFxifzVXuntTvElKfi9b1uXu+/vrR/H++rmZ5RztuoA7ge8A0fhyPknqr64W6O26GXWKzHD3k4DzgWvN7PRUF9QF3AOMAaYAW4GfpaoQM+sF/BH4urtXp6qOllqpK+V95u4Rd59C7P7CpwDHt9bs6FZ1cF1mdgLwXeA44CNAf+A/jmZNZnYRUO7uC5qvbqVpQvqrqwV6e25YnRLuviX+vRx4itg/9M5iu5kNBoh/L09xPQC4+/b4L2EU+F9S1Gdmlk0sNH/v7n+Kr055n7VWV2fps3gtlcA/iI1V943fIB5S/HvZrK7z4kNX7u77gN9w9PtrBnCxmW0gNkR8FrEj9qT0V1cL9PbcsPqoM7OeZpa3/zHwMWDZ4Z91VDW/ifdVwNMprOWA/YEZ90lS0Gfx8cwHgJXufkezTSnts0PVleo+M7MCM+sbf5wLnE1sfP8lYjeIh9T0V2t1rWr2P2UjNk59VPvL3b/r7sPcfSSxvJrn7p8lWf2V6k9/P8SnxRcQ+8R/LfD9VNcTr2k0sRk3i4HlqawLeJTYn+KNxP6iuZrYmN2LwJr49/6dpK6HgaXAEmIBOjgFdZ1G7M/dJUBJ/OuCVPfZYepKaZ8BJwKL4q+/DLgpvn408DZQCjwB5HSSuubF+2sZ8DviM2FS8QXM5P1ZLknpL536LyISiK425CIiIoegQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEP8ferJcFyf58BQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "u0 = uini.copy()\n", "\n", "# total simulation time\n", "TotTime = 7\n", "# Times of outputs\n", "nsteps = 10\n", "\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "\n", "for m in range(int(TotTime/dt)):\n", " u0, u = do_timestep(u0, u)\n", " if m in range(0, int(TotTime/dt), int(TotTime/(dt * nsteps))):\n", " ax.plot(x[1:-1], u.copy()[1:-1])\n", " \n", "plt.show()" ] }, { "cell_type": "markdown", "id": "eP2n76VR4qNr", "metadata": { "id": "eP2n76VR4qNr" }, "source": [] }, { "cell_type": "markdown", "id": "6mIxUuiw4qNr", "metadata": { "id": "6mIxUuiw4qNr" }, "source": [ "### Diffusion in two dimensions" ] }, { "cell_type": "markdown", "id": "lAH8LeIq4qNr", "metadata": { "id": "lAH8LeIq4qNr" }, "source": [ "An example of a two dimensional diffusion process is the diffusion of heat from a circular source across a rectangular plate. The initial temperature is `Tcold` and the temperature of the heat source is `Thot`." ] }, { "cell_type": "code", "execution_count": null, "id": "qm87Vpc04qNr", "metadata": { "id": "qm87Vpc04qNr", "outputId": "c7cf90a4-30a7-4be1-b2ea-bf594e46fda6" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 1\n", "10 2\n", "50 3\n", "100 4\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAERCAYAAABb1k2bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2de3xdV3Xnv8uyZdmyHdmxY+NHkiZxQ0NoTEhDhpRHQwtJBgj9DJ2mQ3gNUx4Fpp1paaGfebTQfGbotKUwMJmGMjSQtJk0MwyB8ialDVACTkgggQQ7xsGKYseyIluWYkmR1vyx99LZd+vceyVbj3t01/fzOZ9z77nndaW1f/d31n6JquI4juNUh2WLfQOO4zjO7HDhdhzHqRgu3I7jOBXDhdtxHKdiuHA7juNUDBdux3GciuHC7TiOUzFaUrhFZIOIfEpEhkXkURH5Vw32FRF5v4gcicsfi4gs5P067Y2IXCsiP4zx+oiIvKDBfg+LyFEReUJEbhKRdQt9v071aUnhBj4CjAGbgdcAN4jIs+rs+2bgVcBFwM8CLwfeshA36Tgi8kvA+4E3AmuBFwL76uz+DeByVT0NOAdYDvzRQtyns7RoOeEWkW7gXwD/UVWPq+rXgTuA19Y55PXAn6pqr6o+Bvwp8IY6536xiPSKyO9Gx/O4iLxKRK4WkR+JyICI/H6y/6UisltEjonIIRH5szn9ss5S4A+B96rqt1R1UlUfi3E4DVU9oKr9yaYJ4Lx6JxYRFZHfEJE9IjIkIu8TkXNF5J9iTN4mIp1x340i8lkRGYxxfJeItFz5duaG5Yt9AyX8NDChqj9Ktt0PvKjO/s+Kn6f71nPnAFuALmAbQeA/CnwZeC5wJnCPiNyqqvuADwIfVNVPisga4MLZfx1nqSIiHcAlwB0ispcQV/8PeJeqPlXnmJ8H/g5YB4wAv9zkMlcSYnMHcC/wfMJT6BHgn4BfA24CfhvoBTbF4y4DfDyLJUor/iKvAY5m244SHkNnsv9RYE2DPPc4cL2qjgO3AhsJ4jykqg8CDxJSLrbveSKyMbr/b83+6zhLmM3ACuDVwAuAXcBzgP9Q7wBV/XpMlWwH/huwv8k13q+qx2JsPgB8SVX3qepR4PPxehBi9RnAWao6rqp3qQ9EtGRpReE+TnAjKeuAoRnuvw443iBoj6jqRHxtruhQ8vlThB8DgDcRngAeEpHviMjLZ3D/Tvtg8fPfVfXxmAb5M+DqZgfGdMoXCOahEXls1ovV/wbsBb4kIvtE5N0zuH+norSicP8IWC4iO5NtFxGccBkPxs9nsu+sUNU9qvprwBmECqjbYw7ecVDVJwnpiZN1tsuBc+foXoZU9bdV9RzgFcC/F5GXzMW5ndaj5YRbVYeB/wu8V0S6ReRy4Brgk3UO+QQhSLeJyFZCru+v5uJeROQ6EdmkqpPAYNw80egYp+34OPBOETlDRNYDvwV8tmxHEXmNiJwZm7CeBVwPfHUubkJEXi4i58UU4TFCnHqsLlFaTrgjvwGsAp4A/gZ4W8zxISIvEJHjyb5/AXwG+D4hB/h3cdtccCXwYLzeB4FrVfXEHJ3bWRq8D/gO4Unxh8B3CYJMFOnjInJm3PcC4JuE9N43gIeBX5+j+9gJfCWe+5+A/6GqX5ujczsthnj9heM4TrVoVcftOI7j1MGF23Ecp2K4cDuO41QMF27HcZyK4cLtOI5TMRqOVdIt4k1OWoRhVR+qtgEeq62Dx+r8447bcRynYrhwO47jVAwXbsdxnIrhwu04jlMxXLgdx3HmEBE5X0TuS5ZjIvJbcS7dL8cZjb4cByWzeXM/JCJ7ReR7InJxs2u4cDuO48whqvqwqu5S1V2E2YtGgE8B7wa+qqo7CaNC2pjpVxEGCdtJmEP3hmbXcOF2HMeZP14CPKKqjxKGp74pbr+JMMk5cfsnNPAtoEdEntHopK0456TjOM6CcOWVV2p/f3/zHRPuueeeB4F0eOcbVfXGOrtfSxiaGmCzqj4OoKqPi8gZcfs24EByTG/c9ni9e3DhdhynbenvP8zu3XfP6hiRFSdU9ZLm+0kn8ErgPc12LdnWsEOZC7fjOG3O0/N14quAe1XV5gk9JCLPiG77GYSJYiA47B3JcduBvkYn9hy34zhtjBKEezbLjPk1ijQJwB3A6+Pr1wOfTra/LrYuuQw4aimVerjjdhynjTHhnltEZDXwS8Bbks3/FbhNRN4E/AT4lbj9c8DVwF5CC5Q3Nju/C7fjOG3MJLX1jHODqo4Ap2fbjhBameT7KvD22ZzfhdtxnDZmfhz3fOPC7ThOm+PC7TiOUyEUmFjsm5g1LtyO47Qx85Pjnm9cuB3HaXM8VeI4jlMhvHLScRynYrhwO47jVAzPcTuO41QQd9yO4zgVwlMljuPMgI7sfb2R3iaz99VrbVwFPFXiOE4dTKyXxdcrsu31MLEej8dOZtudU8Udt+M4GSbMK+LSAayO6y4KIU/3nUjW5gcnCMPGTRBEnGxf51Rw4XYcJ2LOegVBpLuBtQThPp0g2mvjvp3ZsWNxPUQQ7yME4R4ChglibgLu4n0quON2HIdCsM1drwV6gE3ABsJUJ93A1riPibflui0dMkQQ5T6CWB8ABoDDwGD8PHXhLuAnwyQwutg3MWuWlHDfOsv9r52Xu3DaFUt1dBGEewNBlLcQRHoHsDO+Pp8g3p0bqc2bmGrH/MhYfxDthwkCvocg4H3AQYJ4D8RTWBWbC/hscMe9aNwKvOI18c3N+UNnHa4bYxj4zC0u4M6pY6K9lqC/O+LrnfH1BXE5dwVwEbA5fthNUHgoaiwtBzIAncPQuQcuOwTcD4+Mww8IywGCkK+Nr1cQhLwDF+/ZMS8z4PQAfwlcSPh1+NfAy4BfJzw0Afy+qn4u7v8e4E2Ef92/VdUvNjp/pYV7umDP4pHnZoCVvAIXcOfUMNE+Lb4+n2Cgn09w1y8Azl4NXALsijtcQBBsy5d0ZycdpsiTDBCU+mE49z44dzc8ewTuIvwofDMebiLu4j0b5s1xfxD4gqq+Os72vpog3B9Q1T9JdxSRCwjy8yxCRHxFRH5aVev+Cysp3Kck2DWMuoA7p0Qu2s8miPYrCRWQV68mOOtfJHivXQTR7jyTUE7XABsJRdGK49PQ+TTQDxuPAw/C838SxPs+4Ofg7K/A2XvgcyMhFfPZeNofUIi3MxPmvh23iKwDXgi8AUBVx4AxEal3yDXArao6CvxYRPYClwL/VO+Aygn38JwIds50Ae++ZY5O7Sx51lIr2q8jOOArLiL4p38JPBN4HiCXEgT7fIJgm2h3ZWc9QXCC/cCLoPNh2PUgXPRtuJvwI3AbXH0A7nwAXg3cTviN2BPPZs0InUaclOPeKCK7k/c3quqNyftzCOmQj4vIRcA9wG/Gz94hIq8DdgO/rapPAtuAbyXH98ZtdamUcBeirfN0hSjg1wnDr3HxdhqT1ilaeuR1BCF/8S8SUiH/jiDeZ5pgXwScTUhy98SjrBgmjntqfYLQhuSZ4Vh5Flz2IGz9dpCHj8AVq+Fr34aXE5z3VkLOO20P7jRi1sLdr6qXNPh8OXAx8E5VvVtEPgi8G/gw8D7Cr8X7gD8l5L7LrHhDkauMcM+/aCfcrC7eTlOsjfYOivRIN1G0NxGK5mnAxjcC24HnEgR7I0G0y5x2zgkKZ76ZIPrb4cxzYfXfwHuAP4YXjwTn/TxC7tvaf0/gwt2Yeeny3gv0qurd8f3twLtV9ZDtICIfJfzO2v47kuO3E2o36lJvmISWYkFF24jXmm0TQ6c96IiLNfl7PiGnfcVFceP7COK98V0EP351XJ9HyEqviUtXk8X22xKPtXNdBBvfFlz3W4EL4Io1Ic99DsF1dxE69jTrVt/eWKpkNkuTM6oeBA6IyPlx00uAH4jIM5Ldfhl4IL6+A7hWRFaKyE8RMl7fbnSNlnfcU8K5kKI9dc1OXsEYuOt2EqyDjXWusbbZV6+OL/4dwWmvexfBJb+MIKNbKCohZ1r0bN+nKRz6coIW/AOsey1c8MlQvTUEF38+9LI8HN4yhrvuxsxbq5J3ArfEFiX7gDcCHxKRXfGi+4G3AKjqgyJyG6Fu+Wng7Y1alEAFhPsVr2HmbbPnnFHAUybOdFZQdK7ZQWjyx05CReRWYnrkdIJom2OejWDnpAJu57o8vN/4fPj5b8JDsKUPzrk/VJQOEMR7FG8e2Ji5F25VvY/QADTltQ32vx64fqbnb+lUSZGmWMQuqZ4ycRLMbXcRstRbCemJs1cTmvw9k1gRuR14EdOd9qlirtvy3hcCL4EzN4QMygVw+eqQCT+HYmwU64Lv5FiOezbL4tPSwr24bjvh5s6i3bjT9qwgVEJuougVySUEDX0ehNYjzyVI+1yKtmHi3QOcG5drQgr8pcD54Xa2xqt3U3TKdHLmPse9ELSscLeE254i3IO7bseGYV1LqIPcSezGvisucilFn/a8Y81cspyi4vJ84FzYeGlw/C+Ac5eFlz0U7cxbtrAvOhOzXBaflv1ftozbNtx1tz3WkmR1XHYQXC0XUXRj51mERMVGattozwfLCdLcQ7D754fWw+cBF4d7Oyu5X7t/J8Ud95zRWm7bcNftFLni0wkpiPMhmOupbuznU3SumS+3bdj5N8Y7ehZ0bg01k5uLEQhPp3bWHSfFc9xzSyu5baMV78lZcKy35Fbi2FA7iSP8PYvazjUL0WgrFe+NwEVTORwb83sZzbv5tC/VdNwt3xzQcVoJyxVb3rhzI4VCTg0YZZ1nFkq47Xobwz3EX5TOjdDRXww162mSMnw8bsdZ0tjjaTr29lT3yQ6Y38rIZiTOO7knm10nnay4NarXWgWlVdIfs8GF23FOgqmkmeUgumFmY4/MJ9HldydvmT6fpZPijttx2oLUvbKMpNYvH+VvIcmuvQJYNv0pwclx4Xac9mOy+S6LRivfW8vgwu04bYHliKd00eaIrBlHe6HJrh3vye7R89r18By347QNY/bCyvwwcbqxxRSBOGvOcPKW5F6dEtxxO86SJnevQ/ZmwDb2s3htfe26/TX3ZHNPTntKcCIu3I6z5JkgiN9QfD3WD53DhPlKNh4niLf1sFuIpoHm8k/Eax8P9zIc7m2C2tlwnJxqCnfr9py8rgUf8FrxnpwF5wRBCKM+htl5BwAeJIjnIAvnvBOnTT9wf7iXPeHe+pinybmWDEoYzmI2y+LTksJ97dSrlYt4FznhXq5tspeztBknONcjBGF8GOAQYe6SsZ/ELYdYGPFORfsI8CCM9cH3wy08HO/xSLzn8brnaWfmp8u7iPSIyO0i8pCI/FBE/pmIbBCRL4vInrheH/cVEfmQiOwVke+JyMXNzt+Swg3wmVtoLYd73Vi4J6dtsUE9R+JygDij6/0ElfwBBNe9nyJlMt/CPRiXB8JN3AvsBe4N9/Zocr+tMyhpKzFvY5V8EPiCqj6TMH7kDwkzvX9VVXcCX43vAa4ijHizE3gzcEOzk7escLeW63a37QQsVzzEVEaCR8aB++Ki3yYo+SHmt7LyaeB4XB4GHoH+b8NDwF3wyGR4OUiRj/eKyXrMrXCLyDrghcDHAFR1TFUHCTOD3hR3uwl4VXx9DfAJDXwL6MkmFp5Gywo3tJDrdrftJIwTUhCHCa77BwC7Cab3bgiu+x6CbB5k7p23ifYg8Ehc/k9w2l8CHg630xevPoynSeozL8O6nkMIj4+LyHdF5C9FpBvYrKqPA8T1GXH/bYRQMnrjtrq0tHC3hOu+TgB3207AcsUnCLLZRxDu/SPAVwg29yffJqjoP8S95jJtkop2L+HX4ovwk2PwuXAz3xgJyZp9BLc9QpGbd3JOKlWyUUR2J8ubs5MuJ0xpcYOqPofw2/lu6iN1bqwuLd8c8DO3wCsYg5sX4+rhB8NneHdyxgmieJBgle4Czt4D3EbwW6v/Bjb2EMQ7zsZOD2Ho15NpJmiikYr2PcA3QorkTmA3HLw/CPb3KbLfnttugs76r9OvqvkM7im9QK+q3h3f304Q7kMi8gxVfTymQp5I9t+RHL+dWH1Sj5Z23JA43evKfpTmGU+ROCWY6x4hiPceQin7nNVYfoRQFI/dQMh3/wMhD21N9uyReyYOPG+n3R/PdXc477HPhLzI54E7Q93kvYTn9EOExmueJmmAEv5As1manVL1IHBARM6Pm15CeDC7A3h93PZ64NPx9R3A62LrksuAo5ZSqUfLO24Ijnf4NQTxvrnhE8Tc4SkSpwHmYgcIY15/kzCj+p0PwBWrgT8G3gpc8EnY+HzC03I/YUKx4xRjdzcbBtYE3tqHx4pIvhic9m7gE8DdcOdocNr7CD8kJwjd3d1tN0CZrz/QO4FbRKST8C95I8Eo3yYibwJ+AvxK3PdzwNWE/NpI3LcholpfCLtFFkglZ8awTdZ7cyfz1xB+5VSFaCulSIZVF+GRozosRqx2EGR3JcWM6q8mzDj54gsJ81BeA/w8cOaG+OZcwsS+pzN9QmFbpwNGmdM+Qshnx4rInxwL6ZHPA3fD1x4Nkv4Vgul/jMXLbVcpVi+5WHT3XbM7RtZwT5NUybxTCcdtdN8SJut9BWPAPLjv6LI/c4s7bac5EwRZXUFIl1xASGa+nOi89xNyKQ8BVw/AeR+HjZcSJNbmp6znvFOn3U9oqfJwcNl7CR5tN3BncNp7CKJ9hJB3txSJu+0mzJ/jnlcq5bhT5tZ9t6bLTqmSi1kMFitW02nMVhB6UPQAzyMI+cXAlovim5cSrPnFhNnYuYhinsq0wjLtEXkcuD/0iLyX8CPwJeAHoSLyXkJ6ZDdBtPcRJN863CwGVYrVS3aJ7v7q7I6Rje64T5rp7vtkBLwQbHfZzskwQRDvIYJ47yHM1XsXIV1xBDjnfrh8D0F0XwB8B3h2H2zoCzt3UEw3ZgzHk/cREunfJzjtu4CHQ5O/fQTh3keoiLQW44sp2pWkgn+sygo3RKEtFfAZ4ILtzBGpeHcRcsxDBBE9DDwb6B2BS74L595PcNybCfa8mzCxLxRToFnLhQGCgO8hKPO9oUfkbkI77e/H8/cRqi1H4zUrqEOLh1LJLqWVFm4jFXBumXlPSxdsZ64wsTThnKRoLjhAEPNHgGdOwtbdoW1JN9C5kaKWcxmFiMQTjfUXg1n1EUx7H8FlDxL03BoLek77JKngH21JCLfhQuwsNhPJepzQHG+IopflXuAs4OvELEl/SLFA0anCtNvGGLHhYx+l6ERv5xyN1/Dxtk8Sa8ddMZaUcDtOK2BtvE1MRwlCe5Qg0vsJjQGXUYh2nuCz50abBOEIhYMfpshju8s+RSraqsSF23HmCRPwDoqmg4PAaoKLtuyItUzpSI6ztQ2BZMPJmli7YM8hnuN2HCfHBHacINTW9hsKsZ7JsenQrC7ac8QklZxN2YXbcRaA3EWnadV6AwblRtDFep5wx+04TjNyAXZBXkQ8x+04jlNBXLgdx3EqhDcHdBzHqRjec9JxHKeCeKrEcRynQlS0crLlpy5zHMeZN+Zh6jIAEdkvIt8XkftEZHfc9gci8ljcdp+IXJ3s/x4R2SsiD4vIy5qd3x234zjtzfw57l9Q1f5s2wdU9U/SDSJyAWGopWcRhrD5ioj8tGr9WYzdcTuO075Y5eRslrnnGuBWVR1V1R8TxiK7tNEBLtyO47Q3E7NcZoYCXxKRe0Tkzcn2d4jI90Tkf4nI+rhtG2HkX6M3bquLC7fjOO3LyeW4N4rI7mR5c8mZL1fVi4GrgLeLyAuBGwizRe8CHgf+NO5bNtVbw6n4PMftOE77cnKtSvqbzTmpqn1x/YSIfAq4VFX/0T4XkY8Cn41ve4EdyeHbCQNI1sUdt+M47c0c57hFpFtE1tprwjTRD4jIM5Ldfhl4IL6+A7hWRFaKyE8RJrX7dqNruON2HKd9mZ9hXTcDnxIRCBr716r6BRH5pIjsIvj8/cBbAFT1QRG5DfgB8DTw9kYtSuykjuM47csctxRR1X3ARSXbX9vgmOuB62d6DRdux3Hal4r2nHThdhynvfFBphzHcSqE4lOXOY7jVAof1tVxHKeCeI7bcRynQrjjdhzHqRie43Ycx6kg7rgdx3EqxFJsxz2sWjZqleO0HB6rzknhs7w7juNUjKXouB3HcZY8nuN2HMepEO64HcdxKobnuB3HcSpIBR13S8yAIyI3i8jjInJMRH4kIv8m+axTRG4Xkf0ioiLy4pM9l+PMBSJyrYj8UESGReQREXlB3H6ZiHxZRAZE5LCI/G0260l+nq+JyAkROR6XhxfuWzhAkSqZ+8mC55WWEG7gvwBnq+o64JXAH4nIc5PPvw5cBxycg3M5zkkjIr8EvB94I7AWeCGwL368HrgROBs4CxgCPt7klO9Q1TVxOX9ebtppzBxPXbYQtIRwq+qDqjpqb+NybvxsTFX/XFW/zgx+7xqdK0dE3iAi3xCRD4jIoIjsE5Hnx+0HROQJEXl9sv/VIvIDERkSkcdE5HdO4Ws71eQPgfeq6rdUdVJVH1PVxwBU9fOq+reqekxVR4APA5fPxUVF5A+ig785xt/3ReSnReQ9MU4PiMhLk/3fEON5SER+LCKvmYv7WHKc3CzvTYkZgu+LyH0isjtu2xCfyPbE9fq4XUTkQyKyV0S+JyIXNzt/Swg3gIj8DxEZAR4iTF3/uQU61/OA7wGnA38N3Ar8HHAeweV/WETWxH0/BrxFVdcCFwJ3nuw9OtVDRDqAS4BNsZD1isiHRWRVnUNeCDzY5LT/RUT6o4F4cZN9XwF8kuDsvwt8kVCGtwHvBf4i3mc38CHgqhirzwfua/oF25H5TZX8gqruSmaEfzfwVVXdCXw1vge4ijBB8E7gzcANzU7cMsKtqr9BePR8AfB/gdHGR8zZuX6sqh+Pk3P+b2AHwVGNquqXCEPQnBf3HQcuEJF1qvqkqt57svfoVJLNwArg1YTY2gU8B/gP+Y4i8rPAfwLe1eB8vwecQxDeG4HPiEjp02HkLlX9oqo+DfwtsAn4r6o6TjAcZ4tIT9x3ErhQRFap6uOq2uwHpH1ZuBz3NcBN8fVNwKuS7Z/QwLeAnkZ1I9BCwg2gqhMxJbIdeNsCnetQ8vqpeGy+zRz3vwCuBh4VkX8QkX92KvfoVI6n4vq/RzHsB/6MEBNTiMh5wOeB31TVu+qdTFXvVtWhaBJuAr6Rnysjj8v+ZDZwu7c1qjoM/CrwVuBxEfk7EXnmDL9je3FyqZKNIrI7Wd5c58xfEpF7ks83q+rjAHF9Rty+DTiQHNsbt9WlVZsDLqdOXnoxz6Wq3wGuEZEVwDuA2wgO3WkDVPVJEeklFMpSROQs4CvA+1T1k7O9BDAnY66o6heBL8Y0zh8BHyU8JTgpJ9cBpz9Jf9TjclXtE5EzgC+LyEMN9i37n9eNMWgBxy0iZ8TmVWtEpENEXgb8Gkn+WERWikhXfNspIl0iMu3LzuRcp3CfnSLyGhE5LT6aHqNlGgc5C8jHgXfGWFsP/BbwWQAR2UaItY+o6v9sdBIR6RGRl8VYXh4rD19IyFufEiKyWUReGXPdo8BxPFbrMw+tSlS1L66fAD4FXAocshRIXD8Rd++l1gBuB/oanX/RhZvwy/I2ws0/CfwJ8Fuq+ulkn4cJj4LbCIH9FKG5FSLy+yLy+Vmc61R4LbBfRI4RHkOvm6PzOtXhfcB3gB8BPyRUEl4fP/s3hJz1f07aZh+3A7NYXUFwwoeBfuCdwKtUdS7aci8DfptQ+AeAFwG/MQfnXXrMQ+WkiHSLyFp7DbwUeAC4A7BWaq8HTJfuAF4XW5dcBhy1lErda6g2dOSO4zhLlkt6RHfPMoEkn+WeRqkSETmH4LIhpGr/WlWvF5HTCenVM4GfAL+iqgMxe/Bh4EpgBHijqu5udA+tmuN2HMeZf+ZhkClV3QdcVLL9CPCSku0KvH0213DhdhynvWmR3pCzwYXbcZz2xYd1dRzHqRhLcVjXbhGvuWwRfE7Fxnistg6Vi1V33I7jOBVC8Ry34zhO5XDH7TiOUyGWYo7bcRxnSeOtShzHcSqGC7fjOE7F8FSJ4zhOBfFWJY7jOBXCUyWO4zgVxIXbcRynQniO23Ecp2J4qsRxHKeCVLByshWmLnOctqLjJBdnHpiHqcuMOO/td0XE5iT9KxH5sYjcF5ddcbuIyIdEZK+IfE9ELm52bnfcjjPPzER0cwdVZgLz81TwCb/1mN8c928S5iVdl2x7l6renu13FbAzLs8Dbojrurjjdpx5oMwpL2uw5MxkX3fjc8Q8OG4R2Q78c+AvZ7D7NcAnNPAtoMdmg6+HC7fjzBHNxDrfx5bOOku9dEkjEXdmiQ3rOpsFNorI7mR5c8mZ/xz4XaY/PF0f0yEfEJGVcds24ECyT2/cVhdPlTjOHJCL5rKS7bmgNyPdx0p/B4Xpm0j2ST8HT6PMGAXGZn1Uf5NZ3l8OPKGq94jIi5OP3gMcJPwu3wj8HvBeoGziiYYTgywp4Z6t4/Dgdk6VemLcka3LhLxsfyONzUlqBXlZsq0shvP9nSbMfauSy4FXisjVQBewTkRuVtXr4uejIvJx4Hfi+15gR3L8dqCv0QUqnSpp9hjZLE/otfbOyVKWEkm3dwAr4vYVJUsXwXZ1xWVFndedybZ8sXOXlYF69+lMZ65T3Kr6HlXdrqpnA9cCd6rqdZa3FhEBXgU8EA+5A3hdbF1yGXBUVR9vdI3KOe5Gj5t5gDZyMfXcStm+jpNST7BtvSx73ZG9LtuvDIvBSWrTI/benDfZ65SyFItTsMD9b24RkU2E1Mh9wFvj9s8BVwN7gRHgjc1OJKr1UymtNAFrs8JStr0eaYFI3+fb022LTeUmYF1gFipW68VaLswrqC/c+bb83HmaxJzeZLZOX48nx5YJPcn7+aZKsXqxiN41y2PWwD2NctwLQcs77lyw6xWU/DPbllIm1JOEQjZB4Uo6ks88V+gY9UR7RbbujPvYOk1n5IJez2zkbjsV6HQ9Sahb62B6HZvFstEoL97OVLDjZGsLdy7K+ftmrqaMZq4lfeTMBdwDvn0pE+0yF20x2UUh3iuS7fk6Nx9QazDqCfZ4XMbicePZepLaWJ5JpWY7UtGhSlpTuFOXsCLZlrqXFdQ6GnvdURU9SskAACAASURBVHJM7rRztzKevZ5IthnuvtuXZqKdxqQtVqm4jKJysTPbpyxWoTbVkQu1ifU4cCIen27rSI7JcfGejgv3HFHmsu11WUHIa9xTMU/PlxaEVKRtOUFtgeiMa7t26s496NuHeqJt8ZU766466+7sfSPnXc9pn0jWJ+K5TpR8VpY2MVy8a6noqK6tJdx5zi9tUmVNo1ZT21wqX8oKRT33cqLOsgIYjdc9wfQcmAd9e9BItNMcdirSaYyuTd7Xi9M0/50/HeaxmsfsSFwPx/OMUMRses9mSEjO7eJdUMXv3zLCnYp2XqljbsUKwNq47o7buqktHHmBgPKCYME/XLK2zzrje8sfWkHwoF/azES0TayXEWKwi9oY7aSI1bXUxnCePilLleRpkVSsR+I1h+K5UvHuiO+hVrBdvKdjPd6rRksId5lod8Xt5l564jZbbyAE7tq4mIBb4UlTK/aPOUEI0NSpjBCCfwg4BhyNxw/FezkWzzESzzEet7t4L11mI9qWBllJGALOYnUl02PWRPxUTMYJQmwOU4j2sXiOIWrTLsPJOY1cvJ1q/h0WXbgbiba5lA1x26ZkvRY4Pa43xHWZ804LQ5nTPgoMEoL+CKGQHY7nOhLPMRjPYwXDcuAu3kufvL11LtoWo/WMxWlx3UOIzfRpcSW1FZYSL6KTtRWRo0x/KhyM5zoaz3OUWvE+Fu/XnLdRr5NOu8aw57hPgtQN5PlCe9Q0sd5MCPathAKwLe6zmeB0NgCdK6lVb3sWjRZm1TisS1R7bBQGCIXgcDzHwXjeQ/EUh+NpBuJ9WkGwSiBo78BfapT1C6hXOZ6KdirSmwgheHpcb2K6eIvFappviReUCegcD0v3ODAME6OFyx6K5zgSD10Z1wPUpl0G4/0PU7h4q1SFIl7bOWXirUpOkrQdbOq0TbRXE8S6GziLINJnxX12AOuWAVsIpaaHWkuTJsqtltES2IPQOQRbBmDLYdg8GUZ12UAY8aWbMM7iCoKYQyHeQ9QWAKhmnsxpTlrRbTFqOe1UtM1h29Ngvt4ArFpGbY6vLF8C0/Mkw9AxAuuGYN0gPDUZYnE1YW2Vomm6BUJMHktOmZM3k21Xqlh2F024O5J12gIkF+0dcds5cb2TUFi2rSRYmq2EwmAlxEqS2RBr6AqhENhz5yDh+fIIMADr+2D9YXhsPPw4PBIPPUAoqDZU1wDTuyTbuh0dy1KiWYcvS2tYeKWinYr0lrjeGtfrzV1bUJv1tkqZZnm9ND+yAVYNwLYhWD1e/IDYD0tZ57Nj8VLDFENKT5Ts144xbE2Dq8aiCLfFZ95szyoic9E2sb6AYKzXbyaUlnMoxLuHogSZm7FUCdRWz9vzpuVJ+uJF98O2w7C6PxTMvfG+9jO9mZYF/4pkezsG/lIhr5C013nnGtPZdXGdpkd6KMR6W/xslcVkmtMz8bYnw7RtYNqkJKZJGCGo7yZCzMaKmPWD0HW0vIIzFehG45yAx6s77lmwLFmnj55WEWkFYGfc9uy4XnVefPFMCvHuoUiXWC1l6mQsSlMHM0Rw3Afjsfvj+fbA+m7oenS6e7HTpMFvziUtME616aBWBFPRTttpm+M2p23+YUdcd2wlxGRaIZPmVcxgrKS2V491OBilthlJItqsDsuqFbClv7xVSt46ZTw5dbrk6ZV2Mh+e454heYrEKn3SmnmLb9PkZxMKR+fPxBe2wVTdrPkmYiKxB1hDod4AJ2LN5HFgEI7HROFmiuT2vnjhh0KB2LE3HJn3tEy7yadibvu2U+AvVdJWJKnjtnpvi9U8PTIl2tsJT4PbKaz4OmAjtTWVjRx32gRqMB5jyW27kU7o6IRNMZeXuulcrMeoFfNO3HVDNc3WojjuvEIy7V1mNfBWAXkBscXIzxBKx8VxJ/vA1H3VGRTWu5tCuO0rPk0oBceBo7CmPyybBkIB209ta5QO6JyArT+ubZaV92LLXUy9cZGd1qVebjvvtZu2JllJkfUw8d5GJtpnURgLC02rqcwdtzWrgsIppI57A0G0004LJvbxnjf1TR+AypY0XjspYnYZ7e26T27mssVnQYW7I1tbwchdjJmTncQ84XmE0nExwSHvIgTy+cCaTkLexOzNGoouEGXCPUqwL4PAQVh1CLbthdVjRQm1AhSbEG7uDU+qVq9paUcb38ScTDo4ULsE/lIkFfB6zVQtTWJ9CSwr0rGVQrQ3AWfHDyyPklaidwGr0nlwSkzGU5NFas9E224mzdEBHWOwob/WWIwmSzeFqFudfbu7bu85OUPS3HaaJjHzYRWTZxFie6oi0tIjuygcd+cG4EKCmm+naBOYuu0uQkGwJUmXsIXQ+K8H1j8E5w/UJqvjo+r6IdhytOhdOUTRGaKLwol7rru6lLUksSUV77Qb+2nUOu5VG+LG7YSYPZvCcW+g1nEvX03Im6yJy3KC2YAQUU8Dx2HVcVjVD+tGQsg2aoUyBqtOwIbjtY1RhuN9pgNT1XPd7fjUOF8/WCLSAewGHlPVl4vITwG3EiLgXuC1qjoWZ3v/BPBcQju3X1XV/Y3OvWitSspq6q0C3upyttnz6DOpTY9cAHSeQfieG4FzCYJtBcGeP5dT62Kepnj+PA701+6/5kE474nwn7T8YmyBsvVo4dOtUcpwslue93aqQd4l3LaVxWie404rJtevpDZfsoMg4ibaVg+z3pqUbCHE3mnUDj+VGowThBr0HujshzP6YcXY9F40aZvvEVg/AiOTtQZjJC5lg1t1UNTTpE+K7fDUOM+Vk78J/JCQQAB4P/ABVb1VRP4n8Cbghrh+UlXPE5Fr436/2ujECybc6RgK6Rja9qDYTe2j5w7im3MIgW8VkecTnfZzCcF/PrUFwZqULKdcuK1QpJWXK4v91jwNZw0Uoj0IDIAchS19Ic04QNFV3lxYPgyssdQDf6lhLrtMvC2yzHGn7bjpJjz4pc1LrI+B5bjXWe2kPR1ujGdq4LjpoXiS7IL1B0NvHLPEmWibq+gZKAzGEEV9ZvqUaOmSDoonxnaL1/nq8i4i24F/DlwP/Ps4QfAVwL+Ku9wE/AFBuK+JrwFuBz4sIqIN5pVcVMdtj5828I41wd5M7BFp+UIT73OIOe0LCUF/flxbhWRPcsZUtI1UuK2fpu2bOPKNu2ForBjAZAA4DKcfhA2TRYustdSOzmaPn061SAeTSrflue58jBLzzquWEaz3Oor22luorZBcZxvPotxolOW4TyNElpmMjrDPul6YGKl9MkxHSxuC7kFYOzk9LW5GI+0ab0u7pktO4ntuFJHdyfsbVfXGbJ8/B36X8C+AoGSDqvp0fN9LeD4jrg8AqOrTInI07t9f7wYWrTmguW4LJqtk30B8rtgS3+SNY3kmoXRYeiTt7566F/tqaWGw9fJsydMo/bDlR0GwtxAGLtkMPAan9xd1S4cpukCndUWeLqkeZc1UU7edpvTSoVtXQ9G8xFqLbKS2K+V6y2efRWE0yupj8qdDS+mVVLSv74XxsSIHYs0FbfSpHlg7UDzJWhmz4Rpmmi5Z6pxkqqS/0WTBIvJy4AlVvUdEXmyb61y+2WelLGpzwNTF2KOciTenUXRHM9ey6gxqG8emTjtNfeSiba/LRBumF5R+WHUQNh0rutPHFgFr+8OmPmqHm8hHkXNan7L8dkrquBuJd027bBNvq7DZQPwwjdnTmV4xaVEEtVa6TNSjwejpC4b8WDz1UYrH1rXQPVDcZzoBST55Q7O/0VIW8XlKlVwOvFJErqboaPvnQI+ILI+uezvFSBq9BGvaKyLLCeo3MP20BQuiM6mbyackSwuDuYNOGwjCAt9EfKoAWH7Qmv6llTsW5HmeO9+WHtOTLFaothTNt05jqkB2rywKrN172m0//Y7pd3dan3z87TRGy5x3N3GUvzRv0p2tl6+mcNi2TmPXKiit74GlT8yRp4sdG2O0c3VMw1AItr3ugo6VhWhbXZLFalkuv13Je5I2W5qhqu9R1e2qejZwLXCnqr4G+Hvg1XG31wOfjq/viO+Jn9/ZKL8NC+i481+IMtdtQTaVNzHxXktMJFpe0ALeHgLzx80yxw2F4zbSWnz7tVhDKBw9IdHePVn7I9IFa0drh5mwAlH2HdslT7hUyJur1kudWJV2TR4izZ/Y49iUSOdi3cP0/Hb+BLicWkk1tx07kdED3SPF9dIelfH6q0fD/aa96q3c5a47zXW3Cwvc5f33gFtF5I+A7wIfi9s/BnxSRPYSnPa1zU60KDnutBCkjsZibuqFCXg31NqYNC1S5qyh/Kul6RLDxHsN4bnTzn96WK89Fi5r05usha6jxSw76XdIv5tXUlaXXMhy192ZrKdZ8HS9qpPaJn+pq+6iPIahEG5jgsKgHKcwLT0hpdc9Of368b46jxYVkiuT72DfK//O7ch8mitV/Rrwtfh6H3BpyT4ngF+ZzXkXfZCpfITALiiCz+pkuqAIeNuYttPOxbrR10rFO89527mtBPRA17EiQRjF2xx2PjZQOwd/1clTBR0lr8sqLGuaa5i1tXXN02BeB5MGdzPhfjrbd03teVccm34fcZ2Wr/y7wPTmq/ZZu7huG9OlaixKl/e8+VUq2itsQ9rPeErNZ+KwZ/KVGom3LbFg2D2sLu7J7rWsQKTfrV2Cv8rkP7b1BDyvn5n63+d23D7ohCJW663LlpSy+pl0HV/nyevkfRqj1gO0Xv1L/qTYLqm+KpbTRXPcML0gTAVVOhGfva5b8Qin9jXSc2QFguXTB2TurG36561Jlh5lrjTtPGavpV5PHQuQGmedPiF2UF+401Te8jrHrKw9d/rjYcnreG95Y4D0+6TfsV3xYV1Pgnxwm6lfd+sJMJ7tMMXT+YZT4OkG6+W1XcsSG1333p0lTcP/82S9nWYarzPdb6J2nV836Ukjy2BisnZeybLv0M71MlUsu4sq3DAtBMPamrFCMZU6J6jttp7mAfMKx2bUKyBpe26K9VjydmL670kVf7Gd+uT/13TyDJLXOhkm9p02a8HU9DLWCiRdPx13KIvjPMedLqPZ+xPUxGo66La9J4h28naa0Wj3GK7qsK4L+oRfL2gs5mw4yqnRmtLxKac6JJQJOCXrZpQVjrRAHA+v7R5Ginuye7Vy2ui7Oa1N7rby/1sqbPnMR1Oxms9YYINh14hr2Xq2S14G4joNRis38X0ao+aB6ol1/r6KTnS22LCus1lagUVz3OkYOalROUGcqMYGwLY45XhcbGPuQMoqHMvIa+uhKAh2zmQqeBsHYpzQQ228KBej1BaMVvmnOrPHxudI3+evU8G2pTMd5Gmc2vCcilnrkXsa03tEGvV68p5ITppui+fVY8XclCWzfOQGo8xoNBPwpU4Vv++CC3fuWNIplmxygqkXNoD8MBTTsg9SFAZb8hr5MvHOBTsXbTunnf9IWNv1B+NNDhVann+H2fSuclqXNGVsSyrWY8m6O5+NPQ3LpybDeNrWWYZhQmx1xR3SqsFGwp3Gpp0jBqVd14xF+j55OrQJddKYTY1Gu5qO+RodcL5ZMOG2ipH0vaXk8un1poQyHfz6+GSYamzqA+uIkLc2MUy8U9J0Sl4wbGYcG6d7EI5NFvdik7WeCLeTlY8irZl9R6dapKJd1t05NdijUEyJZINeDyXrYcIkCFPd1a0fQlrsLN9dJtzWS9Im/rB1sqTXG6UIzri2h8VUtK3c2XfKf6jaCW9V0oB0uEhbW8DkhmUYGBuFTpuFfSBZr+mnmJa9n6LBahr0Jt5l6ZLcaZtoJ2LNofi6N7w9nFx/GIZHi3KRinbquBvlEZ3WJR0ZL4/RVLDTeJ0YhQ4boe8YxQ/9BkKQrBsJkyBMG3LYSPsnmLWxlJ1dyYT6CCE24/qpkeJ6Nna8PaWeCOXIzNAY01Ph+Y9Tu1LFH6tFyXGXuW2LfdPILQOE+IxTQ7KZMLHvqkOEwbTKRgM0YbZtac4b6ou2pWBMtB8NheJwvJmDxcdDTBnvKXNjBaLMdTutixmJeqTClqePLV6HgXUmmpsIwWGzsa+N28/ITQYUsVpvWFdL4aWpu35CjB4MiwWjXd/MTnyfxqg5b/sO6Typzf5GSxl33DNkgtBPIHUyaTrbvMUWE80+gsHuIwwfsm0vhXuxx04rBD3J67KeaHmrERPtfuAxwg/CI2HdRxDuvuQ+BsMmu08TbfsuaUF3qkX6NGgDLeUVkalopzG7zhx2LtpTMxiMhZlrpgmzpULKDEjq6wepFe1eeHKyMDZHCGF8hBoRtxjNU3v16mXa8SnRc9yzIHWnaW7b0tmHgc2TsL6P4GL2EwrGfsJs7OsfYrpoW/vYldROppA7bhPtUYr0SCra++GJsVA+DhDWUbyfmCzKid1vKt5VbA/qhP+difWyZFv6Q2ypBlvS1PJTk7DKRPswtaP02ZAJHSNh5pqp3LU98eWTBVvOO30qtKiLon1srHgaTNcm3kdD9YzNr5DXx9iTrpXBNLedN2td6rjjbkL6x7FeWmMUBcLGgz9CMVHB+sMUor2PojCcPxAm9gWmFwSrvW82WbB5e0u9PALshScHgmAfiDdxAHgMJvpC2cgdt+UPrY2sFXafBaea5IJdlt+2FIlNNjMAbDPHYaJtI/StSE48MRJmrpmK1aMUIwamjjt9Kkwd98HgtC0QDybrI0yJuA4Us+6l9Zb2pFBWL2O32I5U8XsvmuNexvTHT2s0dZCg1Y+Nw7bDwB6KAftsLIbznggT+9aIdhyKdSqVYi1N8s41lje0nHZ02k8OwKPhJY8SfiyiiB+Me9sUlPZ0YM2s8kLgVIM0z90oXWItr0copgBbS4iH1cDqcVg/SO0Ej2UzR4+PhZlrOm2usbTSssxxR4Px1EjtHKgWkIcIgm3BeXhqfuua+sq0pWJaqd4oTdIOseyOe4ak4yVYG9MVhMBaTTEJby9hFNXV/bC+G3iIZPCeePBZA2FiX/rjspFiwPpuapsImmhbO1g75tFwtSfGgkjvJwj2HqYE/MnjRcr7MEU9kBWEtKLHKyirSZ4uSVMIY4TQS9MkXQSj0UUQyW7COO2r0oHR0l+E1KEMA2tHwiQIqw5SO30ZTMWqHisuaM1RraWTibWl8uJ7q1O3Xg9pC8U0x13WnrsdnbfnuGdA6mgg9jyj8MxDhAJwKK4fIY5/9mgsEFYQ0uT40FiY2HeVNRO06cds3OI0x20dGaxWpzdEeh9FTvtRgmjvC8tTh8KPyGNxF2t/Ys12zW1bYUi/q1NNUiEbo/bpcJgQWccIsWlj/1lWZEs/dHQmJ8sbf1uzwXXElMpkMZ52Oo6qWWMLNHMLqeO2yvNDhHTeQJHyPhIvU89x5+2520msc+a6rIpIF/CPFPna21X1P4vIXwEvIsgIwBtU9T4REeCDwNWE//YbVPXeRtdYtOaA5rpPEOJ1hCDiRwgxfIBQEPbGfXfshU57pkutjzWJ2nQMeo6F6campUuWU/voeTzU3qTuxXLa+5hy3WO9QccPUDyRHkkum1b01B0Yzml5zFCk/zsT67Q+Jn06XEHhuNNRiDuATX1JhiQXbcuQ2JR8aRPuFckxExQW2QLObLQFogXlYzDRX2i5dX04Sm2spk0By9z2BNPd91JnnlIlo8AVqnpcRFYAXxeRz8fP3qWqt2f7XwXsjMvzgBviui6L0hwwdd3LCAHVSXAIKwjBt4IgmOk4wlt/DKusEFggHyEYbZuNvXsyTDfWlcwKYs9CqeAPU0R56rij637qUJHuNrdtBcIOz9twt2NzqqWKmYvccY8QQmoorgeYPskCRPFOa98tZ5HOxp5XYlo72dSlH6Po1GuOOxPviYFCtG3zkeQQa3NuT7butmuZ6+8fJ/o9Ht/af7bR5L/XAJ+Ix31LRHpE5Bmq+ni9AxZ9kClzOsOE4B8kfMuD8f3+ZP9xYHMvrE+d9gDBedhs7DY1pU03ZqXJbEbaJOAohSL3ERT6QMhp9zLVoGQqhZgWBisEeUsSp5rklZRQGAZz3unTYQeFeKdzKdjxE8CGfliVt3W1piiWKjHHvYw4aw5FE6XcaJhrsHTJ4SKnbYsVidRglKVJ0qfEdnXbMH85bhHpAO4BzgM+oqp3i8jbgOtF5D8BXwXeraqjwDaC3Bi9cVtrCXf6aJoWEisQA3FbH0UwmXkZBrYcha1HQcx5bCYI9gaKQmHTjaWO27qPpWOPDDLVTnuir6jreYyi7qdepWSaK2zHoF+KWCWlvYbatMkJyh12Wmee1kVuOA7rU+FeS22aJJnYt+ZRNM1xp2nBaDh0oMj0mVibnufD/FhzW/Mu+fAM7cxJpko2isju5P2NqnpjzXlVJ4BdItIDfEpELgTeQ5CUTuBGwqzv7wWkzq3VZdEcd1lFJYRAg0K8bV8bKMfG0RkEtvTB6QcJKruJUAg2UEzsm06/bvY4Ds1qY49wKJzsicmioj5tIpv2vD9KUVeUDtLjKZKlQVmuG2qHaFhGiIEOQigtK9kv7xU8Mgk9A9Cd5re74joVbfs1MGU18Tb1jWp8bLJ44Mw7Th6lVrQtXq2fZlmfg3Y3HifxnftV9ZKZ7KiqgyLyNeBKVf2TuHlURD4O/E583wvsSA7bTvCPdVn0GXDSlIm5muG4bYDpBcE66lhb1Q2TcHo/rO2HbqviTwtEJ0V0ZgVheLToqWnuxfoxlBUEMz+5c2l317KUSXtTpliMDjJ9RMG0MtO6xQ8Baydh7QB0D0BHOsl7WrOZOpmk148NGDVIEGKrJD9KKA+5YA/F7SbaeaVkbjralflIlYjIJmA8ivYq4BeB91veOrYieRXwQDzkDuAdInIroVLyaKP8NiyycKcOZ4ygsSfiZ0PUuoGygmC9jNcSjHbXKKwdDe1p61XUp2nDExQdFY5Q62LyQpD2kvQUydIld925Ay8zGMeSY1Phtt6KlilZS1H90jUKq0eh82ih2cso5oi0WJ2gEN/UeI9k69Rh5047z23nZqOdY3ieWpU8A7gp5rmXAbep6mdF5M4o6gLcB7w17v85QlPAvYR/2xubXWDRHXc98a7pKUxtrb4VBEtx9xCeK+wJdHU8Z1nT2GYFoawAWEWki3Z7UfYklY9Hk4p3Gqcm3FY5OETRiMRiNM2STBnueNE8W5LHbB6jaSebdCCssp6SniIpsKnL5vScqt8DnlOy/Yo6+yvw9tlcY9GFG6aLd95L2AIsHS8tHSvCeluuIBSMDkLL97KafhtQPm0imxeIPPDTJlQu2kuftIWJvTfyx2qLq25qR9yzOpnVhFhKG5CsJhiUlRRtwKFw3XaN9AdgjNpRCe1HoUys7X2abUl7SnrFZC1VLMMtIdxQXjE0SQjq/BHUJn7qouh+XJYuLEsb5h3Z0iUdcH40eZ0PgemivfRJ49G6whvj2X4QJ1WgNr/dTRFX5rDzRiQr4/Gd1D4dmrM3o5E2MkkdeNrUL30ytPjNu7eX5bXbOY69y/scUK9W3wI6dSDWacfSIVYQ7H3qYuzcUASwuZATyfs84FN37aLdfjQTbzMWULTSyM2BmQnrKm9Gw4xFmt+Gwmjk50rFO38SzNdpRWRZDNv9gsexDzI1R6SFZTJ7nzvvToruxx3ZGooCkQdpXiCsEDZyKOmPSRX/0c7J0Ui8U/KWGp0UYpo2cEo7Sc42VlNRziscx0rep80TXbTrU8WUUcsJNxQBlbtvC0LLA5rrTjtDLKPWweSUBbOJdv5Zes30vpz2YibinYq2iWc6UOAo0x12GrswveNPHpu2Th14er28DsZFuznuuOeB3H0vY3ohsqZZaQHIC4KRB20u0vWCPD3GaU/qiXdqJix+rOuACWr+NFhmMPK5L/OYTFtXla1naj7Sczue4543GrnvdJRXqF8I6p2zTMjT7ek2x6lXB5PvY0LemRwzTq1Yl3WXz+tjysYSafa0CG5AZksV/yYtL9xGPQHPP4PiF7Sei6n33gPcaUZuGHL3nYqvfZ6n8aCx206vZdcpe0qs98SYH5ufzymYj3bcC0FlhNtoFHxWUGa672zO7ThGIxNhn+cCnjpuqJ/OS5uvQvlTYb7NBfvkUao5yXflhDslD8hmj7HNjnec2VCWriv7PHXmtp8JuZE/LZaJd7p9Jmm+/DxOOe64FxkPUmehydN0ZcLZQa2IG6lwN6ogqyfG9cQ9/8ypj7cqcZw2pp6A18uBlx1T75yNtrlgnxou3I7jNHTgjdqAzwavn5k7vDmg4zhTpELaKA9u1OtzMNNrOCeP57gdx5lGWSV6jgv14uCpEsdxZkQVhWIpU8X/R70hPRzHcZY8luOezdIMEekSkW+LyP0i8qCI/GHc/lMicreI7BGR/y0inXH7yvh+b/z87GbXcOF2HKdtsVTJbJYZMApcoaoXAbuAK0XkMuD9wAdUdSfwJPCmuP+bgCdV9TzgA3G/hrhwO47T1kzOcmmGBo7HtzYopAJXALfH7TcRJgwGuCa+J37+kjihcF1cuB3HaVvmyXEjIh0ich/wBPBl4BFgUFWfjrv0Atvi623AAYD4+VHg9Ebn98pJx3Halkn44hBsnOVhXSKyO3l/o6remO6gqhPALhHpAT4F/EzJeTSuy9y1lmybwoXbcZy2RVWvnOfzD4rI14DLgB4RWR5d9XagL+7WC+wAekVkOXAaMNDovJ4qcRzHmUNEZFN02ojIKuAXgR8Cfw+8Ou72euDT8fUd8T3x8ztVtaHjlkafd4s0PNhZOIZVG1ZWtDseq61Du8eqiPwsobLRhl2/TVXfKyLnALcCG4DvAtep6qiIdAGfBJ5DcNrXquq+htdw4a4G7V4YmuGx2jp4rM4/nipxHMepGC7cjuM4FaNhqsRxHMdpPdxxO47jVAwXbsdxnIrhwu04jlMxXLgdx3Eqhgu34zhOxXDhdhzHqRgu3I7jOBXDhdtxHKdiuHA7juNUDBdux3GciuHC7TiOUzFcuB3HcSqGC7fjOE7FcOF2HMepGC7cJk+PHgAAAWVJREFUjuO0BSLy6yJyX1wmk9d/VrLvX4jI5dm248nrq0Vkj4icuRD3nuPjcTuO01aIyDbgm6p6VoN97gOeq6oTybbjqrpGRF4C3Ai8VFUfmf87ns7yxbio4zjOInIh8P16H4rIzwA/SkU7+ewFwEeBqxdLtMGF23Gc9uPZwAMNPr8K+ELJ9pXAp4EXq+pD83FjM8Vz3I7jtBsNHTfwMsqFexz4JvCm+bip2eDC7ThOu1HXcYvIaqBHVftKPp4E/iXwcyLy+/N4f03xVInjOG2DiCwDdgL1Uh2/APx9veNVdUREXg7cJSKHVPVj83CbTXHhdhynnTgP6FXV0TqfXwXc3ugEqjogIlcC/ygi/ar66bm+yWZ4c0DHcZyIiNwLPE9Vxxf7Xhrhwu04jlMxvHLScRynYrhwO47jVAwXbsdxnIrhwu04jlMxXLgdx3Eqhgu34zhOxXDhdhzHqRgu3I7jOBXj/wPsnBEqyqUR2wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# plate size, mm\n", "w = h = 10.\n", "# intervals in x-, y- directions, mm\n", "dx = dy = 0.1\n", "# Thermal diffusivity of steel, mm2.s-1\n", "D = 4.\n", "\n", "Tcool, Thot = 300, 700\n", "\n", "nx, ny = int(w/dx), int(h/dy)\n", "\n", "dx2, dy2 = dx*dx, dy*dy\n", "dt = dx2 * dy2 / (2 * D * (dx2 + dy2))\n", "\n", "u0 = Tcool * np.ones((nx, ny))\n", "u = u0.copy()\n", "\n", "# Initial conditions - circle of radius r centred at (cx,cy) (mm)\n", "r, cx, cy = 2, 5, 5\n", "r2 = r**2\n", "for i in range(nx):\n", " for j in range(ny):\n", " p2 = (i*dx-cx)**2 + (j*dy-cy)**2\n", " if p2 < r2:\n", " u0[i,j] = Thot\n", "\n", "def do_timestep(u0, u):\n", " # Propagate with forward-difference in time, central-difference in space\n", " u[1:-1, 1:-1] = u0[1:-1, 1:-1] + D * dt * (\n", " (u0[2:, 1:-1] - 2*u0[1:-1, 1:-1] + u0[:-2, 1:-1])/dx2\n", " + (u0[1:-1, 2:] - 2*u0[1:-1, 1:-1] + u0[1:-1, :-2])/dy2 )\n", "\n", " u0 = u.copy()\n", " return u0, u\n", "\n", "# Number of timesteps\n", "nsteps = 101\n", "# Output 4 figures at these timesteps\n", "mfig = [0, 10, 50, 100]\n", "fignum = 0\n", "fig = plt.figure()\n", "for m in range(nsteps):\n", " u0, u = do_timestep(u0, u)\n", " if m in mfig:\n", " fignum += 1\n", " print(m, fignum)\n", " ax = fig.add_subplot(220 + fignum)\n", " im = ax.imshow(u.copy(), cmap=plt.get_cmap('hot'), vmin=Tcool,vmax=Thot)\n", " ax.set_axis_off()\n", " ax.set_title('{:.1f} ms'.format(m*dt*1000))\n", "fig.subplots_adjust(right=0.85)\n", "cbar_ax = fig.add_axes([0.9, 0.15, 0.03, 0.7])\n", "cbar_ax.set_xlabel('$T$ / K', labelpad=20)\n", "fig.colorbar(im, cax=cbar_ax)\n", "plt.show()" ] } ], "metadata": { "colab": { "include_colab_link": true, "provenance": [] }, "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.9.12" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }