5. 1D diffusion equation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "59321108",
   "metadata": {},
   "source": [
    "Step 4: Diffusion Equation in 1-D\n",
    "-----\n",
    "***\n",
    "The one-dimensional diffusion equation is:\n",
    "\n",
    "$$\\frac{\\partial u}{\\partial t}= \\nu \\frac{\\partial^2 u}{\\partial x^2}$$\n",
    "\n",
    "The first thing you should notice is that —unlike the previous two simple equations we have studied— this equation has a second-order derivative. We first need to learn what to do with it!\n",
    "\n",
    "### Discretizing $\\frac{\\partial ^2 u}{\\partial x^2}$\n",
    "\n",
    "The second-order derivative can be represented geometrically as the line tangent to the curve given by the first derivative.  We will discretize the second-order derivative with a Central Difference scheme: a combination of Forward Difference and Backward Difference of the first derivative.  Consider the Taylor expansion of $u_{i+1}$ and $u_{i-1}$ around $u_i$:\n",
    "\n",
    "$u_{i+1} = u_i + \\Delta x \\frac{\\partial u}{\\partial x}\\bigg|_i + \\frac{\\Delta x^2}{2} \\frac{\\partial ^2 u}{\\partial x^2}\\bigg|_i + \\frac{\\Delta x^3}{3!} \\frac{\\partial ^3 u}{\\partial x^3}\\bigg|_i + O(\\Delta x^4)$\n",
    "\n",
    "$u_{i-1} = u_i - \\Delta x \\frac{\\partial u}{\\partial x}\\bigg|_i + \\frac{\\Delta x^2}{2} \\frac{\\partial ^2 u}{\\partial x^2}\\bigg|_i - \\frac{\\Delta x^3}{3!} \\frac{\\partial ^3 u}{\\partial x^3}\\bigg|_i + O(\\Delta x^4)$\n",
    "\n",
    "If we add these two expansions, you can see that the odd-numbered derivative terms will cancel each other out.  If we neglect any terms of $O(\\Delta x^4)$ or higher (and really, those are very small), then we can rearrange the sum of these two expansions to solve for our second-derivative. \n",
    "\n",
    "$u_{i+1} + u_{i-1} = 2u_i+\\Delta x^2 \\frac{\\partial ^2 u}{\\partial x^2}\\bigg|_i + O(\\Delta x^4)$\n",
    "\n",
    "Then rearrange to solve for $\\frac{\\partial ^2 u}{\\partial x^2}\\bigg|_i$ and the result is:\n",
    "\n",
    "$$\\frac{\\partial ^2 u}{\\partial x^2}=\\frac{u_{i+1}-2u_{i}+u_{i-1}}{\\Delta x^2} + O(\\Delta x^2)$$\n",
    "\n",
    "### Back to Step 4\n",
    "\n",
    "We can now write the discretized version of the diffusion equation in 1D:\n",
    "\n",
    "$$\\frac{u_{i}^{n+1}-u_{i}^{n}}{\\Delta t}=\\nu\\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\\Delta x^2}$$\n",
    "\n",
    "As before, we notice that once we have an initial condition, the only unknown is $u_{i}^{n+1}$, so we re-arrange the equation solving for our unknown:\n",
    "\n",
    "$$u_{i}^{n+1}=u_{i}^{n}+\\frac{\\nu\\Delta t}{\\Delta x^2}(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})$$\n",
    "\n",
    "The above discrete equation allows us to write a program to advance a solution in time. But we need an initial condition. Let's continue using our favorite: the hat function. So, at $t=0$, $u=2$ in the interval $0.5\\le x\\le 1$ and $u=1$ everywhere else. We are ready to go!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5c52a7f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e895d0fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f69b35e6860>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjqUlEQVR4nO3de3jU9Z0v8PdnJpnc7zOQEAi5kAiSgGIECRDQdlsutlptXZW11apUt+3TbU/7tLvnOXtp97S7Z7e77R7rhXqstbZqt9rWCtqbQoAAEhBIuCUkJBCSmMmF3G8z8z1/JIOAuUyS38x35jfv1/PwkGR+M/Nm8uOdyed3E6UUiIgo9Fl0ByAiImOw0ImITIKFTkRkEix0IiKTYKETEZlEhK4nttvtKjs7W9fTExGFpMOHD7cppRzj3aat0LOzs1FRUaHr6YmIQpKINEx0G0cuREQmwUInIjIJFjoRkUmw0ImITIKFTkRkEix0IiKTYKETEZkEC520qHP24rm953CiqUt3FCLT0HZgEYWfYZcHfzz5Pn5+sAHlte2Xv37DgmRsXZWF25fNQ4zNqjEhUWgTXRe4KC4uVjxSNDxc6OjHy4fO45VDjWjrHUJmcgzuW7kAGwvTsbu6Db842IBaZx8SoyNw903zsXVVFhbNSdAdmygoichhpVTxuLex0MlfDtV34Ml3zmJXtRMC4LbFc7B11UKUFjhgtcjl5ZRSOHiuAz8/eB5vVTVjxK2wMicVj63PxW2L5+r7BxAFIRY6BVzVxS7c/VQ5kmIice/KLNx78wLMS46Z8n5tvUP41eFG/PxgAy52DuAnD63E+oJxz0NEFJZY6BRQ7b1D+OQT+6CUwutfXgt7fNS0H6N/2IW7nixH06UB/O7La7EwLc4PSYlCz2SFzr1cyFAjbg+++IsjaOsdwjMPFM+ozAEg1haB7Q8UQ0Sw7YXD6BtyGZyUyHxY6GSo7+48hQN1HfjeXUUomp80q8fKSovFE/ffiJrWHnzjV8eg67dJolDBQifDvHq4ET/ZV4+H1mTjrhXzDXnMdfkOfHPjYuysbMGTu2oNeUwis2KhkyGON17C3/66Eqtz0/B3m5cY+tjbSnPxieXz8O9/OIN3zrQa+thEZsJCp1lz9gzhCz87DEd8FJ64/0ZEWo1drUQE/+fuZVicnoivvPQezrX1Gfr4RGbBQqdZ8W4E7ewfxjMP3IS0GW4EnUqMzYrtD9wEq0Ww7YUK9HIjKdGHsNBpVv75jZN491wH/vXuZSjMnN1G0KksSI3FE/evQK2zF1//JTeSEl2LhU4z9lZVM366vwGPrM3BHTdkBuQ51yyy4+82L8FbJ1rw0/L6gDwnUahgodOMKKXwgz/VYNGceHxr0+KAPvfDa3OwOjcNT+6qxZDLHdDnJgpmLHSakV3VTpxu6cEXSnMRYfBG0KmICP761jy09gzh10cuBvS5iYIZC51m5KldtchIig7YqOVaaxfZUZiZiGfK6uD2cJZOBLDQaQYON3Tg3XMdeGRdLmwRelYhEcHj6xfhXFsffn+iRUsGomDDQqdpe2pXHZJjI3HvzQu05thYmI7stFg8tauWe7wQgYVO01Tzfg/+dOp9fG51NuKi9F7wymoRfGF9Hiovdl11BSSicMVCp2l5encdYiKt+FxJtu4oAIC7VmRiTkIUnuJ5XohY6OS7i5cG8NujF3HvygVIjbPpjgMAiIqw4uG1Odh7tg3HGy/pjkOkFQudfPbsnjoAwCPrcjUnudr9q7KQEB2Bp3fzXTqFNxY6+aSjbxgvv3sBd9yQiUwfLiUXSAnRkfjs6oV4s6oFdc5e3XGItGGhk09+Wl6PgRE3HlsfXO/OvR4syUGk1YLtZXW6oxBpw0KnKfUNufDT/fX46JK5yJ+boDvOuBwJUbineD5ePdKIlq5B3XGItGCh05RePnQBl/pH8PiGPN1RJrVtXR7cHoXn9p3THYVICxY6TWrY5cGze+qwMicVNy1M0R1nUllpsbh92Tz8/EADuvpHdMchCjgWOk3qt0cvorlrMOjfnXs9tj4PfcNu/OxAve4oRAE3ZaGLyHMi0ioiVRPcniQivxORYyJyQkQeMj4m6aCUwjNldVicnoANBQ7dcXxy/bxEbLjOgZ/sq8fgCE+tS+HFl3fozwPYOMntXwRwUim1HMAGAN8XkeA46oRm5eiFSzjb2ovPr82BiOiO47PPr8lBe98wdvGC0hRmpix0pVQZgI7JFgGQIKP/4+PHluUFH01gx/FmRFoFH1+arjvKtJTkpSE1zoY3jjfrjkIUUEbM0J8AsARAE4BKAF9RSnnGW1BEtolIhYhUOJ1OA56a/MXjUdhZ2YzSfAeSYiJ1x5mWCKsFGwvT8edTrRgY5tiFwocRhf5xAEcBzANwA4AnRCRxvAWVUtuVUsVKqWKHIzRmsuHqaOMlNHUNYsuyDN1RZuT2ogwMjLg5dqGwYkShPwTgNTXqLIBzAAJ7kUky3I7jzbBZLfjo9XN1R5mRlTmpsMfb8EYlxy4UPowo9PMAPgIAIjIXwHUAePx1CLs8bilwIDE6tMYtXt6xy9unWtE/zE06FB582W3xJQD7AVwnIo0i8rCIPCYij40t8h0AJSJSCeDPAL6plGrzX2Tyt/cudKK5axC3h+i4xWtL0TwMjLjxzmlur6HwMOUlZ5RS901xexOAjxmWiLR743gzbBEWfGTJHN1RZmV07BKFHZVNIbstgGg6eKQoXcU7bllf4EBCiI5bvKwWwabCdLx9uhV9Qxy7kPmx0Okqh8934v3uoZAft3htWZaBwREP3j7NvV3I/FjodJUdl8ctobl3y7Vuzk6FIyEKO7m3C4UBFjpd5h233HqdA/FRU25eCQlWi2Azxy4UJljodFlFQydae4awZdk83VEMtWXZPAy5PPgzxy5kcix0umzH8SZERVjwkcWhvXfLtYoXpmBOQhR2HG/SHYXIr1joBABwexR2VrXgtsVzEGeScYuXxSLYXJSBd8440cuxC5kYC50AAIfqO+DsGTLt/tq3L8vAsMuDP596X3cUIr9hoRMAYGdlM6IjLbjNZOMWrxVZKUhPjMYOnlKXTIyFTqPjlsrRcUuszVzjFi+LRbCpKB27qp3oGeT1RsmcWOiEd891oK13CFuKzLV3y7U+GLtwbxcyJxY6YUdlE2Iirbh1sbnPUX/jghRkJEXzSkZkWiz0MOdye/BWVQtuW2LecYuXd2+Xsmonujl2IRNioYe50XHLMG4vMufeLdfasiwDw24P/nSSe7uQ+bDQw9zOqmbE2qzYcJ0592651o0LkpGZHMNzu5ApsdDDmFIKu844sWaRHTE2q+44ASEiuG3xHJTXtmPYNe61zIlCFgs9jJ1r60Nj5wBKC8y9MfRapQUO9A+7UdHQoTsKkaFY6GGsrHr00mzr88Or0FfnpSHCIiir5pUSyVxY6GGsrKYN2WmxyEqL1R0loOKjInDTwpTLP9CIzIKFHqaGXG7sr20Pu3GLV2mBAyebu+HsGdIdhcgwLPQwdbi+EwMjbpSG2bjFa/3YD7I9NXyXTubBQg9Tu2uciLQKVuel6Y6ixfUZiUiLs3HsQqbCQg9TZdVtuGlhiunOfe4ri0WwLt+OPTVt8HiU7jhEhmChh6HWnkGcau4O2/m5V2mBA+19wzjZ3K07CpEhWOhhaM/Y7nrhOj/3WptvBwDs5tiFTIKFHobKapxIi7Ph+oxE3VG0mpMQjSUZiZyjk2mw0MOMx6Owt6YN6/LtsFhEdxztSgvsOHK+k9caJVNgoYeZk83daO8bDvv5udf6fAdG3AoHatt1RyGaNRZ6mPHOi9eF+fzc66bsFMREWlHG/dHJBFjoYaas2onrMxLhSIjSHSUoREVYsTovjXN0MgUWehjpHXLhcEMnxy3XKM23o769H+fb+3VHIZoVFnoY2V/bDpdHobTArjtKUPH+gNvNsQuFOBZ6GCmrdiLWZkXxwlTdUYJKjj0O81NiOHahkDdloYvIcyLSKiJVkyyzQUSOisgJEdltbEQySlmNE6tz02CL4M/xK4kISgsc2F/bjhE3r2JEocuX/9nPA9g40Y0ikgzgSQCfVEotBfAZQ5KRoRra+9DQ3o91+Ry3jKc0347eIReONHTqjkI0Y1MWulKqDMBk1+q6H8BrSqnzY8u3GpSNDOQdJ3CD6PhKFtlhtQh3X6SQZsTv3gUAUkRkl4gcFpHPTrSgiGwTkQoRqXA6+R8nkHZXt2F+Sgxy7HG6owSlxOhI3LggmZelo5BmRKFHALgJwBYAHwfwv0SkYLwFlVLblVLFSqlih4PvFANl2OXB/to2lBY4IMLD/SdSWuBAVVMX2nt5FSMKTUYUeiOAt5RSfUqpNgBlAJYb8LhkkCPnO9E3HL5XJ/JVaYEDSgF7z/JdOoUmIwr9twDWiUiEiMQCWAXglAGPSwYpq3bCahGULArPqxP5qigzCcmxkTydLoWsKS9XIyIvAdgAwC4ijQD+AUAkACilnlZKnRKRtwAcB+AB8KxSasJdHCnw9tS0YUVWMhKjI3VHCWpWi2DtotGrGCmlOJ6ikDNloSul7vNhmX8D8G+GJCJDtfcOoaqpC1/76LibNegapQUOvHG8GadberAkzM8XT6GHR5iY3L7adigFrOPuij7xbmfYW8M5OoUeFrrJlZ9tQ0J0BIoyk3RHCQnpSdHIdcRhXy0LnUIPC93kymvbcUtuGqy8OpHP1uTZ8e65Dp4GgEIOC93ELnT043xHP0ryuHfLdJTkpaF/2I1jFy7pjkI0LSx0E9s/dlm1NYt4/pbpuCU3DSKjv90QhRIWuontq22DPT4K+XPidUcJKSlxNlyfkYh9PMCIQgwL3aSUUiivbUdJXhr3p56BNYvseO/8JQwMu3VHIfIZC92kzrb2wtkzhDU8OnRGSvLSMOz2oKJhshONEgUXFrpJeccFJXmcn8/EzdmpiLAI9p3lHJ1CBwvdpMpr27EgNQYLUmN1RwlJcVERuDErGfu5PzqFEBa6Cbk9Cgfq2lGSy3fns7E6z47Ki13oGhjRHYXIJyx0EzrR1IXuQRfPrjhLa/LS4FHAwTqOXSg0sNBNyDv3Xc0DimblhqxkREdauD86hQwWugmV17ahYG485iRE644S0qIirLg5OxXlnKNTiGChm8yQy41D9R3cu8UgaxbZUf1+L1p7BnVHIZoSC91k3jt/CYMjHp6/xSDe13E/xy4UAljoJlNe2w6LAKtyWehGWDovCYnRESjn/ugUAljoJlN+tg1FmUlIiuHl5oxgtQhuyU3j+dEpJLDQTaRvyIWjFy6hhGdXNNSaRXY0dg7gQke/7ihEk2Khm8i79R1weRTn5wbzvp48+yIFOxa6ieyvbYfNakHxwlTdUUxl0Zx4zEmI4v7oFPRY6Cay72wbVixMRozNqjuKqYgISvLSUF7bDqWU7jhEE2Khm0Rn3zBONndz/3M/Kcmzo613CNXv9+qOQjQhFrpJHKhrh1Lg+c/9xHteHB41SsGMhW4S+2rbEGezYtn8ZN1RTGl+SiyyUmN5fnQKaix0kyivbcfKnFREWvkt9Zc1i9JwsK4dLrdHdxSicfF/vwm0dA2iztmHNdz/3K9K8uzoGXKhqqlbdxSicbHQTcA71+Xpcv3L+/pyjk7BioVuAvvOtiMlNhJL0hN1RzE1e3wUFqcn8LwuFLRY6CFOKYXy2jaszkuDxSK645je6rw0HKrvwOCIW3cUog9hoYe4s629aO4axLp8h+4oYaE034EhlweH6jt0RyH6EBZ6iNtd7QQAlBaw0ANhVW4qbFYLysZed6JgMmWhi8hzItIqIlVTLHeziLhF5NPGxaOplNW0Ic8Rh8zkGN1RwkKsLQI356SgrJobRin4+PIO/XkAGydbQESsAP4VwO8NyEQ+Ghxx42BdO9+dB1hpvgNn3u9BSxcvS0fBZcpCV0qVAZhqYPhlAK8CaDUiFPnm3XMdGHJ5UMr5eUB5t1eU1XDsQsFl1jN0EckE8CkAT88+Dk3HnhonbFYLVuXydLmBtCQjAY6EKOyp4diFgosRG0V/AOCbSqkp9+MSkW0iUiEiFU4n393MVll1G27OSUGsLUJ3lLAiIliXb8feGifcHp5Ol4KHEYVeDOBlEakH8GkAT4rIneMtqJTarpQqVkoVOxwcE8xGS9cgzrzfw3GLJusLHOjsH0HVxS7dUYgum3WhK6VylFLZSqlsAL8C8NdKqd/M9nFpct75LTeI6rF27Lw53H2Rgokvuy2+BGA/gOtEpFFEHhaRx0TkMf/Ho4nsqWmDI2H0UHQKvLT4KBRmJnKOTkFlyuGrUuo+Xx9MKfXgrNKQT9wehb01Tty6eA5EeLi/LqX5Dmwvq0PP4AgSoiN1xyHikaKhqOpiFzr7R7Ce4xatSgsccHkULx5NQYOFHoLKqp0Q+WCOS3qsyEpBnM3KOToFDRZ6CCqrcaJwXhLS4qN0RwlrtggLVufZUVbjhFLcfZH0Y6GHmJ7BERw5fwnr8vnuPBiUFthxoWMADe39uqMQsdBDTXltO9wexd0Vg0QpTwNAQYSFHmLKqp2Is1mxIitFdxQCkG2PQ1ZqLOfoFBRY6CFEKYWyGidW59lhi+C3LliUFtixv7Ydwy6P7igU5tgKIaS+vR8XOgawvoDz82BSmu9A37Abhxs6dUehMMdCDyF7xua0vNxccFmdl4YIi1z+/hDpwkIPIWXVTmSlxiLbHqc7Cl0hIToSK7JSuGGUtGOhh4hhlwf7a9tRynFLUCotsKPqYjfaeod0R6EwxkIPEYcbOtE37ObpcoOUdzfSvTxZF2nEQg8Re2qciLAIVuel6Y5C4yicl4SU2EiOXUgrFnqIKKtxYkVWCs/qF6QsFsHafAf21LTxNACkDQs9BLT1DqHqYjfn50GuNN8OZ88QTjX36I5CYYqFHgK8c1ke7h/cvN+f3TxqlDRhoYeA359ogSMhCoXzknRHoUnMTYxGYWYi/nCyRXcUClMs9CDXN+TCO2dasakwHRYLr04U7DYVZuC985dw8dKA7igUhljoQe7t060YHPFgc1GG7ijkgy1j36c3K5s1J6FwxEIPcjsrm+FIiMLN2am6o5APsu1xWDovETtY6KQBCz2I9Q9/MG6xctwSMjYXjY5dmjh2oQBjoQcxjltCk3fsspPv0inAWOhBbGdlM+zxHLeEmmx7HK7PSGShU8Cx0INU/7ALb5/muCVUbVmWgSMcu1CAsdCDlHfcsmUZxy2haDPHLqQBCz1IcdwS2nI4diENWOhBiOMWc+DYhQKNhR6E3jnt5N4tJuD9/r1ZxVMBUGCw0IPQjsom2OOjsDKH45ZQlmOPw5KMROw43qQ7CoUJFnqQ8Y5bNhbO5bjFBLYUpXPsQgHDQg8y3nHLlqJ5uqOQATh2oUBioQcZ794tHLeYQ64jHku4twsFCAs9iHDcYk5bitJxuKETzV0cu5B/TVnoIvKciLSKSNUEt28VkeNjf8pFZLnxMcPDO6edGBhxc+8Wk/ngICOOXci/fHmH/jyAjZPcfg7AeqXUMgDfAbDdgFxhaXTcYsOqnDTdUchAuY54LE5P4NiF/G7KQldKlQHomOT2cqVU59inBwDMNyhbWBkYdo+NW3gwkRndviyDYxfyO6Nn6A8DeHOiG0Vkm4hUiEiF08kL6V7pnTOtHLeY2OW9XTh2IT8yrNBF5FaMFvo3J1pGKbVdKVWslCp2OHgF+yvt4LjF1Dh2oUAwpNBFZBmAZwHcoZRqN+Ixw0nvkAtvn2rFx5dy3GJmW4oyUNHQicbOft1RyKRmXegikgXgNQAPKKWqZx8p/Lxy6AIGRty4p3iB7ijkR3fdNB9Wi+CF/Q26o5BJ+bLb4ksA9gO4TkQaReRhEXlMRB4bW+TvAaQBeFJEjopIhR/zmo7L7cFze89hZU4qli9I1h2H/CgzOQZbijLwi4Pn0T04ojsOmVDEVAsope6b4vZHADxiWKIws6OyGRcvDeDbdyzVHYUCYFtpLl4/1oSX3z2PbaV5uuOQyfBIUY2UUvjxnjrkOeJw63VzdMehACjMTEJJXhqe21uPYZdHdxwyGRa6Rvvr2lF1sRuPrsuFhRtDw8ajpblo6R7EjkqeVpeMxULX6MdldbDH23DnjZm6o1AAbShwIH9OPLaXnYNSSnccMhEWuibV7/fgnTNOfG51NqIjrbrjUACJCB4tzcWp5m7sO8u9fMk4LHRNnt1Th+hIC/7qloW6o5AGd9wwD46EKGzfU6c7CpkIC12D1u5B/Oa9JtxTvAApcTbdcUiDqAgrHizJRlm1E6eau3XHIZNgoWvwfHk9XB4PHl6bozsKabR1VRZibVb8mO/SySAs9ADrG3LhxQMN2FiYjoVpcbrjkEbJsTbcU7wArx9t4lkYyRAs9AD7ZcUFdA+68Oi6XN1RKAg8vDYHHqXwfHm97ihkAiz0AHK5Pfh/e8/h5uwU3JiVojsOBYEFqbHYVJSBXxw4jx6eDoBmiYUeQG9WtaCxc4Dvzukq29blomfIhVcOXdAdhUIcCz1AlFLYXlaHHHscPrpkru44FESWL0jGypxUPLf3HEbcPB0AzRwLPUAOnutA5cUuPLIuh4f504d8oTQXTV2DvAAGzQoLPQDcHoXv/+EM0uJsuHsFL7lKH3brdXOQ54jDD/9cg8ERt+44FKJY6AGwvawOh+o78a1Ni3mYP43LYhH84yeXos7Zh39587TuOBSiWOh+VnWxC//xxzPYVJiOT9/Ed+c0sXX5DjxYko3ny+uxu5oXUafpY6H70eCIG3/zylGkxNrw3U8VQYSzc5rctzYtRv6ceHzjv4+hs29YdxwKMSx0P/qXN0/jbGsv/v0zy3nOFvJJdKQVP7j3BnT2D+NvX6vk6XVpWljoflJW7cTz5fV4sCQbpQUO3XEohCydl4T/8bHr8NaJFvzqcKPuOBRCWOh+0Nk3jK//9zHkz4nHtzYt1h2HQtCj63KxMicV//S7k7jQ0a87DoUIFrrBlFL4u19XorN/GP/5lzdwrxaaEatF8B/3LIcA+OorR+H2cPRCU2OhG+zVIxfxZlULvvYX16EwM0l3HAph81Ni8e07l6KioRNP767VHYdCAAvdQBc6+vGPr5/AyuxUbCvl+Vpo9u68IRNblmXgP/9YjcrGLt1xKMix0A3i9ih87ZdHIQC+f89yWHl4PxlARPC/7yyEPT4KX3nlPQwM8yhSmhgL3QDdgyN49IUKHKrvxD/dsRQLUmN1RyITSY614fv3LEedsw9bnz2A1p5B3ZEoSLHQZ6nO2Ys7f7QPZdVOfOfOQtzFc7WQH6xZZMeTW1fgVHMPPvl/9+HYhUu6I1EQYqHPwq4zrbjjR/twqX8ELz6yCg/cslB3JDKxzUUZePXxElgtgs88sx+/fo/7qNPVWOgzoJTCM7tr8fnnD2F+Sixe/9Ia3JKbpjsWhYHr5yXi9S+twYqsZHz1lWP47s5T3KWRLmOhT9PgiBtffeUovvfmaWwqzMCrj6/G/BTOzClw0uKj8LOHV+Gzqxdie1kdHnr+ELr6efk6YqFPS3PXAD7z9H785mgTvv6xAjxx/42ItUXojkVhKNJqwbfvKMT37irC/to23PnkPpxt7dEdizRjG/ngZFM3XjzYgN++dxEigh9/thh/cT0vI0f63bcyC4vmxOPxFw9j83/txZaiDPzVLVlYkZXCs3uGIdF1Nrfi4mJVUVGh5bl9MTjixo7jzXjxYAPeO38JUREWfGL5PDy+IQ95jnjd8Yiu0tw1gKd21eK1IxfRO+TC4vQEbL1lIe68YR4SoiN1xyMDichhpVTxuLex0K9W5+zFLw6ex6+ONOJS/whyHXHYumoh7l6RieRYngKXglvfkAuvH2vCiwcacKKpG3E2K+64MRNbV2Vh6TyeisIMZlXoIvIcgNsBtCqlCse5XQD8EMBmAP0AHlRKHZkqlO5C93gUGjsHcKqlG6eau3G6uQenW7pR396PCIvg40vTsfWWLKzOTeOvrhRylFI41tiFFw804HfHmjDk8iAzOQZLMhKwOD0Ri8f+zrHH8ajmEDPbQi8F0AvghQkKfTOAL2O00FcB+KFSatVUoYwudKUUXB6FIZcHfUMuXOofwaX+YVwaGEFX/wg6xz7u6B1GTWsPzrT0oG/sMGoRIDstDovTE3DDgmR8akUm5iREG5aNSKeu/hH85uhFHG7oxOmWbtQ6+y7v6hgVYUHB3ATkz42HIz4KSbGRSIqJRHKMDcljHyfFRCIhOgKRVsvYH+GbHI0mK/QpN4oqpcpEJHuSRe7AaNkrAAdEJFlEMpRSzTOLO7l3zrTin984iSGXB0MuD4ZdHgy53BhyeTDV9CjCIkiOtSHXEYdP3zQfizMSsSQjEQVz47m3CplWUmwkPleSjc+VZAMY3T50trUXp1t6cLq5G6dbelB+th0d/cMYdnl8eswIi1wud1uEBRYRiAAWkcsfez8X4PIPgMs/BuSDv2b7wyEUf7T85c0L8Mg640/gZ0SLZQK4cMXnjWNf+1Chi8g2ANsAICsra0ZPlhQTicXpiYiKsCAq0gKb1YKoSCuiIrwfWxBri0BK7AfvMJJjI5Eca0Oczcp3FhT2oiOtKMxMGvf0zoMjblzqH0HXwNW/4fYOueDyeDDiVhh2eTDi9v5RGHZ7oJSCxwMoKHgU4FEKSo3+5uw97sn7fss7FVBXfnGG1GwfQBN7fJRfHteIQh+vIcd9lZVS2wFsB0ZHLjN5shVZKVixNWUmdyWiKURHWpGeZEV6EkeOociIA4saASy44vP5AJoMeFwiIpoGIwr9dQCflVG3AOjy1/yciIgmNuXIRUReArABgF1EGgH8A4BIAFBKPQ1gJ0b3cDmL0d0WH/JXWCIimpgve7ncN8XtCsAXDUtEREQzwpNzERGZBAudiMgkWOhERCbBQiciMgltZ1sUESeAhhne3Q6gzcA4RgnWXEDwZmOu6WGu6TFjroVKKcd4N2gr9NkQkYqJTk6jU7DmAoI3G3NND3NNT7jl4siFiMgkWOhERCYRqoW+XXeACQRrLiB4szHX9DDX9IRVrpCcoRMR0YeF6jt0IiK6BgudiMgkgq7QRWSjiJwRkbMi8q1xbhcR+a+x24+LyApf7+vnXFvH8hwXkXIRWX7FbfUiUikiR0XE0Ctj+5Brg4h0jT33URH5e1/v6+dc37giU5WIuEUkdew2f75ez4lIq4hUTXC7rvVrqly61q+pculav6bKFfD1S0QWiMg7InJKRE6IyFfGWca/65dSKmj+ALACqAWQC8AG4BiA669ZZjOANzF6paRbABz09b5+zlUCIGXs403eXGOf1wOwa3q9NgB4Yyb39Weua5b/BIC3/f16jT12KYAVAKomuD3g65ePuQK+fvmYK+Drly+5dKxfADIArBj7OAFAdaD7K9jeoa8EcFYpVaeUGgbwMkYvQn2lyxelVkodAJAsIhk+3tdvuZRS5UqpzrFPD2D0yk3+Npt/s9bX6xr3AXjJoOeelFKqDEDHJIvoWL+mzKVp/fLl9ZqI1tfrGgFZv5RSzUqpI2Mf9wA4hdHrK1/Jr+tXsBX6RBec9mUZX+7rz1xXehijP4W9FIA/iMhhGb1QtlF8zbVaRI6JyJsisnSa9/VnLohILICNAF694sv+er18oWP9mq5ArV++CvT65TNd65eIZAO4EcDBa27y6/plxEWijeTLBacnWsbni1XPgM+PLSK3YvQ/3NorvrxGKdUkInMA/FFETo+9wwhEriMYPfdDr4hsBvAbAPk+3tefubw+AWCfUurKd1v+er18oWP98lmA1y9f6Fi/piPg65eIxGP0B8jfKKW6r715nLsYtn4F2zt0Xy44PdEy/rxYtU+PLSLLADwL4A6lVLv360qpprG/WwH8GqO/XgUkl1KqWynVO/bxTgCRImL35b7+zHWFe3HNr8N+fL18oWP98omG9WtKmtav6Qjo+iUikRgt858rpV4bZxH/rl9GbxiYzR+M/sZQByAHH2wYWHrNMltw9UaFd329r59zZWH0uqol13w9DkDCFR+XA9gYwFzp+OAAspUAzo+9dlpfr7HlkjA6B40LxOt1xXNkY+KNfAFfv3zMFfD1y8dcAV+/fMmlY/0a+3e/AOAHkyzj1/XLsBfXwG/SZoxuHa4F8D/HvvYYgMeueNF+NHZ7JYDiye4bwFzPAugEcHTsT8XY13PHvjnHAJzQkOtLY897DKMb00omu2+gco19/iCAl6+5n79fr5cANAMYwei7ooeDZP2aKpeu9WuqXLrWr0lz6Vi/MDoGUwCOX/F92hzI9YuH/hMRmUSwzdCJiGiGWOhERCbBQiciMgkWOhGRSbDQiYhMgoVORGQSLHQiIpP4/+dyBAJJD3DNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nx = 41\n",
    "dx = 2 / (nx - 1)\n",
    "nt = 20 # the number of timesteps\n",
    "nu = 0.3 # diffusion coefficient\n",
    "sigma = 0.2 # sigma is a numerical parameter, we will learn it later\n",
    "dt = sigma * dx**2 / nu  # dt is defined using sigma \n",
    "\n",
    "u = np.ones(nx) # a numpy array with nx elements all equal to 1\n",
    "u[int(0.5/dx):int(1/dx+1)] = 2 # setting u = 2 between 0.5 and 1 \n",
    "\n",
    "un = np.ones(nx) # our placeholder array to advance the solution in time\n",
    "\n",
    "for n in range(nt): # iterate through 0 to nt-1, in total nt steps\n",
    "    un = u.copy() # copy the existing values of u into un\n",
    "    for i in range(1, nx-1): # skip first and last component!\n",
    "        u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])\n",
    "        \n",
    "plt.plot(np.linspace(0,2,nx),u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66a743b1",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}