{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Step 1: 1D Linear Convection\n", "\n", "The first step in this journey begins with a simple equation called the 1-D Linear Convection equation. It is the most basic model that we can use to learn something useful about CFD.\n", "Here it is:\n", "\n", "$$ \\frac{\\partial u}{\\partial t} + c \\frac{\\partial u}{\\partial x} = 0 $$\n", "\n", "This equation can be uderstood as a representation of a *wave* that propagates whatever initial shape it had at a constant speed `c` without changing it. Assuming initial condutions to be $ u(x,0) = u_0(x) $ we can find an exact solution to this equation of the type $ u(x,t) = u_0(x-ct) $. Which represents the behavior of a wave being propagated with no change in shape. \n", "\n", "But we must discretize the above equation before we can work with it on our computer. To do so we perform a discretization in both space and time using the Forward Difference scheme for the tiem derivative and the backward difference scheme for the spatial derivative. \n", "\n", "This means that from now on we will consider space to have been discretized into **x** points that we index with i= 0 to **N** and step time in discrete intervals of size $ \\Delta t$.\n", "\n", "Performing a taylor series approximation with a forward difference in time and eliminating the higher order terms we can obtain a discretization for the time derivative:\n", "\n", "$$ \\frac{\\partial u}{\\partial t} \\approx \\frac{u_i^{n+1} - u^n_i}{\\Delta t} $$\n", "\n", "And similarly for the spatial with a backwards difference:\n", "\n", "$$ \\frac{\\partial u}{\\partial x} \\approx \\frac{u_{i}^{n} - u^n_{i-1}}{\\Delta t} $$\n", "\n", "Combining and substituting into our equation:\n", "\n", "$$\\frac{u_i^{n+1} - u^n_i}{\\Delta t} + c \\frac{u_{i}^{n} - u^n_{i-1}}{\\Delta t} = 0$$\n", "\n", "Our initial condition information will be enough to solve for everything in the above equation except $u_i^{n+1}$ so we can rearrange to solve as such:\n", "\n", "$$ u_i^{n+1} = u^n_i - c \\frac{\\Delta t}{\\Delta x}(u_i^n - u^n_{i-1}) $$\n", "\n", "Now that we have fully defined the equation, discretized it and found an expression that we can easily solve for we are going to implement it in python. \n", "\n", "### Libraries" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Adding inline command to make plots appear under comments\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import time, sys\n", "%matplotlib inline " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Key Variables\n", "\n", "Before beginning to calculate we must first define the grid of points that will characterize our spatial domain. We begin by giving it a `grid_length` parameter that defines the length of our spatial domain. Next, `grid_points` will define the number of points in the grid and `dx` will be the resulting distance between any pair of adjacent grid points. `nt` will represent the number of timesteps we want to calculate and `dt` represents the amount of time each timestep will cover. Finally `c` will represent the wave speed of our system. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "grid_length = 10\n", "grid_points = 41\n", "dx = grid_length / (grid_points - 1) \n", "nt = 400\n", "dt = 0.025\n", "c = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initial Conditions\n", "\n", "Next we set up our initial conditions for the wave. $ u = 1 $ everywhere except:\n", "\n", "$$ u = 2 \\Rightarrow 0.5 <= x <= 1 $$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 1. 2. 2. 2. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n", " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n" ] } ], "source": [ "u = np.ones(grid_points)\n", "u[int(.5/ dx):int(1 / dx + 1)] = 2\n", "print(u)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting initial conditions\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAHGpJREFUeJzt3X2UZHV95/H3p58bBgGdkeV58ARNMMkEHDckuskoZgWDkqwkEY0YDId9UFezObs+rEo2MYkm6olZEzhzyGQ06hiDGAjHB1wTnWRdPA5EeUwIaoDhaUZQhKGqu7rqu3/ce2vKsru6pvveur+e/rzO6cNM1e17v1XMqU/9ft97708RgZmZGcBY3QWYmVk6HApmZtblUDAzsy6HgpmZdTkUzMysy6FgZmZdDgVLgqS3Sbqq7joOF5I+I+k1dddha49DwZYk6fWS9kiak7Sz77ltkjqSnsh/9kr6hKTnDtjfZkkhaaL/uYj4vYi4tIKXsSKSXixpt6THJe2X9CVJL6u7rsVI+i1JH+l9LCLOi4gPVXCsnZLetcLffaqkT0k6IOkeSa8suz5bPYeCDfIA8C5gx1LPR8QG4CjgbOCfgL+XdM6I6ls1SeOLPHYh8FfAh4GTgOOAdwIvHW11h50/AebJ3s9XAVdIena9JdkPiAj/+GfgD1kw7Ox7bBuwd5FtPwjsWWI/m4EAJhZ57reAj/Rt9xrgXuDbwP/s2XYMeAvwDeAR4BPAU3ue/yvgIeAxYDfw7J7ndgJXAJ8GDgAv6qtD+TH/+4D3Ywx4O3APsI8sPI5ernbgBKDRV+uZ+TaT+d9fC9wJfAf4HHBqz7bPBj4PPAo8DLwNOJfsg7YFPAF8Pd/2i8Clq6l3kdd9WX6c+fxYf3MI/4aOzH/vmT2P/QXw7rr/ffvn+388UrCyXQOcJenIEvb1fOBZwDnAOyX9SP74G4BfAH6W7IP2O2TfQgufAU4Hng7cDHy0b7+vBH6XbITzD33PPQs4Gbh6QF2/lv+8AHgGsIEsDAfWHhEPAP8PeHlfLVdHREvSBWQf9P8B2AT8PbALQNJRwP8BPpu/5h8CvhARnwV+D/jLiNgQEVvKqrd/JxGxney9/IP8WC/Na7te0neX+Lk+//VnAgsRcVfPLr9OFnSWEIeCle0Bsm/bx5Swr/8VEY2I+DrZB0jxgfefyL7N7o2IObJRxoVFryIidkTE4z3PbZF0dM9+r42I/xsRnYho9h3zafl/HxxQ16uA90fENyPiCeCtwCv6eiVL1f4x4CIASQJekT9WvK7fj4g7I2KB7MP+JySdCpwPPBQR74uIZv76vjLozSup3mVFxPkRccwSP+fnm20Avtf3q4+RBbMlxKFgZTuRbDriuyXs66GePz9J9sECcCrwqeLbKNl0Sxs4TtK4pHdL+oak7wH/mv/Oxp593TfgmI/k/z1+wDYnkE3FFO4BJsjmyper/ZPAT0k6HvgZoEM2Iihe1wd6XtejZAF7Itno5RsDahpkNfWW5QngKX2PPQV4vOTj2Co5FKxsvwjcHBEHKjzGfcB5fd9IZyLifrLpmAuAFwFHk82ZQ/bhWhh0a+B/zvf/8gHbPED2AV44BVggm+cfKCK+A9wA/Epe68cjoqjnPuA/9r2u2Yj4cv7cM5ba7TKHXXG9wxwrP/31iSV+PpNvdhcwIen0nl/dAty+ghqsQg4FW5KkCUkzwDgwLmlmsdNJlTlR0uXApWTz4oNM5/sqfg713+GVwO/m0ypI2pTPx0M2HTFH9o3/CLIpmKHlH9D/DXiHpEskPUXSmKTnS9qeb7YL+A1Jp0nawME5/YUhD/Mx4GLgQg5OHRWv663FGTmSjpb0S/lz1wPHS3qTpGlJR0n6yfy5h4HNA97H1dbb62H6wimy0183LPFzXr7NAbJ+029LOlLS88jC+y9WUINVyKFgg7yd7GyZtwC/mv/57T3PnyDpCbKpga8CPwZsi4gbltnvE/m+ip8XHmJdHwCuA26Q9DhwI1B8QH6YbHrkfuCO/LlDEhFXk32Tfy3Zt+yHyc7AujbfZAfZh9lu4FtAk6z5PazryBrhD+Vz+MVxPwW8B/h4PvV1G1B8qD4O/BzZabEPAf9C1jiG7GwrgEck3bzI8VZbb68/A87Ip7j++hB/978As2RnQO0C/nNEeKSQGB0cuZqZ2XrnkYKZmXU5FMzMrMuhYGZmXQ4FMzPr+oHTC1O3cePG2Lx5c91lmJmtKTfddNO3I2LTctutuVDYvHkze/bsqbsMM7M1RdI9y2/l6SMzM+vhUDAzsy6HgpmZdTkUzMysy6FgZmZdlYWCpJMl/Z2kOyTdLumNi2wjSX8s6W5Jt0g6q6p6zMxseVWekroA/GZE3JwvJXiTpM9HxB0925xHdrfI08nucnkFB+92aWZmI1ZZKETEg+RLGkbE45LuJFtBqjcULgA+nN/D/kZJx0g6Pv/dZEQEH7nxHvY/PreyHUi8bMvx/NDTvfKgmaVtJBevSdoMnAn0ryl7It+/NOLe/LHvCwVJlwGXAZxyyilVlbmkBx5r8o5rb89rOfTfj4CHH2vyngt/vOTKzMzKVXko5Cs9fRJ4U0T0L9w9lIjYDmwH2Lp168gXgHhyLlug6n9fdCYv3XLCIf/+C977RZ5stcsuy8ysdJWefSRpkiwQPhoR1yyyyf1kC5IXTsofS0oj/0CfnRxf0e/PTI7TmHcomFn6qjz7SGRL990ZEe9fYrPrgIvzs5DOBh5LrZ8AdD/QZ6dWFgqzk2M0PVIwszWgyumj5wGvBm6V9LX8sbcBpwBExJXAp4GXAHcDTwKXVFjPihUjhZkVjhRmp8a7+zAzS1mVZx/9AzCwLZufdfS6qmooS3OV00ezk+N850CrzJLMzCrhK5qH0O0prHD6aGZynOaCRwpmlj6HwhCarQ6wupFC041mM1sDHApDKBrNM5Mre7tmJt1TMLO1waEwBDeazWy9cCgModlqI8H0xMpHCs1Wh05n5NfdmZkdEofCEBrzbWYnx9FK7nHBwV7E3EKnzLLMzErnUBhCo9VecZMZsovXiv2YmaXMoTCEZquz4n4CHDyV1Vc1m1nqHApDaLbaK75GAQ42qD1SMLPUORSG0Gi1V3w6KvSEgq9VMLPEORSGUDSaV6r4XU8fmVnqHApDyEYKq+8pePrIzFLnUBhCc9VnH3n6yMzWBofCEBpuNJvZOuFQGMKqRwp5oMy1fPGamaXNoTCExvwqewoeKZjZGuFQGMJqL16b8RXNZrZGOBSWsdDuMN/urGr6aGbCjWYzWxscCsto5jexm51a+Vs1NiamJ8Z8nYKZJc+hsIzi2/1qRgrgNRXMbG1wKCyjucoFdgqzk+MeKZhZ8hwKyyg+yFdznQJkodDwKalmljiHwjKKKZ/VTh/NTI670WxmyXMoLKP4IF/t9NHMpBvNZpY+h8IyGmX1FNxoNrM1wKGwjGZJ00eznj4yszXAobCMRkmN5hmffWRma4BDYRmN+fziNZ+SambrgENhGaVNH7mnYGZrgENhGd1G8ypucwHFdQoOBTNLm0NhGc1WGwmmxlf3Vk1PjtNsdeh0oqTKzMzKV1koSNohaZ+k25Z4/mhJfyPp65Jul3RJVbWsRmM+W2BH0qr2U0w/zS34qmYzS1eVI4WdwLkDnn8dcEdEbAG2Ae+TNFVhPSvSWOWqa4VZr6lgZmtAZaEQEbuBRwdtAhyl7Cv4hnzbharqWalGa3WrrhWKU1odCmaWsjp7Ch8EfgR4ALgVeGNELDq3IukySXsk7dm/f/8oa8zWZ17lNQpw8Ipon5ZqZimrMxReDHwNOAH4CeCDkp6y2IYRsT0itkbE1k2bNo2yRpqt1a26Vuiu0+yrms0sYXWGwiXANZG5G/gW8MM11rOoxny7u8byanikYGZrQZ2hcC9wDoCk44BnAd+ssZ5FuadgZuvJRFU7lrSL7KyijZL2ApcDkwARcSXwO8BOSbcCAt4cEd+uqp6VarbaPP2o6VXvx9NHZrYWVBYKEXHRMs8/APz7qo5flkbJjWaPFMwsZb6ieRnFxWurVQTLnJfkNLOEORSWUVpPwSMFM1sDHArLmGt1Spk+ciiY2VrgUBhgod1hvt1hZmL1oTA9kd/mwo1mM0uYQ2GAZn7zutlV3jYbYGxMTE+M+ToFM0uaQ2GA4lt9GY1m8EI7ZpY+h8IAxbf6MhrNkC+04+kjM0uYQ2GA4lt9GY1myNdp9noKZpYwh8IAZU8fzXikYGaJcygMUEwfldlTcKPZzFLmUBigmD6aLm2kMOZGs5klzaEwQOkjBU8fmVniHAoDlN1onpn09JGZpc2hMEBjPr94rcyRgkPBzBLmUBig4Uazma0zDoUBuhevlXCbC/BIwczS51AYoNlqMyaYGi/nbcp6Ch06nShlf2ZmZXMoDNCYz9ZSkFTK/orbZcz5qmYzS5RDYYBGq5xV1wqzk2Pd/ZqZpcihMEBZq64VilNbHQpmliqHwgDNVru0axTg4PSRz0Ays1Q5FAZozJc9fTTe3a+ZWYocCgOU3lOY8kjBzNLmUBig2eowU+L0UXek4FAws0Q5FAZottrMTJT3Fs14+sjMEudQGKBRUaPZIwUzS5VDYYDSG83uKZhZ4hwKA5R+nUL3lFRf0WxmaXIoDFD2dQpuNJtZ6hwKS1hod2i1o9Tpo+m8ae1Gs5mlqrJQkLRD0j5Jtw3YZpukr0m6XdKXqqplJZoL5S6wAzA2JmYmx9xTMLNkVTlS2Amcu9STko4B/hR4WUQ8G/ilCms5ZMW3+ZnJct+iGa+pYGYJqywUImI38OiATV4JXBMR9+bb76uqlpXoLrBT4kgB8oV2PH1kZomqs6fwTOBYSV+UdJOki5faUNJlkvZI2rN///6RFNddirPERjN49TUzS1udoTABPAf4eeDFwDskPXOxDSNie0RsjYitmzZtGklxxbf5MnsKcHD1NTOzFE3UeOy9wCMRcQA4IGk3sAW4q8aaurojhbKnj6bG3Wg2s2TVOVK4Fni+pAlJRwA/CdxZYz3fpwiFMm+IB54+MrO0VTZSkLQL2AZslLQXuByYBIiIKyPiTkmfBW4BOsBVEbHk6aujNlfRSGFmcpxHD8yXuk8zs7JUFgoRcdEQ2/wh8IdV1bAajYrOPvJ1CmaWMl/RvITGfPkXrxX78/SRmaXKobAEN5rNbD1yKCyhe/HaVLlvkUcKZpYyh8ISGvNtxgRT4+Xf5qLZ6tDpRKn7NTMrg0NhCc1WtsCOpFL3W1whPbfgC9jMLD0OhSWUvRRnwWsqmFnKHApLaLTaTE+UHwrFXVcdCmaWIofCEspeda1QXPfgO6WaWYocCktozLdLPx0VetdpdiiYWXocCktotCoKhSmHgpmly6GwhEarU/rN8MCNZjNLm0NhCXOtNrMlL8UJ7imYWdocCkuoevrIIwUzS9FQd0mV9M7FHo+I3y63nHQ05tul3yEVDo4U3FMwsxQNe+vsAz1/ngHOJ6EFcarQaFUTCrOePjKzhA0VChHxvt6/S3ov8LlKKkpEVdcpdE9J9W0uzCxBK+0pHAGcVGYhKWm1O7TaUUlPYXoiv6LZIwUzS9CwPYVbgeK2nuPAJuCw7Sc0K1pLAWBsTF59zcySNWxP4fyePy8AD0fEQgX1JKHZyqZ2qrhOAbymgpmla9iewj1VF5KSKkcKxX49fWRmKfJ1CosovsXPVHDxWrZfjxTMLE0OhUUU3+KrGikUq6+ZmaXGobCIRtXTR1PjbjSbWZIcCovoTh+50Wxm64xDYRHNEUwfudFsZilyKCyiueDpIzNbnxwKi2jM59cpVDVSmBjz9JGZJcmhsIhRNJodCmaWIofCIprdRnM1b8/spKePzCxNDoVFNObbjAmmxqu7eK3Z6tDpxPIbm5mNkENhEcWqa5Iq2X9xS+453z7bzBJTWShI2iFpn6TbltnuuZIWJF1YVS2HqlHRWgqF7kI7nkIys8RUOVLYCZw7aANJ48B7gBsqrOOQNStada3gUDCzVFUWChGxG3h0mc3eAHwS2FdVHStRdShMT3qhHTNLU209BUknAr8IXDHEtpdJ2iNpz/79+yuvrTHfrux0VOhZktMjBTNLTJ2N5j8C3hwRy3ZbI2J7RGyNiK2bNm2qvLCi0VyVol/hUDCz1Ay78loVtgIfz8/w2Qi8RNJCRPx1jTUB0Gh1OHp2srL9u6dgZqmqLRQi4rTiz5J2AtenEAiQ3RDv3zxlurL9F/0K9xTMLDWVhYKkXcA2YKOkvcDlwCRARFxZ1XHLMKrpI48UzCw1lYVCRFx0CNv+WlV1rERzRNcpuKdgZqnxFc2LaLTaTE9UFwqePjKzVDkUFjGykYJvc2FmiXEo9Gm1O7TaUWlPYXrCF6+ZWZocCn2aFa+lADA2JmYmx9xTMLPkOBT6NLprKVQXCpCFjs8+MrPUOBT6NPOlOKscKRT79/SRmaXGodCnuVD99BFkIxGPFMwsNQ6FPsW395nJat+amQkvyWlm6XEo9GmMoNEM2VXNzZZPSTWztDgU+rjRbGbrmUOhT3N+RD0FN5rNLEEOhT6jnT5yKJhZWhwKfbqhUPn00Zinj8wsOQ6FPkXzt8o1msE9BTNLk0OhTzGlU/kpqe4pmFmCHAp9GvNtxgRT49WHwtxCh04nKj2OmdmhcCj0KVZdy9eOrkzRs5jz7bPNLCEOhT6NitdSKBRnN7mvYGYpcSj0ac63K28yg0PBzNLkUOhTTB9Vrbhi2s1mM0uJQ6FP1UtxFrpLcnqkYGYJcSj0abTazEyMYKSQn/LqUDCzlDgU+jRancpvhgfuKZhZmhwKfZrzbWYrvnANDl4x7Z6CmaXEodBnVI3mom/hkYKZpcSh0GfU1ym4p2BmKXEo9Bn5dQqePjKzhDgU+jQXRj195NtcmFk6HAo9Wu0OrXaMZKQwPeFTUs0sPQ6FHs0RrboGIImZyTGHgpklpbJQkLRD0j5Jty3x/Ksk3SLpVklflrSlqlqGVZwJNIrrFMAL7ZhZeqocKewEzh3w/LeAn42IHwN+B9heYS1Dac5n8/ujGCkUx3Gj2cxSMlHVjiNit6TNA57/cs9fbwROqqqWYTVGOH0E2YjEIwUzS0kqPYVfBz6z1JOSLpO0R9Ke/fv3V1ZEt6cwNZq3ZXZy3D0FM0tK7aEg6QVkofDmpbaJiO0RsTUitm7atKmyWro9hVFOHzkUzCwhlU0fDUPSjwNXAedFxCN11gKjD4UZh4KZJaa2kYKkU4BrgFdHxF111dGrOT/inoIbzWaWmMpGCpJ2AduAjZL2ApcDkwARcSXwTuBpwJ9KAliIiK1V1TOMUTeaZ6fcUzCztFR59tFFyzx/KXBpVcdfiW4ojOw6hTFPH5lZUmpvNKekmMpxo9nM1iuHQo+5hdFevDYz5Z6CmaXFodCjMd9mfExMjmskx5udHGduoUOnEyM5npnZchwKPRqtNjMTY+SN78oV01TFCMXMrG4OhR6jWnWt0F1ox30FM0uEQ6HHqFZdKzgUzCw1DoUejdZoVl0rFLfodrPZzFLhUOhR1/SRL2Azs1Q4FHo0W54+MrP1zaHQo9HqjHT6qLhFt0cKZpYKh0KPrNE8urdkesI9BTNLi0Ohx6gbzUX/wtNHZpYKh0IPN5rNbL1zKPSo7ToFTx+ZWSIcCj3qmz7ybS7MLA0OhVyr3WGhEyMNhemJ7O13T8HMUuFQyDVHvMAOgKTsTqkOBTNLhEMhV3xbnx7hSAFgxquvmVlCHAq55vxoF9gpzE56oR0zS4dDIdddn3nUI4UpL8lpZulwKOS6oTA12rdkdnLc1ymYWTIcCrliCmeU1ylAPn3kUDCzRDgUcs2FeqaPZqfGafo6BTNLhEMh16xppDA94UazmaXDoZCrq9GcjRQcCmaWBodCrlHDxWsAs75OwcwS4lDIudFsZuZQ6GrWeZ2CewpmlgiHQq7RajM+JibHNdLjzk6OM7fQodOJkR7XzGwxDoVcM1+fWRp9KADMLfi0VDOrn0Mh12iNdn3mQtHDcF/BzFJQ2aegpB2S9km6bYnnJemPJd0t6RZJZ1VVyzBGvepaYdahYGYJqfKr8U7g3AHPnwecnv9cBlxRYS3LGvWqa4WZKS/JaWbpmKhqxxGxW9LmAZtcAHw4IgK4UdIxko6PiAerqOdLd+3nXdffseTzD3y3wTM2baji0AMVQfTanV/trsRmZraYX3nuyVz6755R6TEqC4UhnAjc1/P3vfljPxAKki4jG01wyimnrOhgG6YnOP24pT/0Tz9uAz93xnEr2vdqbD31WF5+1kk0WgsjP7aZrS0bN0xXfow6Q2FoEbEd2A6wdevWFZ27+ZxTj+U5pz6n1LrKcOyRU7zvl7fUXYaZGVDv2Uf3Ayf3/P2k/DEzM6tJnaFwHXBxfhbS2cBjVfUTzMxsOJVNH0naBWwDNkraC1wOTAJExJXAp4GXAHcDTwKXVFWLmZkNp8qzjy5a5vkAXlfV8c3M7ND5HEgzM+tyKJiZWZdDwczMuhwKZmbWpazfu3ZI2g/cs8Jf3wh8u8Ry1gK/5vXBr3l9WM1rPjUiNi230ZoLhdWQtCcittZdxyj5Na8Pfs3rwyhes6ePzMysy6FgZmZd6y0UttddQA38mtcHv+b1ofLXvK56CmZmNth6GymYmdkADgUzM+taN6Eg6VxJ/yzpbklvqbueqkk6WdLfSbpD0u2S3lh3TaMgaVzSP0q6vu5aRiVfyvZqSf8k6U5JP1V3TVWS9Bv5v+nbJO2SNFN3TVWQtEPSPkm39Tz2VEmfl/Qv+X+PLfu46yIUJI0DfwKcB5wBXCTpjHqrqtwC8JsRcQZwNvC6dfCaAd4I3Fl3ESP2AeCzEfHDwBYO49cv6UTgvwJbI+JHgXHgFfVWVZmdwLl9j70F+EJEnA58If97qdZFKAD/Frg7Ir4ZEfPAx4ELaq6pUhHxYETcnP/5cbIPihPrrapakk4Cfh64qu5aRkXS0cDPAH8GEBHzEfHdequq3AQwK2kCOAJ4oOZ6KhERu4FH+x6+APhQ/ucPAb9Q9nHXSyicCNzX8/e9HOYfkL0kbQbOBL5SbyWV+yPgfwCdugsZodOA/cCf59NmV0k6su6iqhIR9wPvBe4FHiRbsfGGeqsaqeN6Vqh8CDiu7AOsl1BYtyRtAD4JvCkivld3PVWRdD6wLyJuqruWEZsAzgKuiIgzgQNUMKWQinwO/QKyMDwBOFLSr9ZbVT3yhcpKv6ZgvYTC/cDJPX8/KX/ssCZpkiwQPhoR19RdT8WeB7xM0r+STQ++UNJH6i1pJPYCeyOiGAVeTRYSh6sXAd+KiP0R0QKuAX665ppG6WFJxwPk/91X9gHWSyh8FThd0mmSpsgaU9fVXFOlJIlsnvnOiHh/3fVULSLeGhEnRcRmsv+/fxsRh/03yIh4CLhP0rPyh84B7qixpKrdC5wt6Yj83/g5HMaN9UVcB7wm//NrgGvLPkBlazSnJCIWJL0e+BzZ2Qo7IuL2msuq2vOAVwO3Svpa/tjbIuLTNdZk1XgD8NH8C883gUtqrqcyEfEVSVcDN5OdYfePHKa3u5C0C9gGbJS0F7gceDfwCUm/TraEwC+Xflzf5sLMzArrZfrIzMyG4FAwM7Muh4KZmXU5FMzMrMuhYGZmXQ4FMzPrciiYmVmXQ8FslSQ9V9ItkmYkHZnf6/9H667LbCV88ZpZCSS9C5gBZsnuRfT7NZdktiIOBbMS5LeY+CrQBH46Ito1l2S2Ip4+MivH04ANwFFkIwazNckjBbMSSLqO7JbdpwHHR8Tray7JbEXWxV1Szaok6WKgFREfy9cD/7KkF0bE39Zdm9mh8kjBzMy63FMwM7Muh4KZmXU5FMzMrMuhYGZmXQ4FMzPrciiYmVmXQ8HMzLr+Pxaki4oDzG32AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(np.linspace(0,grid_length,grid_points), u);\n", "plt.xlabel('x')\n", "plt.ylabel('u')\n", "plt.title('1D Linear Convection t=0');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Applying the discretization\n", "\n", "For every element of the array `u` we must apply the equation $ u_i^{n+1} = u^n_i - c \\frac{\\Delta t}{\\Delta x}(u_i^n - u^n_{i-1}) $ and store the results in a new temporary array called `new_u` which will become the solution for the next time step. We will repeat this operation for every time step and see how far the wave has moved.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "un = np.ones(grid_points)\n", "\n", "for n in range(nt): #Runs however many timesteps you set earlier\n", " un = u.copy() #copy the u array to not overwrite values\n", " for i in range(1,grid_points):\n", " u[i] = un[i] - c * dt/dx * (un[i]-un[i-1]) " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAGKBJREFUeJzt3XmUHXWd9/H3NzvZSEhiCAkkYQiEzaiERUVBQRZBccEZcWfk4HOOOuqoj7v4uI7jcmRGhcnBGBwVdABHZBRwD8ugBFQgBCSQQBrInpDO1unl+/xxb5dtm+50kq6upPv9Oueee2/V71Z9q3NzP7eqfvW7kZlIkgQwqOoCJEn7DkNBklQwFCRJBUNBklQwFCRJBUNBklQwFFSJiPhoRFxVdR39RUT8LCLeWnUd2v8ZCipExLsiYlFENEXEgk7zTo+ItojYXL81RMQPI+LEbpY3IyIyIoZ0npeZn8/MS0rYjD0SEWdHxMKIaIyINRHx24h4ZdV17UxEfCoivttxWmaem5lXl7CuBRHx2T18bZfvp/r8MyLioYjYGhG/jojpe12w9pqhoI6eAj4LzO9qfmaOBsYApwAPAbdFxBl9VN9ei4jBO5l2IfBfwHeAacBk4JPAK/q2un6ny/dTREwEbgA+ARwELAJ+0KfVaecy05u3v7pR+4+8oNO004GGnbT9OrCoi+XMABIYspN5nwK+26ndW4EngLXAxzq0HQR8GHgUWAf8EDiow/z/AlYCzwALgWM7zFsAXAH8FNgCnNmpjqiv84Pd/D0GAR8HHgdWUwuPA3dVO3AIsK1Trc+ttxlaf/6PwBJgA3ALML1D22OBnwPrgVXAR4FzgB1AM7AZ+FO97W+AS/am3p1s96X19eyor+snvfh+uhS4s8PzUfW/1eyq3/8D/eaegvbWDcDzImJULyzrVOAo4AzgkxFxdH36u4FXAadR+6DdAHyjw+t+BswCngXcC3yv03LfAHyO2h7O7Z3mHQUcClzXTV1vq99eAhwOjKYWht3WnplPAf8LvLZTLddlZnNEXEDtg/41wCTgNuAagIgYA/wCuLm+zUcAv8zMm4HPAz/IzNGZOae36u28kMycR+1v+a/1db2iXttNEbGxi9tNXfwNOzsW+FOHdW2hFvrH9vD1KomhoL31FLVv2+N6YVn/LzO3ZeafqH1gtH/g/R9q32YbMrOJ2l7Ghe3nKjJzfmY2dpg3JyIO7LDcH2fmHZnZlpnbO61zQv3+6W7qeiPw1cx8LDM3Ax8BXt/pXElXtX8fuAggIgJ4fX1a+3Z9ITOXZGYLtQ/759SPrZ8PrMzMr2Tm9vr2/a67P14v1btLmXl+Zo7r4nZ+DxczmtqeXUfPUAtuVchQ0N6aSu1wxMZeWNbKDo+3UvvgAJgO/Kj92yi1wy2twOSIGBwR/xIRj0bEJmB5/TUTOyxrRTfrXFe/n9JNm0OoHYpp9zgwhNq5h13Vfj3w/IiYArwYaKO2R9C+XZd32K711AJ2KrW9l0e7qak7e1NvX9kMjO00bSzQ2Md1qBNDQXvr1cC99d3/sqwAzu30jXREZj5J7XDMBcCZwIHUjplD7cO1XXdDAT9cX/5ru2nzFLUP8HaHAS3UjvN3KzM3ALcC/1Cv9drMbK9nBfCOTtt1QGbeWZ93eFeL3cVq97jenqyr3v11cxe3n/VwuYvpsHdSP/z4d/XpqpChoEJEDImIEcBgYHBEjNhZd9KomRoRlwGXUDsu3p3h9WW133b3fXcl8Ln2LosRMal+PB5qhxuaqH3jH0ntEEyP1T+g/xn4RERcHBFjI2JQRJwaEfPqza4B3hcRMyNiNH85pt/Sw9V8H3gLcCF/OXTUvl0fiYhj69t1YES8rj7vJmBKRLw3IoZHxJiIOLk+bxUwo5u/497W29EqOoVT1rq/ju7idm57u128n34EHBcRr623+SRwX2Y+tAc1qhcZCuro49R6gHwYeFP98cc7zD8kIjZT2/W/GzgeOD0zb93FcjfXl9V+e+lu1nU5cCNwa0Q0AncB7R+Q36F2eORJ4MH6vN2SmddR+yb/j9S+Za+i1mPmx/Um84H/pNazaRmwndrJ7566kdqJ8JX1Y/jt6/0R8EXg2vqhrweAc+vzGoGXUesWuxJ4hNqJY6j1tgJYFxH37mR9e1tvR98Cjqkf4vrv3Xxtl++nzFxDbe/sc9Q6DpxM7XyLKhZ/2ZOVJA107ilIkgqlhUJEHFq/dP3BiFgcEe/ZSZuIiH+LiKURcV9EPK+seiRJu/Y3JxF7UQvw/sy8t34hzj0R8fPMfLBDm3OpHWudRe2Y4hX85VixJKmPlbankJlPZ+a99ceN1PqWT+3U7ALgO1lzFzCu3p9bklSBMvcUChExg9qYL52vyJzKX19Y1FCf9ldXl0bEpdTGSmHUqFEnzJ49u6xSJalfuueee9Zm5qRdtSs9FOr9pK8H3puZm/ZkGfUxWOYBzJ07NxctWtSLFUpS/xcRj++6Vcm9jyJiKLVA+F5m3rCTJk9Su5y/3bT6NElSBcrsfRTULnxZkplf7aLZjcBb6r2QTgGeyczuBiaTJJWozMNHLwTeDNwfEX+sT/sotXFYyMwrqY1x/3JgKbVBuS4usR5J0i6UFgqZeTt/PSjZztok8M6yapAk7R6vaJYkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFQwFSVLBUJAkFUoLhYiYHxGrI+KBLuYfGBE/iYg/RcTiiLi4rFokST1T5p7CAuCcbua/E3gwM+cApwNfiYhhJdYjSdqF0kIhMxcC67trAoyJiABG19u2lFWPJGnXqjyn8HXgaOAp4H7gPZnZtrOGEXFpRCyKiEVr1qzpyxolaUCpMhTOBv4IHAI8B/h6RIzdWcPMnJeZczNz7qRJk/qyRkkaUKoMhYuBG7JmKbAMmF1hPZI04FUZCk8AZwBExGTgKOCxCuuRpAFvSFkLjohrqPUqmhgRDcBlwFCAzLwS+AywICLuBwL4UGauLaseSdKulRYKmXnRLuY/BZxV1volSbvPK5olSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUMBQkSQVDQZJUKC0UImJ+RKyOiAe6aXN6RPwxIhZHxG/LqkWS1DNl7iksAM7pamZEjAO+CbwyM48FXldiLZKkHigtFDJzIbC+myZvAG7IzCfq7VeXVYskqWeqPKdwJDA+In4TEfdExFu6ahgRl0bEoohYtGbNmj4sUZIGlipDYQhwAnAecDbwiYg4cmcNM3NeZs7NzLmTJk3qyxolaUAZUuG6G4B1mbkF2BIRC4E5wJ8rrEmSBrQq9xR+DJwaEUMiYiRwMrCkwnokacArbU8hIq4BTgcmRkQDcBkwFCAzr8zMJRFxM3Af0AZclZlddl+VJJWvtFDIzIt60OZLwJfKqkGStHu8olmSVDAUJEkFQ0GSVDAUJEkFQ0GSVDAUJEkFQ0GSVDAUJEkFQ0GSVDAUJEkFQ0GSVDAUJEkFQ0GSVDAUJEkFQ0GSVOjR7ylExCd3Nj0zP9275UiSqtTTH9nZ0uHxCOB8/OlMSep3ehQKmfmVjs8j4svALaVUJEmqzJ6eUxgJTOvNQiRJ1evpOYX7gaw/HQxMAjyfIEn9TE/PKZzf4XELsCozW0qoR5JUoZ6eU3i87EIkSdXzOgVJUsFQkCQVDAVJUsFQkCQVDAVJUsFQkCQVDAVJUsFQkCQVDAVJUsFQkCQVDAVJUsFQkCQVDAVJUqG0UIiI+RGxOiIe2EW7EyOiJSIuLKsWSVLPlLmnsAA4p7sGETEY+CJwa4l1SJJ6qLRQyMyFwPpdNHs3cD2wuqw6JEk9V9k5hYiYCrwauKIHbS+NiEURsWjNmjXlFydJA1SVJ5q/BnwoM9t21TAz52Xm3MycO2nSpD4oTZIGpp7+RnMZ5gLXRgTARODlEdGSmf9dYU2SNKBVFgqZObP9cUQsAG4yECSpWqWFQkRcA5wOTIyIBuAyYChAZl5Z1nolSXuutFDIzIt2o+3byqpDktRzXtEsSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkQmmhEBHzI2J1RDzQxfw3RsR9EXF/RNwZEXPKqkWS1DNl7iksAM7pZv4y4LTMPB74DDCvxFokST0wpKwFZ+bCiJjRzfw7Ozy9C5hWVi2SpJ7ZV84pvB34WVczI+LSiFgUEYvWrFnTh2VJ0sBSeShExEuohcKHumqTmfMyc25mzp00aVLfFSdJA0xph496IiKeDVwFnJuZ66qsRZJU4Z5CRBwG3AC8OTP/XFUdkrQ/eGZbM89sbS59PaXtKUTENcDpwMSIaAAuA4YCZOaVwCeBCcA3IwKgJTPnllWPJO2Plq7ezNV3Luf6exu4+IUz+ODZs0tdX5m9jy7axfxLgEvKWr8k7a/a2pLf/Hk1375jObc9spZhQwZxwZxDOO/4Q0pfd6XnFCRJf9G4vZnr7mng6juXs3zdViaPHc4Hzz6K1594KBNGD++TGgwFSarYhi07uPK3j/Ldux5ny45WTpg+nvefdRTnHHcwQwf37alfQ0GSKtK4vZmrblvGt25fxpYdLbzi2Yfw9lNnMufQcZXVZChIUh/btqOVq/93OVf+9lE2bm3m3OMO5p9fdiSzJo+pujRDQZL6SlNLK9f+fgVf//VS1jQ2cdqRk/jAWUdx/LQDqy6tYChIUskyk5/c9zRf/NlDPLlxGyfNPIhvvvF5nDjjoKpL+xuGgiSV6L6GjXz6Jw+y6PENHDNlLF94zfG8aNZE6tdn7XMMBUkqwepN2/nXWx7munsamDh6GF987fFceMKhDB60b4ZBO0NBknrR9uZWvnX7Mr7x66W0tCbvOO1w3vWSIxgzYmjVpfWIoSBJvSAzufmBlXzup0to2LCNs46ZzMfOO5rpE0ZVXdpuMRQkaS89sqqRT/1kMXcsXcfsg8fw/UtO5gVHTKy6rD1iKEjSHtq0vZnLf/EIV9+5nJHDBvPpC47lDScdxpA+vgq5NxkKkrSb2tqS6+9t4Is3P8y6LU28/sRD+cBZR/XZ+ERlMhQkaTfc17CRy25czB+e2MhzDxvH/LfN5dnTqhuWorcZCpLUA2sam/jyLQ/zw3tWMGHUcL78ujm85rlTGbSPdzHdXYaCJHWjqaWV+bcv5xu/XkpTSyuXnDqTd58xi7H7SRfT3WUoSNJOZCa3LF7J53/6EE+s38qZRz+Lj513DDMn7l9dTHeXoSBJnTz41CY+fdNi7npsPUdOHs1/vv0kXjRrUtVl9QlDQZLqVm3aztd+8WeuvXsF4w4YymcuOJaL9vMuprvLUJA04G3cuoMrfvsoC+5YTmtb8rYXzOC9ZxzJgSP753mD7hgKkgasLU0tfPuOZfzHwsfY3NTCq54zlfedeSSHTRhZdWmVMRQkDTg7Wtq45vdP8O+/eoS1m3dw5tGT+cDZRzL74LFVl1Y5Q0HSgNHU0sqP//AU//arR2jYsI2TZx7Ef7x5NidMH191afsMQ0FSv9e4vZnv/+4J5t+xjFWbmjhu6lg+9+rjefE+/GM3VTEUJPVbqxu38+07lvPdux6ncXsLL/i7CXzpwjn79C+fVc1QkNTvLFu7hXkLH+P6extobm3j5cdN4R2nHd6vxigqi6EgqV9oamnlFw+u5tq7n+D2pWsZOngQF54wjUtfdDgz+vlVyL3JUJC0X3tkVSM/uHsFN/zhSdZv2cEhB47gn146izeechjPGjOi6vL2O4aCpP3OlqYW/uf+p/nB3Su45/ENDB0cvOyYyfzDiYdx6hETGdzPRi7tS4aCpP3C+i07+OWSVdyyeBW3PbKGppY2jnjWaD728qN59fOmMrEf/MDNvsBQkLTPWrF+K7c+uIpbF6/k7uXraUuYOu4ALjrpMM5/9hROmD7eXkS9zFCQtM/YuHUHi5Zv4PfL13PbI2tZ8vQmAGYfPIZ3veQIzjr2YI49ZKxBUCJDQVJlVm3azu+Xrefu5ev5/bL1PLSyEYBhgwfxnMPG8fHzjuZlx0xm+gR7D/UVQ0FS6VrbksfXbeGhlY2129ObWLJyEyvWbwNg5LDBnDB9POc/ewonzjiIOYeOY8TQwRVXPTAZCpJ6RWaydvMOGjZspWHDNho2bOOxNZt5eFUjf17VyPbmNgAGBcyYOIrjpx7IW58/g5NmHsQxU8YOqN8s2JeVFgoRMR84H1idmcftZH4AlwMvB7YCb8vMe8uqR9KeaW1LNm7dwbotO1i3eQfrtjSxfssO1m7ewdrNTTy5YRsNG7by5MZtxQd/u4mjhzP74DG86eTpHHXwGGYfPJZZk0e7F7APK3NPYQHwdeA7Xcw/F5hVv50MXFG/lwaszKQtO9yTZEImtGbS2pa0tSWt2eE+obU1aW5ro6U1aW5to6UtaWlto7k1aWlro6m5je0trX9zv725je3NrWxuamHz9pbivrGphc1NzWze3sIz25ppy7+tNQLGjxzG1HEHcOTkMbx09rOYNn4k08YfwLTxI5k6/gBGD/dgxP6mtH+xzFwYETO6aXIB8J3MTOCuiBgXEVMy8+ky6rn5gZW8/4d/LGPR2o/s5LOt56/dxYuzw9Lb2+ZfNyg+5NvnZWb9fi8K20vDhwxizIghjB4+hNH1+6njDmDMiDGMGj6Y8SOHMWHUMA4aPZyJo4Zx0OhhTBg1nPEjh3rIpx+qMsanAis6PG+oT/ubUIiIS4FL6083R8TDe7jOicDaPXzt/sptHhjc5oFhb7Z5ek8a7Rf7dpk5D5i3t8uJiEWZObcXStpvuM0Dg9s8MPTFNle57/ckcGiH59Pq0yRJFakyFG4E3hI1pwDPlHU+QZLUM2V2Sb0GOB2YGBENwGXAUIDMvBL4KbXuqEupdUm9uKxaOtjrQ1D7Ibd5YHCbB4bStzmyym4PkqR9iv3JJEkFQ0GSVBgwoRAR50TEwxGxNCI+XHU9ZYuIQyPi1xHxYEQsjoj3VF1TX4iIwRHxh4i4qepa+kr9ws/rIuKhiFgSEc+vuqYyRcT76u/pByLimojol7+5GRHzI2J1RDzQYdpBEfHziHikfj++t9c7IEIhIgYD36A2tMYxwEURcUy1VZWuBXh/Zh4DnAK8cwBsM8B7gCVVF9HHLgduzszZwBz68fZHxFTgn4C59THVBgOvr7aq0iwAzuk07cPALzNzFvDL+vNeNSBCATgJWJqZj2XmDuBaasNs9FuZ+XT7AIOZ2Ujtg2JqtVWVKyKmAecBV1VdS1+JiAOBFwPfAsjMHZm5sdqqSjcEOCAihgAjgacqrqcUmbkQWN9p8gXA1fXHVwOv6u31DpRQ6GpIjQGhPgbVc4HfVVtJ6b4G/F+gbVcN+5GZwBrg2/XDZldFRL/9RZrMfBL4MvAEtSFxnsnMW6utqk9N7nA910pgcm+vYKCEwoAVEaOB64H3ZuamquspS0S0D9N+T9W19LEhwPOAKzLzucAWSjiksK+oH0O/gFoYHgKMiog3VVtVNeqDifb6NQUDJRQG5JAaETGUWiB8LzNvqLqekr0QeGVELKd2ePClEfHdakvqEw1AQ2a27wVeRy0k+qszgWWZuSYzm4EbgBdUXFNfWhURUwDq96t7ewUDJRTuBmZFxMyIGEbtxNSNFddUqvqPGH0LWJKZX626nrJl5kcyc1pmzqD27/urzOz33yAzcyWwIiKOqk86A3iwwpLK9gRwSkSMrL/Hz6Afn1jfiRuBt9YfvxX4cW+vYL8YJXVvZWZLRLwLuIVab4X5mbm44rLK9kLgzcD9EdH+QxIfzcyfVliTyvFu4Hv1LzyP0TdDxlQiM38XEdcB91LrYfcH+ulwF10MFfQvwA8j4u3A48Df9/p6HeZCktRuoBw+kiT1gKEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEg7aWIODEi7ouIERExqj7W/3FV1yXtCS9ek3pBRHwWGAEcQG0soi9UXJK0RwwFqRfUh5i4G9gOvCAzWysuSdojHj6SescEYDQwhtoeg7Rfck9B6gURcSO1IbtnAlMy810VlyTtkQExSqpUpoh4C9Ccmd+v/x74nRHx0sz8VdW1SbvLPQVJUsFzCpKkgqEgSSoYCpKkgqEgSSoYCpKkgqEgSSoYCpKkwv8HZwsQBwmnFNgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(np.linspace(0,grid_length,grid_points), u);\n", "plt.ylim(1,2);\n", "plt.xlabel('x')\n", "plt.ylabel('u')\n", "plt.title('1D Linear Convection t=10');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Results\n", "\n", "As we can see above the wave has indeed been transported from its initial position at around 2 towards 10 and is currently halfway out of the plot. However it has lost its square wave shape. Because I was having a hard time understanding how that happened I thought an animation of each timestep would be a good way to see how the wave has been transported over time and how its shape changed.\n", "\n", "\n", "## Animating the wave moving\n", "\n", "To animate the moving of the wave and be able to interact a bit with it in this jupyter notebook I followed a [fantastic guide](http://tiao.io/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/) by Louis Tiao that outlines the process in detail. The following code was taken from that guide and aims to set up the figure object that contains the animation while also re-initializing the U array back to its initial square-wave configuration." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAFwxJREFUeJzt3Xu0ZGV95vHvQzfaAi1E6SFcFDAjIiQhajveEtNRR8GgjDNmIk5EvKETNZq4xnuEFdExS82oY5R0kGkJiFFEIY5GWCrBxPHSIHI1SkCluXUDolxUbPnNH3s3XV285/Q5fU71Pt39/axV61Tt/dZ+f3tX1X5qX2qfVBWSJI3baegCJEkLkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA2KWkrw5yclD17G9SPL5JC/cSn1dnmTF1uhrEpKcn+Slc3j+HUkeNp819dN9RJKLk9ye5E/me/oaznYREElelWR1kp8nWTU2bkWSe/oPxx1J1iT5RJLHTjO9A5JUksXj46rqnVW1xR/S+ZbkGUku6D+c65L8U5JnD11XS5ITkpw2Oqyqjqiqj06gr1VJThzr69CqOn+++1qIWmFSVbtV1dUT6O71wJeramlVfWAC05/WXIKz/6x/OcldSb6T5GnTtL3Pe2or1vn2JJcmWZ/khMb45yf5QZI7k3wmyYO2tM5R20VAANcDJwKnTDW+qnYDlgKPB74DfCXJU7dSfXOWZFFj2HOBTwKnAvsBewFvA561davTDm5/4PKpRrbeuwvIGcC3gAcDbwHOTLJs2JKarqIL4v87PiLJocDfAC+gWwfcBXxoXnqtqu3mRhcSq8aGrQDWNNp+EFg9xXQOAApY3Bh3AnDaWLsXAj8EbgbeMtJ2J+CNwL8BtwCfAB40Mv6TwI3Aj4ELgENHxq0CPgx8DrgTeNpYHen7/B/TLI+dgLcCPwDW0gXJ7purHdgH+OlYrY/q2+zcP34xcCXwI+ALwP4jbQ8FzgNuBW4C3gwcDtwN/AK4A/h23/Z84KVzqbcx38f1/dzd9/UP/fDvb1iO/ev4SeA04HbgUuAg4E1939cCTx+Z5u7AR4AbgOvo3muLplnuzdcd+DzwqrH23wb+c3//icA3+/fEN4EnjrQbXVYn0L8Px9+zwDuAXwI/6+f/g32bAv79yPycCqzrl/dbgZ36cccC/wy8p399rwGOmGJevzTW10E03rsz6O9fgP8F3AZc3S+HY/vXYS3wwin6b87rDNcXBwE/B5aODPsK8IpZvKce2b8ut9GF5LPnu86x6ZwGnDA27J3Ax0Ye/1pf59It6WOTac91AgvpxuwC4inAPcCujXH3ftga4+79YI60+1vgAcBh/Rvukf341wBfo/t2f3+6lD9jZFovptuquT/wPuDikXGr6FYST6Jb4SwZq+Pgvu8Dp1keL6b75vEwYDfgLODvZlj7l4CXjUzr3cBJ/f2j+uk+km6F9Fbgq/24pXQr0dcBS/rHjxtfdiPTPZ+NK70trrcx76uAE8eGfZ9NA+JnwDP6eTiVbkX4FmBn4GXANSPP/XT/+u0K/DvgG8DLp+h7ytcdOAb4l5G2h9CtXO4PPIhuhfyCvqaj+8cPbiyrTZYlY+/Z0bYjbUYD4lTg7P71OQD4LvCSftyxdCvDlwGLgP9Ot5WeKeZ3k75ovHdn0N964EV9fyfSfQn46365PJ0uxHebSf/9sEv65dq6fahv8xzgyrHnfRD43zN5T/Xvk6vovgDdj26dcjvwiPmsc6x9KyDOBt4wNuwO4DGzWX82a57rBBbSjdkFxIYV7L6NcZt82MbGncB9A2K/kfHfAJ7X378SeOrIuL37D15runv009rwjXkVcOo08/qkvv2Sadp8EfjjkceP2ND/DGp/KfCl/n7ovsk9uX/8efoPd/94J7rN2v3pVmrfmqKee5fdyLB7PzRzqbfR1yo2HxDnjYx7Vv+hWtQ/Xtr3twfdZvvPgQeMtD+abr97q+8pX/d+unfSb3HRfbM8pb//AuAbY9P6f8CxjWW1ybJkFgFBtxK+GzhkZNzLgfP7+8cCV42M26V/7q9OMb+b9MXYe3eG/X1vZNxv9P3tNTLsFuC3ZtL/TG/98v7a2LB3MLYOmeo9BfwO3R6AnUaGncHYCnyudY5NoxUQX2Rsq4duK3fFXPqqqu3mGMSW2JfuTXjbPEzrxpH7d9F9+4VuhfnpJLcluY1uxfFLYK8ki5K8K8m/JfkJ3coLYM+RaV07TZ+39H/3nqbNPnSb8xv8gG4ltdcMav8U8IQkewNPptva+srIfL1/ZL5upQuRfYGH0O1a2RJzqXdL3DRy/6fAzVX1y5HH9NPfn+7b4g0j8/w3dFsSLVO+7lV1O91+5Of1bY8GTu/vj88//eN9t2TmprEn3fyML+vRfu5dzlV1V393Nst69L07k/7GXwuqanzYXF7rljuAB44NeyDdVsBM7ANcW1X3jAybxOu1OXOdjyntyAHxHOCiqrpzgn1cS7fvdo+R25Kqug54Pt2umg37Zw/on5OR59c00/7Xfvr/ZZo219OtrDZ4KN2m/E3t5iMdV/0IOBf4w77Wj1f/1aTv9+Vj8/WAqvpqP26qUymnm5851bsFfc3GtXRbEHuOzO8Dq+rQadpP9bpD9y3z6CRPoNv98uV++Pj8Q7cMruO+7qT7Zr/Br46Nn27+b6bbohlf1q1+ttRo/5Pu7z7z2p/SfMcUt5P6ZpcDD0uydOSphzH1Affxfq4HHpJkdD063XxtaZ2bc3lf94ZpPoxu19x3Z/j8KW0XAZFkcZIldJuyi5IsaZ2ims6+SY6n24Xy5s1M+v79tDbcZru8TgLekWT/vv9lSY7qxy2lW+ncQvdBf+dsJtyvrP8M+PMkL0rywCQ7JfntJCv7ZmcAf5rkwCS79X38fVWtn2E3H6PbZ/7c/v7ofL2pP3uCJLsn+YN+3GeBvZO8Nsn9kyxN8rh+3E3AAdMsx7nWO+ompg6qWamqG+jC8r0jy/nXkvzuFE+Z7nWH7uDt/sBf0M3fPSPDD+pPWVyc5A/pjlF8ttHHxcCTkzw0ye50B9dHTTn//VbSJ/oal/Z1/hnd7ot5txX6u8+8VndK825T3F7Rt/ku3XI8vv98Pwf4Tbqt55n083W6rdjXJ9k53W9sngV8fD7rBOinv4Runb24r3fD2WGnA89K8jtJdqV7X53Vb63OyXYREHQHSX9Kd+bIH/X33zoyfp8kd9Btin2Tbh/niqo6dzPTvaOf1obbU2ZZ1/uBc4Bzk9xOd+Byw8ryVLrN0euAK/pxs1JVZ9J9w38x3beZm+iOw5zdNzkF+Du6M6SuoTso++pZdHEO8HDgxqr69ki/nwb+Evh4v3vsMuCIftztwH+k+6DcCHwP+L3+qZ/s/96S5KJGf3Otd9RHgEP63Tyf2cJpjDqG7kDkFXQHjs9k6t17073uVNXP6Q7AP42R4K2qW4Aj6Q7w30J3WuORVXXzeAdVdR7w93QHOS/kviHyfuC5SX6UpPXbhFfTbYVcTXfG0seY+jTx+TDJ/jY3r9N5HrCc7jV9F/Dcqlo3RdtN3lNVdTfd+/wIuq2kDwHHVNV3JlDn39Ktg46mO5Hip3THUKiqy4FX0AXFWrovn388y+k3ZeNeA0mSNtpetiAkSfNsYgGR5CHpfsJ+RX8g5jWNNknygSRXJbkkyaMnVY8kaXbucyB3Hq0HXldVF/VnCVyY5LyqumKkzRF0+7gfTreP9sOM7KuVJA1nYlsQVXVDVV3U37+d7lzw8fODj6L7QU1V1deAPfrz7iVJA5vkFsS9khxAdy2fr4+N2pdNf1Czph92w9jzj6O7Fgq77rrrYw4++OBJlSpJ26ULL7zw5qqa1YUIJx4Q/fnsnwJeW1U/2ZJpVNVKYCXA8uXLa/Xq1fNYoSRt/5KM/0p/syZ6FlOSnenC4fSqOqvR5Dq6SzNssB/z+2tOSdIWmuRZTKH7YcmVVfVXUzQ7BzimP5vp8cCP+1+tSpIGNsldTE+i+6XfpUku7oe9me5aJVTVSXSXFngm3SVz76K73K8kaQGYWEBU1T+z6YXnWm0KeOWkapAkbTl/SS1JajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKlpYgGR5JQka5NcNsX43ZP8Q5JvJ7k8yYsmVYskafYmuQWxCjh8mvGvBK6oqsOAFcB7k9xvgvVIkmZhYgFRVRcAt07XBFiaJMBufdv1k6pHkjQ7Qx6D+CDwSOB64FLgNVV1T6thkuOSrE6yet26dVuzRknaYQ0ZEM8ALgb2AX4L+GCSB7YaVtXKqlpeVcuXLVu2NWuUpB3WkAHxIuCs6lwFXAMcPGA9kqQRQwbED4GnAiTZC3gEcPWA9UiSRiye1ISTnEF3dtKeSdYAxwM7A1TVScDbgVVJLgUCvKGqbp5UPZKk2ZlYQFTV0ZsZfz3w9En1L0maG39JLUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNEwuIJKckWZvksmnarEhycZLLk/zTpGqRJM3eJLcgVgGHTzUyyR7Ah4BnV9WhwB9MsBZJ0ixNLCCq6gLg1mmaPB84q6p+2LdfO6laJEmzN+QxiIOAX0lyfpILkxwzVcMkxyVZnWT1unXrtmKJkrTjGjIgFgOPAX4feAbw50kOajWsqpVVtbyqli9btmxr1ihJO6zFA/a9Brilqu4E7kxyAXAY8N0Ba5Ik9Ybcgjgb+O0ki5PsAjwOuHLAeiRJIya2BZHkDGAFsGeSNcDxwM4AVXVSVV2Z5B+BS4B7gJOraspTYiVJW9fEAqKqjp5Bm3cD755UDZKkLecvqSVJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVLTjP4fRJK3tYZX1V/MbzmSpIVipv8w6M6R+0uAI/Hfg0rSdm1GAVFV7x19nOQ9wBcmUpEkaUHY0mMQuwD7zWchkqSFZabHIC4Fqn+4CFgGePxBkrZjMz0GceTI/fXATVW1fgL1SJIWiJkeg/jBpAuRJC0s/g5CktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTRMLiCSnJFmb5LLNtHtskvVJnjupWiRJszfJLYhVwOHTNUiyCPhL4NwJ1iFJ2gITC4iqugC4dTPNXg18Clg7qTokSVtmsGMQSfYFngN8eAZtj0uyOsnqdevWTb44SdKgB6nfB7yhqu7ZXMOqWllVy6tq+bJly7ZCaZKkmf5P6klYDnw8CcCewDOTrK+qzwxYkySpN1hAVNWBG+4nWQV81nCQpIVjYgGR5AxgBbBnkjXA8cDOAFV10qT6lSTNj4kFRFUdPYu2x06qDknSlvGX1JKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUNLGASHJKkrVJLpti/H9LckmSS5N8Nclhk6pFkjR7k9yCWAUcPs34a4DfrarfAN4OrJxgLZKkWVo8qQlX1QVJDphm/FdHHn4N2G9StUiSZm+hHIN4CfD5qUYmOS7J6iSr161btxXLkqQd1+ABkeT36ALiDVO1qaqVVbW8qpYvW7Zs6xUnSTuwie1imokkvwmcDBxRVbcMWYskaVODbUEkeShwFvCCqvruUHVIktomtgWR5AxgBbBnkjXA8cDOAFV1EvA24MHAh5IArK+q5ZOqR5I0O5M8i+nozYx/KfDSSfUvSZqbwQ9SS5IWJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpomFhBJTkmyNsllU4xPkg8kuSrJJUkePalaJEmzN8ktiFXA4dOMPwJ4eH87DvjwBGuRJM3SxAKiqi4Abp2myVHAqdX5GrBHkr0nVY8kaXYWD9j3vsC1I4/X9MNuGG+Y5Di6rQyAn0+122oHtCdw89BFLBAui41cFhu5LDZ6xGyfMGRAzFhVrQRWAiRZXVXLBy5pQXBZbOSy2MhlsZHLYqMkq2f7nCHPYroOeMjI4/36YZKkBWDIgDgHOKY/m+nxwI+r6j67lyRJw5jYLqYkZwArgD2TrAGOB3YGqKqTgM8BzwSuAu4CXjTDSa+c92K3XS6LjVwWG7ksNnJZbDTrZZGqmkQhkqRtnL+kliQ1GRCSpKZtKiCSHJ7kX/vLc7xx6HqGkuQhSb6c5Ioklyd5zdA1DSnJoiTfSvLZoWsZWpI9kpyZ5DtJrkzyhKFrGkqSP+0/H5clOSPJkqFr2lpalzpK8qAk5yX5Xv/3VzY3nW0mIJIsAv6a7hIdhwBHJzlk2KoGsx54XVUdAjweeOUOvCwAXgNcOXQRC8T7gX+sqoOBw9hBl0uSfYE/AZZX1a8Di4DnDVvVVrWK+17q6I3AF6vq4cAX+8fT2mYCAvgPwFVVdXVV3Q18nO5yHTucqrqhqi7q799OtxLYd9iqhpFkP+D3gZOHrmVoSXYHngx8BKCq7q6q24atalCLgQckWQzsAlw/cD1bzRSXOjoK+Gh//6PAf9rcdLalgJjq0hw7tCQHAI8Cvj5sJYN5H/B64J6hC1kADgTWAf+n3+V2cpJdhy5qCFV1HfAe4Id0l+/5cVWdO2xVg9tr5LdmNwJ7be4J21JAaEyS3YBPAa+tqp8MXc/WluRIYG1VXTh0LQvEYuDRwIer6lHAncxgN8L2qN+/fhRdaO4D7Jrkj4atauGo7vcNm/2Nw7YUEF6aY0SSnenC4fSqOmvoegbyJODZSb5Pt8vxKUlOG7akQa0B1lTVhq3JM+kCY0f0NOCaqlpXVb8AzgKeOHBNQ7tpwxWz+79rN/eEbSkgvgk8PMmBSe5Hd8DpnIFrGkSS0O1nvrKq/mroeoZSVW+qqv2q6gC698OXqmqH/ZZYVTcC1ybZcNXOpwJXDFjSkH4IPD7JLv3n5ansoAfsR5wDvLC//0Lg7M09YZu4mitAVa1P8irgC3RnJJxSVZcPXNZQngS8ALg0ycX9sDdX1ecGrEkLw6uB0/svUVcz80vYbFeq6utJzgQuojvr71vsQJfdmOJSR+8CPpHkJcAPgP+62el4qQ1JUsu2tItJkrQVGRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaENEdJHpvkkiRLkuza/w+CXx+6Lmmu/KGcNA+SnAgsAR5Adz2k/zlwSdKcGRDSPOgvbfFN4GfAE6vqlwOXJM2Zu5ik+fFgYDdgKd2WhLTNcwtCmgdJzqG75PiBwN5V9aqBS5LmbJu5mqu0UCU5BvhFVX2s/9/pX03ylKr60tC1SXPhFoQkqcljEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqen/A6szcbwReQ0BAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Imports for animation and display within a jupyter notebook\n", "from matplotlib import animation, rc \n", "from IPython.display import HTML\n", "\n", "#Generating the figure that will contain the animation\n", "fig, ax = plt.subplots()\n", "fig.set_size_inches(11, 7)\n", "ax.set_xlim(( 0, grid_length))\n", "ax.set_ylim((1, 2))\n", "line, = ax.plot([], [], lw=2)\n", "plt.xlabel('x')\n", "plt.ylabel('u')\n", "plt.title('1D Linear Convection time evolution from t=0 to t=10');\n", "#Resetting the U wave back to initial conditions\n", "u = np.ones(grid_points)\n", "u[int(.5/ dx):int(1 / dx + 1)] = 2" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "#Initialization function for funcanimation\n", "def init():\n", " line.set_data([], [])\n", " return (line,)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "#Main animation function, each frame represents a time step in our calculation\n", "def animate(j):\n", " x = np.linspace(0, grid_length, grid_points)\n", " un = u.copy() #copy the u array to not overwrite values\n", " for i in range(1,grid_points):\n", " u[i] = un[i] - c * dt/dx * (un[i]-un[i-1]) \n", " line.set_data(x, u)\n", " return (line,)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "anim = animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=nt, interval=20)\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion\n", "\n", "The above animation can be played forwards or backwards in time and illustrates the behaviour of our numerical linear convection simulation over each time step. This behaviour is not exactly what we expected since the initial conditions are not being moved along the wave since the shape is not perfectly preserved. Truncation error is the most likely culprit as to why there are loses over time of the original condition and will be a topic discussed on the first bonus lecture.\n", "\n", " Next we will look at nonlinear convection using the same set of tools." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }