\u001b[0m in \u001b[0;36mtest_convrate_sincos\u001b[0;34m()\u001b[0m\n\u001b[1;32m 17\u001b[0m \u001b[0;31m# print('rates sin(x)*cos(t) solution:', \\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 18\u001b[0m \u001b[0;31m# [round(r_,2) for r_ in r])\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m0.002\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 21\u001b[0m \u001b[0mtest_convrate_sincos\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mAssertionError\u001b[0m: "
]
}
],
"source": [
"def test_convrate_sincos():\n",
" n = m = 2\n",
" L = 1.0\n",
" u_exact = lambda x, t: np.cos(m*np.pi/L*t)*np.sin(m*np.pi/L*x)\n",
"\n",
" r = convergence_rates(\n",
" u_exact=u_exact,\n",
" I=lambda x: u_exact(x, 0),\n",
" V=lambda x: 0,\n",
" f=0,\n",
" c=1,\n",
" L=L,\n",
" dt0=0.1,\n",
" num_meshes=6,\n",
" C=0.9,\n",
" T=1)\n",
"# print('rates sin(x)*cos(t) solution:', \\\n",
"# [round(r_,2) for r_ in r])\n",
" assert abs(r[-1] - 2) < 0.002\n",
" \n",
"test_convrate_sincos()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Doing `py.test -s -v wave1D_u0.py` will run also this test function and\n",
"show the rates 2.05, 1.98, 2.00, 2.00, and 2.00 (to two decimals).\n",
"\n",
"\n",
"## Visualization: animating the solution\n",
"\n",
"\n",
"Now that we have verified the implementation it is time to do a\n",
"real computation where we also display evolution of the waves\n",
"on the screen. Since the `solver` function knows nothing about\n",
"what type of visualizations we may want, it calls the callback function\n",
"`user_action(u, x, t, n)`. We must therefore write this function and\n",
"find the proper statements for plotting the solution.\n",
"\n",
"### Function for administering the simulation\n",
"\n",
"The following `viz` function\n",
"\n",
"1. defines a `user_action` callback function\n",
" for plotting the solution at each time level,\n",
"\n",
"2. calls the `solver` function, and\n",
"\n",
"3. combines all the plots (in files) to video in different formats."
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"def viz(\n",
" I, V, f, c, L, dt, C, T, # PDE parameters\n",
" umin, umax, # Interval for u in plots\n",
" animate=True, # Simulation with animation?\n",
" tool='matplotlib', # 'matplotlib' or 'scitools'\n",
" solver_function=devito_solver, # Function with numerical algorithm\n",
" ):\n",
" \"\"\"Run solver and visualize u at each time level.\"\"\"\n",
"\n",
" def plot_u_st(u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" plt.plot(x, u, 'r-',\n",
" xlabel='x', ylabel='u',\n",
" axis=[0, L, umin, umax],\n",
" title='t=%f' % t[n], show=True)\n",
" # Let the initial condition stay on the screen for 2\n",
" # seconds, else insert a pause of 0.2 s between each plot\n",
" time.sleep(2) if t[n] == 0 else time.sleep(0.2)\n",
" plt.savefig('frame_%04d.png' % n) # for movie making\n",
"\n",
" class PlotMatplotlib:\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" if n == 0:\n",
" plt.ion()\n",
" self.lines = plt.plot(x, u, 'r-')\n",
" plt.xlabel('x'); plt.ylabel('u')\n",
" plt.axis([0, L, umin, umax])\n",
" plt.legend(['t=%f' % t[n]], loc='lower left')\n",
" else:\n",
" self.lines[0].set_ydata(u)\n",
" plt.legend(['t=%f' % t[n]], loc='lower left')\n",
" plt.draw()\n",
" time.sleep(2) if t[n] == 0 else time.sleep(0.2)\n",
" plt.savefig('tmp_%04d.png' % n) # for movie making\n",
"\n",
" if tool == 'matplotlib':\n",
" import matplotlib.pyplot as plt\n",
" plot_u = PlotMatplotlib()\n",
" elif tool == 'scitools':\n",
" import scitools.std as plt # scitools.easyviz interface\n",
" plot_u = plot_u_st\n",
" import time, glob, os\n",
"\n",
" # Clean up old movie frames\n",
" for filename in glob.glob('tmp_*.png'):\n",
" os.remove(filename)\n",
"\n",
" # Call solver and do the simulaton\n",
" user_action = plot_u if animate else None\n",
" u, x, t, cpu = solver_function(\n",
" I, V, f, c, L, dt, C, T, user_action)\n",
"\n",
" # Make video files\n",
" fps = 4 # frames per second\n",
" codec2ext = dict(flv='flv', libx264='mp4', libvpx='webm',\n",
" libtheora='ogg') # video formats\n",
" filespec = 'tmp_%04d.png'\n",
" movie_program = 'ffmpeg'\n",
" for codec in codec2ext:\n",
" ext = codec2ext[codec]\n",
" cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\\\n",
" '-vcodec %(codec)s movie.%(ext)s' % vars()\n",
" os.system(cmd)\n",
"\n",
" if tool == 'scitools':\n",
" # Make an HTML play for showing the animation in a browser\n",
" plt.movie('tmp_*.png', encoder='html', fps=fps,\n",
" output_file='movie.html')\n",
" return cpu"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Dissection of the code\n",
"\n",
"The `viz` function can either use SciTools or Matplotlib for\n",
"visualizing the solution. The `user_action` function based on SciTools\n",
"is called `plot_u_st`, while the `user_action` function based on\n",
"Matplotlib is a bit more complicated as it is realized as a class and\n",
"needs statements that differ from those for making static plots.\n",
"SciTools can utilize both Matplotlib and Gnuplot (and many other\n",
"plotting programs) for doing the graphics, but Gnuplot is a relevant\n",
"choice for large $N_x$ or in two-dimensional problems\n",
"as Gnuplot is significantly faster than\n",
"Matplotlib for screen animations.\n",
"\n",
"\n",
"A function inside another function, like `plot_u_st` in the above code\n",
"segment, has access to *and remembers* all the local variables in the\n",
"surrounding code inside the `viz` function (!). This is known in\n",
"computer science as a *closure* and is very convenient to program\n",
"with. For example, the `plt` and `time` modules defined outside\n",
"`plot_u` are accessible for `plot_u_st` when the function is called\n",
"(as `user_action`) in the `solver` function. Some may think, however,\n",
"that a class instead of a closure is a cleaner and\n",
"easier-to-understand implementation of the user action function, see\n",
"the section [wave:pde2:software](#wave:pde2:software).\n",
"\n",
"The `plot_u_st` function just makes a standard SciTools `plot` command\n",
"for plotting `u` as a function of `x` at time `t[n]`. To achieve a\n",
"smooth animation, the `plot` command should take keyword arguments\n",
"instead of being broken into separate calls to `xlabel`, `ylabel`,\n",
"`axis`, `time`, and `show`. Several `plot` calls will automatically\n",
"cause an animation on the screen. In addition, we want to save each\n",
"frame in the animation to file. We then need a filename where the\n",
"frame number is padded with zeros, here `tmp_0000.png`,\n",
"`tmp_0001.png`, and so on. The proper printf construction is then\n",
"`tmp_%04d.png`.\n",
"% if BOOK == \"book\":\n",
"the section [vib:ode1:anim](#vib:ode1:anim) contains more basic\n",
"information on making animations.\n",
"% endif\n",
"\n",
"The solver is called with an argument `plot_u` as `user_function`.\n",
"If the user chooses to use SciTools, `plot_u` is the `plot_u_st`\n",
"callback function, but for Matplotlib it is an instance of the\n",
"class `PlotMatplotlib`. Also this class makes use of variables\n",
"defined in the `viz` function: `plt` and `time`.\n",
"With Matplotlib, one has to make the first plot the standard way, and\n",
"then update the $y$ data in the plot at every time level. The update\n",
"requires active use of the returned value from `plt.plot` in the first\n",
"plot. This value would need to be stored in a local variable if we\n",
"were to use a closure for the `user_action` function when doing the\n",
"animation with Matplotlib. mathcal{I}_t is much easier to store the\n",
"variable as a class attribute `self.lines`. Since the class is essentially a\n",
"function, we implement the function as the special method `__call__`\n",
"such that the instance `plot_u(u, x, t, n)` can be called as a standard\n",
"callback function from `solver`.\n",
"\n",
"### Making movie files\n",
"\n",
"From the\n",
"`frame_*.png` files containing the frames in the animation we can\n",
"make video files.\n",
"% if BOOK == \"book\":\n",
"the section [vib:ode1:anim](#vib:ode1:anim) presents basic information on how to\n",
"use the `ffmpeg` program for producing video files\n",
"in different modern formats: Flash, MP4, Webm, and Ogg.\n",
"% else:\n",
"We use the `ffmpeg` program to combine individual\n",
"plot files to movies in modern formats: Flash, MP4, Webm, and Ogg.\n",
"A typical `ffmpeg` command for creating a movie file\n",
"in Ogg format\n",
"with 4 frames per second built from a collection of\n",
"plot files with names generated by `frame_%04d.png`,\n",
"look like"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Terminal> ffmpeg -r 25 -i frame_%04d.png -c:v libtheora movie.ogg\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The different formats require different video encoders (`-c:v`) to\n",
"be installed: Flash applies `flv`, WebM applies `libvpx`, and MP4\n",
"applies `libx264`:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Terminal> ffmpeg -r 25 -i frame_%04d.png -c:v flv movie.flv\n",
" Terminal> ffmpeg -r 25 -i frame_%04d.png -c:v libvpx movie.webm\n",
" Terminal> ffmpeg -r 25 -i frame_%04d.png -c:v libx264 movie.mp4\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Players like `vlc`, `mplayer`, `gxine`, and `totem`\n",
"can be used to play these movie files.\n",
"\n",
"Note that padding the frame counter with zeros in the `frame_*.png`\n",
"files, as specified by the `%04d` format, is essential so that the wildcard\n",
"notation `frame_*.png` expands to the correct set of files.\n",
"% endif\n",
"\n",
"\n",
"The `viz` function creates an `ffmpeg` command\n",
"with the proper arguments for each of the formats Flash, MP4, WebM,\n",
"and Ogg. The task is greatly simplified by having a\n",
"`codec2ext` dictionary for mapping\n",
"video codec names to filename extensions.\n",
"% if BOOK == \"book\":\n",
"As mentioned in the section [vib:ode1:anim](#vib:ode1:anim), only\n",
"% else:\n",
"Only\n",
"% endif\n",
"two formats are actually needed to ensure that all browsers can\n",
"successfully play the video: MP4 and WebM.\n",
"\n",
"Some animations having a large number of plot files may not\n",
"be properly combined into a video using `ffmpeg`.\n",
"A method that always works is to play the PNG files as an animation\n",
"in a browser using JavaScript code in an HTML file.\n",
"The SciTools package has a function `movie` (or a stand-alone command\n",
"`scitools movie`) for creating such an HTML player. The `plt.movie`\n",
"call in the `viz` function shows how the function is used.\n",
"The file `movie.html` can be loaded into a browser and features\n",
"a user interface where the speed of the animation can be controlled.\n",
"Note that the movie in this case consists of the `movie.html` file\n",
"and all the frame files `tmp_*.png`.\n",
"\n",
"\n",
"### Skipping frames for animation speed\n",
"\n",
"Sometimes the time step is small and $T$ is large, leading to an\n",
"inconveniently large number of plot files and a slow animation on the\n",
"screen. The solution to such a problem is to decide on a total number\n",
"of frames in the animation, `num_frames`, and plot the solution only for\n",
"every `skip_frame` frames. For example, setting `skip_frame=5` leads\n",
"to plots of every 5 frames. The default value `skip_frame=1` plots\n",
"every frame.\n",
"The total number of time levels (i.e., maximum\n",
"possible number of frames) is the length of `t`, `t.size` (or `len(t)`),\n",
"so if we want `num_frames` frames in the animation,\n",
"we need to plot every `t.size/num_frames` frames:"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [],
"source": [
"def plot_u(t):\n",
" skip_frame = int(t.size/float(num_frames))\n",
" if n % skip_frame == 0 or n == t.size-1:\n",
" st.plot(x, u, 'r-', ...)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The initial condition (`n=0`) is included by `n % skip_frame == 0`,\n",
"as well as every `skip_frame`-th frame.\n",
"As `n % skip_frame == 0` will very seldom be true for the\n",
"very final frame, we must also check if `n == t.size-1` to\n",
"get the final frame included.\n",
"\n",
"A simple choice of numbers may illustrate the formulas: say we have\n",
"801 frames in total (`t.size`) and we allow only 60 frames to be\n",
"plotted. As `n` then runs from 801 to 0, we need to plot every 801/60\n",
"frame, which with integer division yields 13 as `skip_frame`. Using\n",
"the mod function, `n % skip_frame`, this operation is zero every time\n",
"`n` can be divided by 13 without a remainder. That is, the `if` test\n",
"is true when `n` equals $0, 13, 26, 39, ..., 780, 801$. The associated\n",
"code is included in the `plot_u` function, inside the `viz` function,\n",
"in the file [`wave1D_u0.py`](${src_wave}/wave1D/wave1D_u0.py).\n",
"\n",
"## Running a case\n",
"\n",
"\n",
"The first demo of our 1D wave equation solver concerns vibrations of a\n",
"string that is initially deformed to a triangular shape, like when picking\n",
"a guitar string:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"I(x) = \\left\\lbrace\n",
"\\begin{array}{ll}\n",
"ax/x_0, & x < x_0,\\\\ \n",
"a(L-x)/(L-x_0), & \\hbox{otherwise}\n",
"\\end{array}\\right.\n",
"\\label{wave:pde1:guitar:I} \\tag{1}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We choose $L=75$ cm, $x_0=0.8L$, $a=5$ mm, and a time frequency\n",
"$\\nu = 440$ Hz. The relation between the wave speed $c$ and $\\nu$ is\n",
"$c=\\nu\\lambda$, where $\\lambda$ is the wavelength, taken as $2L$ because\n",
"the longest wave on the string forms frac{1}{2} a wavelength. There is no\n",
"external force, so $f=0$ (meaning we can neglect gravity),\n",
"and the string is at rest initially, implying $V=0$.\n",
"\n",
"Regarding numerical parameters, we need to specify a $\\Delta t$.\n",
"Sometimes it is more natural to think of a spatial resolution instead\n",
"of a time step. A natural semi-coarse spatial resolution in the present\n",
"problem is $N_x=50$. We can then choose the associated $\\Delta t$ (as required\n",
"by the `viz` and `solver` functions) as the stability limit:\n",
"$\\Delta t = L/(N_xc)$. This is the $\\Delta t$ to be specified,\n",
"but notice that if $C<1$, the actual $\\Delta x$ computed in `solver` gets\n",
"larger than $L/N_x$: $\\Delta x = c\\Delta t/C = L/(N_xC)$. (The reason\n",
"is that we fix $\\Delta t$ and adjust $\\Delta x$, so if $C$ gets\n",
"smaller, the code implements this effect in terms of a larger $\\Delta x$.)\n",
"\n",
"A function for setting the physical and numerical parameters and\n",
"calling `viz` in this application goes as follows:"
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [],
"source": [
"def guitar(C):\n",
" \"\"\"Triangular wave (pulled guitar string).\"\"\"\n",
" L = 0.75\n",
" x0 = 0.8*L\n",
" a = 0.005\n",
" freq = 440\n",
" wavelength = 2*L\n",
" c = freq*wavelength\n",
" omega = 2*pi*freq\n",
" num_periods = 1\n",
" T = 2*pi/omega*num_periods\n",
" # Choose dt the same as the stability limit for Nx=50\n",
" dt = L/50./c\n",
"\n",
" def I(x):\n",
" return a*x/x0 if x < x0 else a/(L-x0)*(L-x)\n",
"\n",
" umin = -1.2*a; umax = -umin\n",
" cpu = viz(I, 0, 0, c, L, dt, C, T, umin, umax,\n",
" animate=True, tool='scitools')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The associated program has the name [`wave1D_u0.py`](${src_wave}/wave1D/wave1D_u0.py). Run\n",
"the program and watch the [movie of the vibrating string](${doc_notes}/wave/html/mov-wave/guitar_C0.8/movie.html).\n",
"The string should ideally consist of straight segments, but these are\n",
"somewhat wavy due to numerical approximation. Run the case with the\n",
"`wave1D_u0.py` code and $C=1$ to see the exact solution.\n",
"\n",
"\n",
"## Working with a scaled PDE model\n",
"\n",
"\n",
"Depending on the model, it may be a substantial job to establish\n",
"consistent and relevant physical parameter values for a case. The\n",
"guitar string example illustrates the point. However, by *scaling*\n",
"the mathematical problem we can often reduce the need to estimate\n",
"physical parameters dramatically. The scaling technique consists of\n",
"introducing new independent and dependent variables, with the aim that\n",
"the absolute values of these lie in $[0,1]$. We introduce the\n",
"dimensionless variables (details are found in Section 3.1.1 in [[Langtangen_scaling]](#Langtangen_scaling))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\bar x = \\frac{x}{L},\\quad \\bar t = \\frac{c}{L}t,\\quad\n",
"\\bar u = \\frac{u}{a}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, $L$ is a typical length scale, e.g., the length of the domain,\n",
"and $a$ is a typical size of $u$, e.g., determined from the\n",
"initial condition: $a=\\max_x|I(x)|$.\n",
"\n",
"We get by the chain rule that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial u}{\\partial t} =\n",
"\\frac{\\partial}{\\partial\\bar t}\\left(a\\bar u\\right)\n",
"\\frac{d\\bar t}{dt} =\n",
"\\frac{ac}{L}\\frac{\\partial\\bar u}{\\partial\\bar t}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial u}{\\partial x}\n",
"= \\frac{a}{L}\\frac{\\partial\\bar u}{\\partial\\bar x}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Inserting the dimensionless variables in the PDE gives, in case $f=0$,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{a^2c^2}{L^2}\\frac{\\partial^2\\bar u}{\\partial\\bar t^2}\n",
"= \\frac{a^2c^2}{L^2}\\frac{\\partial^2\\bar u}{\\partial\\bar x^2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Dropping the bars, we arrive at the scaled PDE"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial^2 u}{\\partial t^2} = \\frac{\\partial^2 u}{\\partial x^2},\n",
"\\quad x\\in (0,1),\\ t\\in (0,cT/L),\n",
"\\label{_auto1} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which has no parameter $c^2$ anymore. The initial conditions are scaled\n",
"as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a\\bar u(\\bar x, 0) = I(L\\bar x)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{a}{L/c}\\frac{\\partial\\bar u}{\\partial\\bar t}(\\bar x,0) = V(L\\bar x),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\bar u(\\bar x, 0) = \\frac{I(L\\bar x)}{\\max_x |I(x)|},\\quad\n",
"\\frac{\\partial\\bar u}{\\partial\\bar t}(\\bar x,0) = \\frac{L}{ac}V(L\\bar x)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the common case $V=0$ we see that there are no physical parameters to be\n",
"estimated in the PDE model!\n",
"\n",
"If we have a program implemented for the physical wave equation with\n",
"dimensions, we can obtain the dimensionless, scaled version by\n",
"setting $c=1$. The initial condition of a guitar string,\n",
"given in ([1](#wave:pde1:guitar:I)), gets its scaled form by choosing\n",
"$a=1$, $L=1$, and $x_0\\in [0,1]$. This means that we only need to\n",
"decide on the $x_0$ value as a fraction of unity, because\n",
"the scaled problem corresponds to setting all\n",
"other parameters to unity. In the code we can just set\n",
"`a=c=L=1`, `x0=0.8`, and there is no need to calculate with\n",
"wavelengths and frequencies to estimate $c$!\n",
"\n",
"The only non-trivial parameter to estimate in the scaled problem\n",
"is the final end time of the simulation, or more precisely, how it relates\n",
"to periods in periodic solutions in time, since we often want to\n",
"express the end time as a certain number of periods.\n",
"The period in the dimensionless problem is 2, so the end time can be\n",
"set to the desired number of periods times 2.\n",
"\n",
"Why the dimensionless period is 2 can be explained by the following\n",
"reasoning.\n",
"Suppose that $u$ behaves as $\\cos (\\omega t)$ in time in the original\n",
"problem with dimensions. The corresponding period is then $P=2\\pi/\\omega$, but\n",
"we need to estimate $\\omega$. A typical solution of the wave\n",
"equation is $u(x,t)=A\\cos(kx)\\cos(\\omega t)$, where $A$ is an amplitude\n",
"and $k$ is related to the wave length $\\lambda$ in space: $\\lambda = 2\\pi/k$.\n",
"Both $\\lambda$ and $A$ will be given by the initial condition $I(x)$.\n",
"Inserting this $u(x,t)$ in the PDE yields $-\\omega^2 = -c^2k^2$, i.e.,\n",
"$\\omega = kc$. The period is therefore $P=2\\pi/(kc)$.\n",
"If the boundary conditions are $u(0,t)=u(L,t)$, we need to have\n",
"$kL = n\\pi$ for integer $n$. The period becomes $P=2L/nc$. The longest\n",
"period is $P=2L/c$. The dimensionless period $\\tilde P$ is obtained\n",
"by dividing $P$ by the time scale $L/c$, which results in $\\tilde P=2$.\n",
"Shorter waves in the initial condition will have a dimensionless\n",
"shorter period $\\tilde P=2/n$ ($n>1$).\n",
"\n",
"# Vectorization\n",
"\n",
"\n",
"\n",
"The computational algorithm for solving the wave equation visits one\n",
"mesh point at a time and evaluates a formula for the new value\n",
"$u_i^{n+1}$ at that point. Technically, this is implemented by a loop\n",
"over array elements in a program. Such loops may run slowly in Python\n",
"(and similar interpreted languages such as R and MATLAB). One\n",
"technique for speeding up loops is to perform operations on entire\n",
"arrays instead of working with one element at a time. This is referred\n",
"to as *vectorization*, *vector computing*, or *array computing*.\n",
"Operations on whole arrays are possible if the computations involving\n",
"each element is independent of each other and therefore can, at least\n",
"in principle, be performed simultaneously. That is, vectorization not\n",
"only speeds up the code on serial computers, but also makes it easy to\n",
"exploit parallel computing. Actually, there are Python tools like\n",
"[Numba](http://numba.pydata.org) that can automatically turn\n",
"vectorized code into parallel code.\n",
"\n",
"\n",
"## Operations on slices of arrays\n",
"\n",
"\n",
"\n",
"Efficient computing with `numpy` arrays demands that we avoid loops\n",
"and compute with entire arrays at once (or at least large portions of them).\n",
"Consider this calculation of differences $d_i = u_{i+1}-u_i$:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"def differences(u):\n",
" n = u.size\n",
" d = np.zeros(n-1)\n",
" for i in range(0, n-1):\n",
" d[i] = u[i+1] - u[i]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All the differences here are independent of each other.\n",
"The computation of `d` can therefore alternatively be done by\n",
"subtracting the array $(u_0,u_1,\\ldots,u_{n-1})$ from\n",
"the array where the elements are shifted one index upwards:\n",
"$(u_1,u_2,\\ldots,u_n)$, see [Figure](#wave:pde1:vec:fig1).\n",
"The former subset of the array can be\n",
"expressed by `u[0:n-1]`,\n",
"`u[0:-1]`, or just\n",
"`u[:-1]`, meaning from index 0 up to,\n",
"but not including, the last element (`-1`). The latter subset\n",
"is obtained by `u[1:n]` or `u[1:]`,\n",
"meaning from index 1 and the rest of the array.\n",
"The computation of `d` can now be done without an explicit Python loop:\n",
"\n",
"`d = u[1:] - u[:-1]`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or with explicit limits if desired:\n",
"\n",
"`d = u[1:n] - u[0:n-1]`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indices with a colon, going from an index to (but not including) another\n",
"index are called *slices*. With `numpy` arrays, the computations\n",
"are still done by loops, but in efficient, compiled, highly optimized\n",
"C or Fortran code. Such loops are sometimes referred to as *vectorized\n",
"loops*. Such loops can also easily be distributed\n",
"among many processors on parallel computers. We say that the *scalar code*\n",
"above, working on an element (a scalar) at a time, has been replaced by\n",
"an equivalent *vectorized code*. The process of vectorizing code is called\n",
"*vectorization*.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"Illustration of subtracting two slices of two arrays.
\n",
"![](\"fig-wave/vectorized_diff.png\")
\n",
"\n",
"\n",
"\n",
"\n",
"**Test your understanding.**\n",
"\n",
"Newcomers to vectorization are encouraged to choose\n",
"a small array `u`, say with five elements,\n",
"and simulate with pen and paper\n",
"both the loop version and the vectorized version above.\n",
"\n",
"\n",
"\n",
"Finite difference schemes basically contain differences between array\n",
"elements with shifted indices. As an example,\n",
"consider the updating formula"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'n' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mu2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'n' is not defined"
]
}
],
"source": [
"for i in range(1, n-1):\n",
" u2[i] = u[i-1] - 2*u[i] + u[i+1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The vectorization consists of replacing the loop by arithmetics on\n",
"slices of arrays of length `n-2`:"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "unsupported operand type(s) for -: 'list' and 'list'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mu2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mu2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;31m# alternative\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for -: 'list' and 'list'"
]
}
],
"source": [
"u2 = u[:-2] - 2*u[1:-1] + u[2:]\n",
"u2 = u[0:n-2] - 2*u[1:n-1] + u[2:n] # alternative"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the length of `u2` becomes `n-2`. If `u2` is already an array of\n",
"length `n` and we want to use the formula to update all the \"inner\"\n",
"elements of `u2`, as we will when solving a 1D wave equation, we can write"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "unsupported operand type(s) for -: 'list' and 'list'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mu2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mu2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;31m# alternative\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for -: 'list' and 'list'"
]
}
],
"source": [
"u2[1:-1] = u[:-2] - 2*u[1:-1] + u[2:]\n",
"u2[1:n-1] = u[0:n-2] - 2*u[1:n-1] + u[2:n] # alternative"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first expression's right-hand side is realized by the\n",
"following steps, involving temporary arrays with intermediate results,\n",
"since each array operation can only involve one or two arrays.\n",
"The `numpy` package performs (behind the scenes) the first line above in\n",
"four steps:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" temp1 = 2*u[1:-1]\n",
" temp2 = u[:-2] - temp1\n",
" temp3 = temp2 + u[2:]\n",
" u2[1:-1] = temp3\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need three temporary arrays, but a user does not need to worry about\n",
"such temporary arrays.\n",
"\n",
"**Common mistakes with array slices.**\n",
"\n",
"Array expressions with slices demand that the slices have the same\n",
"shape. mathcal{I}_t easy to make a mistake in, e.g.,"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'n' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mu2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mNameError\u001b[0m: name 'n' is not defined"
]
}
],
"source": [
"u2[1:n-1] = u[0:n-2] - 2*u[1:n-1] + u[2:n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and write"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'n' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mu2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mNameError\u001b[0m: name 'n' is not defined"
]
}
],
"source": [
"u2[1:n-1] = u[0:n-2] - 2*u[1:n-1] + u[1:n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now `u[1:n]` has wrong length (`n-1`) compared to the other array\n",
"slices, causing a `ValueError` and the message\n",
"`could not broadcast input array from shape 103 into shape 104`\n",
"(if `n` is 105). When such errors occur one must closely examine\n",
"all the slices. Usually, it is easier to get upper limits of slices\n",
"right when they use `-1` or `-2` or empty limit rather than\n",
"expressions involving the length.\n",
"\n",
"Another common mistake, when `u2` has length `n`, is to forget the slice in the array on the\n",
"left-hand side,"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'n' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mu2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mNameError\u001b[0m: name 'n' is not defined"
]
}
],
"source": [
"u2 = u[0:n-2] - 2*u[1:n-1] + u[1:n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is really crucial: now `u2` becomes a *new* array of length\n",
"`n-2`, which is the wrong length as we have no entries for the boundary\n",
"values. We meant to insert the right-hand side array *into* the\n",
"original `u2` array for the entries that correspond to the\n",
"internal points in the mesh (`1:n-1` or `1:-1`).\n",
"\n",
"\n",
"\n",
"Vectorization may also work nicely with functions. To illustrate, we may\n",
"extend the previous example as follows:"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [
{
"ename": "IndexError",
"evalue": "list index out of range",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0mu2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mIndexError\u001b[0m: list index out of range"
]
}
],
"source": [
"def f(x):\n",
" return x**2 + 1\n",
"\n",
"n = 50\n",
"u2 = np.zeros(n+1)\n",
"\n",
"for i in range(1, n-1):\n",
" u2[i] = u[i-1] - 2*u[i] + u[i+1] + f(x[i])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assuming `u2`, `u`, and `x` all have length `n`, the vectorized\n",
"version becomes"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "unsupported operand type(s) for -: 'list' and 'list'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mu2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for -: 'list' and 'list'"
]
}
],
"source": [
"u2[1:-1] = u[:-2] - 2*u[1:-1] + u[2:] + f(x[1:-1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Obviously, `f` must be able to take an array as argument for `f(x[1:-1])`\n",
"to make sense.\n",
"\n",
"## Finite difference schemes expressed as slices\n",
"\n",
"\n",
"\n",
"We now have the necessary tools to vectorize the wave equation\n",
"algorithm as described mathematically in the section [wave:string:alg](#wave:string:alg)\n",
"and through code in the section [The solver function](#wave:pde1:impl:solver). There are\n",
"three loops: one for the initial condition, one for the first time\n",
"step, and finally the loop that is repeated for all subsequent time\n",
"levels. Since only the latter is repeated a potentially large number\n",
"of times, we limit our vectorization efforts to this loop. Within the time\n",
"loop, the space loop reads:"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'Nx' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu_n\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mu_nm1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mC2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu_n\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mu_n\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mu_n\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'Nx' is not defined"
]
}
],
"source": [
"for i in range(1, Nx):\n",
" u[i] = 2*u_n[i] - u_nm1[i] + \\\n",
" C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The vectorized version becomes"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"u[1:-1] = - u_nm1[1:-1] + 2*u_n[1:-1] + \\\n",
" C2*(u_n[:-2] - 2*u_n[1:-1] + u_n[2:])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"u[1:Nx] = 2*u_n[1:Nx]- u_nm1[1:Nx] + \\\n",
" C2*(u_n[0:Nx-1] - 2*u_n[1:Nx] + u_n[2:Nx+1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"The program\n",
"[`wave1D_u0v.py`](${src_wave}/wave1D/wave1D_u0v.py)\n",
"contains a new version of the function `solver` where both the scalar\n",
"and the vectorized loops are included (the argument `version` is\n",
"set to `scalar` or `vectorized`, respectively).\n",
"\n",
"\n",
"## Verification\n",
"\n",
"\n",
"\n",
"We may reuse the quadratic solution $\\uex(x,t)=x(L-x)(1+{\\frac{1}{2}}t)$ for\n",
"verifying also the vectorized code. A test function can now verify\n",
"both the scalar and the vectorized version. Moreover, we may\n",
"use a `user_action` function that compares the computed and exact\n",
"solution at each time level and performs a test:"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"def test_quadratic():\n",
" \"\"\"\n",
" Check the scalar and vectorized versions for\n",
" a quadratic u(x,t)=x(L-x)(1+t/2) that is exactly reproduced.\n",
" \"\"\"\n",
" # The following function must work for x as array or scalar\n",
" u_exact = lambda x, t: x*(L - x)*(1 + 0.5*t)\n",
" I = lambda x: u_exact(x, 0)\n",
" V = lambda x: 0.5*u_exact(x, 0)\n",
" # f is a scalar (zeros_like(x) works for scalar x too)\n",
" f = lambda x, t: np.zeros_like(x) + 2*c**2*(1 + 0.5*t)\n",
"\n",
" L = 2.5\n",
" c = 1.5\n",
" C = 0.75\n",
" Nx = 3 # Very coarse mesh for this exact test\n",
" dt = C*(L/Nx)/c\n",
" T = 18\n",
"\n",
" def assert_no_error(u, x, t, n):\n",
" u_e = u_exact(x, t[n])\n",
" tol = 1E-13\n",
" diff = np.abs(u - u_e).max()\n",
" assert diff < tol\n",
"\n",
" solver(I, V, f, c, L, dt, C, T,\n",
" user_action=assert_no_error, version='scalar')\n",
" solver(I, V, f, c, L, dt, C, T,\n",
" user_action=assert_no_error, version='vectorized')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Lambda functions.**\n",
"\n",
"The code segment above demonstrates how to achieve very\n",
"compact code, without degraded readability,\n",
"by use of lambda functions for the various\n",
"input parameters that require a Python function. In essence,"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"f = lambda x, t: L*(x-t)**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"is equivalent to"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"def f(x, t):\n",
" return L(x-t)**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that lambda functions can just contain a single expression and no\n",
"statements.\n",
"\n",
"One advantage with lambda functions is that they can be used directly\n",
"in calls:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"solver(I=lambda x: sin(pi*x/L), V=0, f=0, ...)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Efficiency measurements\n",
"\n",
"\n",
"The `wave1D_u0v.py` contains our new `solver` function with both\n",
"scalar and vectorized code. For comparing the efficiency\n",
"of scalar versus vectorized code, we need a `viz` function\n",
"as discussed in the section [Visualization: animating the solution](#wave:pde1:impl:animate).\n",
"All of this `viz` function can be reused, except the call\n",
"to `solver_function`. This call lacks the parameter\n",
"`version`, which we want to set to `vectorized` and `scalar`\n",
"for our efficiency measurements.\n",
"\n",
"One solution is to copy the `viz` code from `wave1D_u0` into\n",
"`wave1D_u0v.py` and add a `version` argument to the `solver_function` call.\n",
"Taking into account how much animation code we\n",
"then duplicate, this is not a good idea.\n",
"Alternatively,\n",
"introducing the `version` argument in `wave1D_u0.viz`, so that this function\n",
"can be imported into `wave1D_u0v.py`, is not\n",
"a good solution either, since `version` has no meaning in that file.\n",
"We need better ideas!\n",
"\n",
"### Solution 1\n",
"\n",
"Calling `viz` in `wave1D_u0` with `solver_function` as our new\n",
"solver in `wave1D_u0v` works fine, since this solver has\n",
"`version='vectorized'` as default value. The problem arises when we\n",
"want to test `version='scalar'`. The simplest solution is then\n",
"to use `wave1D_u0.solver` instead. We make a new `viz` function\n",
"in `wave1D_u0v.py` that has a `version` argument and that just\n",
"calls `wave1D_u0.viz`:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"def viz(\n",
" I, V, f, c, L, dt, C, T, # PDE parameters\n",
" umin, umax, # Interval for u in plots\n",
" animate=True, # Simulation with animation?\n",
" tool='matplotlib', # 'matplotlib' or 'scitools'\n",
" solver_function=solver, # Function with numerical algorithm\n",
" version='vectorized', # 'scalar' or 'vectorized'\n",
" ):\n",
" import wave1D_u0\n",
" if version == 'vectorized':\n",
" # Reuse viz from wave1D_u0, but with the present\n",
" # modules' new vectorized solver (which has\n",
" # version='vectorized' as default argument;\n",
" # wave1D_u0.viz does not feature this argument)\n",
" cpu = wave1D_u0.viz(\n",
" I, V, f, c, L, dt, C, T, umin, umax,\n",
" animate, tool, solver_function=solver)\n",
" elif version == 'scalar':\n",
" # Call wave1D_u0.viz with a solver with\n",
" # scalar code and use wave1D_u0.solver.\n",
" cpu = wave1D_u0.viz(\n",
" I, V, f, c, L, dt, C, T, umin, umax,\n",
" animate, tool,\n",
" solver_function=wave1D_u0.solver)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solution 2\n",
"\n",
"There is a more advanced and fancier solution featuring a very useful trick:\n",
"we can make a new function that will always call `wave1D_u0v.solver`\n",
"with `version='scalar'`. The `functools.partial` function from\n",
"standard Python takes a function `func` as argument and\n",
"a series of positional and keyword arguments and returns a\n",
"new function that will call `func` with the supplied arguments,\n",
"while the user can control all the other arguments in `func`.\n",
"Consider a trivial example,"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"def f(a, b, c=2):\n",
" return a + b + c"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to ensure that `f` is always called with `c=3`, i.e., `f`\n",
"has only two \"free\" arguments `a` and `b`.\n",
"This functionality is obtained by"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"import functools\n",
"f2 = functools.partial(f, c=3)\n",
"\n",
"print f2(1, 2) # results in 1+2+3=6"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now `f2` calls `f` with whatever the user supplies as `a` and `b`,\n",
"but `c` is always `3`.\n",
"\n",
"Back to our `viz` code, we can do"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"import functools\n",
"# Call wave1D_u0.solver with version fixed to scalar\n",
"scalar_solver = functools.partial(wave1D_u0.solver, version='scalar')\n",
"cpu = wave1D_u0.viz(\n",
" I, V, f, c, L, dt, C, T, umin, umax,\n",
" animate, tool, solver_function=scalar_solver)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The new `scalar_solver` takes the same arguments as\n",
"`wave1D_u0.scalar` and calls `wave1D_u0v.scalar`,\n",
"but always supplies the extra argument\n",
"`version='scalar'`. When sending this `solver_function`\n",
"to `wave1D_u0.viz`, the latter will call `wave1D_u0v.solver`\n",
"with all the `I`, `V`, `f`, etc., arguments we supply, plus\n",
"`version='scalar'`.\n",
"\n",
"### Efficiency experiments\n",
"\n",
"We now have a `viz` function that can call our solver function both in\n",
"scalar and vectorized mode. The function `run_efficiency_experiments`\n",
"in `wave1D_u0v.py` performs a set of experiments and reports the\n",
"CPU time spent in the scalar and vectorized solver for\n",
"the previous string vibration example with spatial mesh resolutions\n",
"$N_x=50,100,200,400,800$. Running this function reveals\n",
"that the vectorized\n",
"code runs substantially faster: the vectorized code runs approximately\n",
"$N_x/10$ times as fast as the scalar code!\n",
"\n",
"## Remark on the updating of arrays\n",
"\n",
"\n",
"\n",
"At the end of each time step we need to update the `u_nm1` and `u_n`\n",
"arrays such that they have the right content for the next time step:"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"u_nm1[:] = u_n\n",
"u_n[:] = u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The order here is important: updating `u_n` first, makes `u_nm1` equal\n",
"to `u`, which is wrong!\n",
"\n",
"The assignment `u_n[:] = u` copies the content of the `u` array into\n",
"the elements of the `u_n` array. Such copying takes time, but\n",
"that time is negligible compared to the time needed for\n",
"computing `u` from the finite difference formula,\n",
"even when the formula has a vectorized implementation.\n",
"However, efficiency of program code is a key topic when solving\n",
"PDEs numerically (particularly when there are two or three\n",
"space dimensions), so it must be mentioned that there exists a\n",
"much more efficient way of making the arrays `u_nm1` and `u_n`\n",
"ready for the next time step. The idea is based on *switching\n",
"references* and explained as follows.\n",
"\n",
"A Python variable is actually a reference to some object (C programmers\n",
"may think of pointers). Instead of copying data, we can let `u_nm1`\n",
"refer to the `u_n` object and `u_n` refer to the `u` object.\n",
"This is a very efficient operation (like switching pointers in C).\n",
"A naive implementation like"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"u_nm1 = u_n\n",
"u_n = u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"will fail, however, because now `u_nm1` refers to the `u_n` object,\n",
"but then the name `u_n` refers to `u`, so that this `u` object\n",
"has two references, `u_n` and `u`, while our third array, originally\n",
"referred to by `u_nm1`, has no more references and is lost.\n",
"This means that the variables `u`, `u_n`, and `u_nm1` refer to two\n",
"arrays and not three. Consequently, the computations at the next\n",
"time level will be messed up, since updating the elements in\n",
"`u` will imply updating the elements in `u_n` too, thereby destroying\n",
"the solution at the previous time step.\n",
"\n",
"While `u_nm1 = u_n` is fine, `u_n = u` is problematic, so\n",
"the solution to this problem is to ensure that `u`\n",
"points to the `u_nm1` array. This is mathematically wrong, but\n",
"new correct values will be filled into `u` at the next time step\n",
"and make it right.\n",
"\n",
"The correct switch of references is"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"tmp = u_nm1\n",
"u_nm1 = u_n\n",
"u_n = u\n",
"u = tmp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can get rid of the temporary reference `tmp` by writing"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"u_nm1, u_n, u = u_n, u, u_nm1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This switching of references for updating our arrays\n",
"will be used in later implementations.\n",
"\n",
"**Caution:**\n",
"\n",
"The update `u_nm1, u_n, u = u_n, u, u_nm1` leaves wrong content in `u`\n",
"at the final time step. This means that if we return `u`, as we\n",
"do in the example codes here, we actually return `u_nm1`, which is\n",
"obviously wrong. mathcal{I}_t is therefore important to adjust the content\n",
"of `u` to `u = u_n` before returning `u`. (Note that\n",
"the `user_action` function\n",
"reduces the need to return the solution from the solver.)\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"# Exercises\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 1: Simulate a standing wave\n",
"\n",
"\n",
"The purpose of this exercise is to simulate standing waves on $[0,L]$\n",
"and illustrate the error in the simulation.\n",
"Standing waves arise from an initial condition"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x,0)= A \\sin\\left(\\frac{\\pi}{L}mx\\right),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $m$ is an integer and $A$ is a freely chosen amplitude.\n",
"The corresponding exact solution can be computed and reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\uex(x,t) = A\\sin\\left(\\frac{\\pi}{L}mx\\right)\n",
"\\cos\\left(\\frac{\\pi}{L}mct\\right)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**a)**\n",
"Explain that for a function $\\sin kx\\cos \\omega t$ the wave length\n",
"in space is $\\lambda = 2\\pi /k$ and the period in time is $P=2\\pi/\\omega$.\n",
"Use these expressions to find the wave length in space and period in\n",
"time of $\\uex$ above.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Since the sin and cos functions depend on $x$ and $t$, respectively,\n",
"the sin function will run through one period as $x$ increases by $\\frac{2\\pi}{k}$, while the cos function starts repeating as $t$ increases by $\\frac{2\\pi}{\\omega}$.\n",
"\n",
"The wave length in space becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lambda = \\frac{2\\pi}{\\frac{\\pi}{L}m} = \\frac{2L}{m}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The period in time becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"P = \\frac{2\\pi}{\\frac{\\pi}{L}mc} = \\frac{2L}{mc}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**b)**\n",
"Import the `solver` function from `wave1D_u0.py` into a new file\n",
"where the `viz` function is reimplemented such that it\n",
"plots either the numerical *and* the exact solution, *or* the error.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"See code below.\n",
"\n",
"\n",
"\n",
"**c)**\n",
"Make animations where you illustrate how the error\n",
"$e^n_i =\\uex(x_i, t_n)- u^n_i$\n",
"develops and increases in time. Also make animations of\n",
"$u$ and $\\uex$ simultaneously.\n",
"\n",
"\n",
"\n",
"**Hint 1.**\n",
"Quite long time simulations are needed in order to display significant\n",
"discrepancies between the numerical and exact solution.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Hint 2.**\n",
"A possible set of parameters is $L=12$, $m=9$, $c=2$, $A=1$, $N_x=80$,\n",
"$C=0.8$. The error mesh function $e^n$ can be simulated for 10 periods,\n",
"while 20-30 periods are needed to show significant differences between\n",
"the curves for the numerical and exact solution.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The code:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"import sys, os\n",
"sys.path.insert(0, os.path.join(os.pardir, os.pardir, \n",
" 'src-wave', 'wave1D'))\n",
"\n",
"from wave1D_u0 import solver\n",
"#from wave1D_u0v import solver # allows faster vectorized operations\n",
"import numpy as np\n",
"\n",
"def viz(I, V, f, c, L, dt, C, T,\n",
" ymax, # y axis: [-ymax, ymax]\n",
" u_exact, # u_exact(x, t)\n",
" animate='u and u_exact', # or 'error'\n",
" movie_filename='movie',\n",
" ):\n",
" \"\"\"Run solver and visualize u at each time level.\"\"\"\n",
" import scitools.std as plt\n",
" import time, glob, os\n",
"\n",
" class Plot:\n",
" def __init__(self, ymax, frame_name='frame'):\n",
" self.max_error = [] # hold max amplitude errors\n",
" self.max_error_t = [] # time points corresp. to max_error\n",
" self.frame_name = frame_name\n",
" self.ymax = ymax\n",
"\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" if animate == 'u and u_exact':\n",
" plt.plot(x, u, 'r-',\n",
" x, u_exact(x, t[n]), 'b--',\n",
" xlabel='x', ylabel='u',\n",
" axis=[0, L, -self.ymax, self.ymax],\n",
" title='t=%f' % t[n], show=True)\n",
" else:\n",
" error = u_exact(x, t[n]) - u\n",
" local_max_error = np.abs(error).max()\n",
" # self.max_error holds the increasing amplitude error\n",
" if self.max_error == [] or \\\n",
" local_max_error > max(self.max_error):\n",
" self.max_error.append(local_max_error)\n",
" self.max_error_t.append(t[n])\n",
" # Use user's ymax until the error exceeds that value.\n",
" # This gives a growing max value of the yaxis (but\n",
" # never shrinking)\n",
" self.ymax = max(self.ymax, max(self.max_error))\n",
" plt.plot(x, error, 'r-',\n",
" xlabel='x', ylabel='error',\n",
" axis=[0, L, -self.ymax, self.ymax],\n",
" title='t=%f' % t[n], show=True)\n",
" plt.savefig('%s_%04d.png' % (self.frame_name, n))\n",
"\n",
" # Clean up old movie frames\n",
" for filename in glob.glob('frame_*.png'):\n",
" os.remove(filename)\n",
"\n",
" plot = Plot(ymax)\n",
" u, x, t, cpu = solver(I, V, f, c, L, dt, C, T, plot)\n",
"\n",
" # Make plot of max error versus time\n",
" plt.figure()\n",
" plt.plot(plot.max_error_t, plot.max_error)\n",
" plt.xlabel('time'); plt.ylabel('max abs(error)')\n",
" plt.savefig('error.png')\n",
" plt.savefig('error.pdf')\n",
"\n",
" # Make .flv movie file\n",
" fps = 4 # Frames per second\n",
" codec2ext = dict(flv='flv') #, libx64='mp4', \n",
" #libvpx='webm', libtheora='ogg')\n",
"\n",
" filespec = 'frame_%04d.png'\n",
" movie_program = 'ffmpeg'\n",
" for codec in codec2ext:\n",
" ext = codec2ext[codec]\n",
" cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\\\n",
" '-vcodec %(codec)s %(movie_filename)s.%(ext)s' % vars()\n",
" os.system(cmd)\n",
"\n",
"def simulations():\n",
" from numpy import sin, cos, pi\n",
" L = 12 # length of domain\n",
" m = 9 # 2L/m: wave length or period in space (2*pi/k, k=pi*m/L)\n",
" c = 2 # wave velocity\n",
" A = 1 # amplitude\n",
" Nx = 80\n",
" C = 0.8\n",
" P = 2*pi/(pi*m*c/L) # 1 period in time\n",
" #T = 10*P\n",
" # Choose dt the same as the stability limit for Nx=50\n",
" dt = L/50./c \n",
"\n",
" def u_exact(x, t):\n",
" return A*sin(pi*m*x/L)*cos(pi*m*c*t/L)\n",
"\n",
" def I(x):\n",
" return u_exact(x, 0)\n",
"\n",
" V = 0\n",
" f = 0\n",
"\n",
" viz(I, V, f, c, L, dt, C, 10.5*P,\n",
" 0.1, u_exact,\n",
" animate='error',\n",
" movie_filename='error')\n",
"\n",
" # Very long simulation to demonstrate different curves\n",
" viz(I, V, f, c, L, dt, C, 30*P,\n",
" 1.2*A, u_exact,\n",
" animate='u and u_exact',\n",
" movie_filename='solution')\n",
"\n",
"if __name__ == '__main__':\n",
" simulations()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"Filename: `wave_standing`.\n",
"\n",
"\n",
"\n",
"### Remarks\n",
"\n",
"The important\n",
"parameters for numerical quality are $C$ and $k\\Delta x$, where\n",
"$C=c\\Delta t/\\Delta x$ is the Courant number and $k$ is defined above\n",
"($k\\Delta x$ is proportional to how many mesh points we have per wave length\n",
"in space, see the section [wave:pde1:num:dispersion](#wave:pde1:num:dispersion) for explanation).\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 2: Add storage of solution in a user action function\n",
"\n",
"\n",
"Extend the `plot_u` function in the file `wave1D_u0.py` to also store\n",
"the solutions `u` in a list.\n",
"To this end, declare `all_u` as\n",
"an empty list in the `viz` function, outside `plot_u`, and perform\n",
"an append operation inside the `plot_u` function. Note that a\n",
"function, like `plot_u`, inside another function, like `viz`,\n",
"remembers all local variables in `viz` function, including `all_u`,\n",
"even when `plot_u` is called (as `user_action`) in the `solver` function.\n",
"Test both `all_u.append(u)` and `all_u.append(u.copy())`.\n",
"Why does one of these constructions fail to store the solution correctly?\n",
"Let the `viz` function return the `all_u` list\n",
"converted to a two-dimensional `numpy` array.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We have to explicitly use a copy of u, i.e. as `all_u.append(u.copy())`, otherwise we just get a reference to `u`, which goes on changing with the computations."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"#!/usr/bin/env python\n",
"\"\"\"\n",
"1D wave equation with u=0 at the boundary.\n",
"Simplest possible implementation.\n",
"\n",
"The key function is::\n",
"\n",
" u, x, t, cpu = solver(I, V, f, c, L, dt, C, T, user_action)\n",
"\n",
"which solves the wave equation u_tt = c**2*u_xx on (0,L) with u=0\n",
"on x=0,L, for t in (0,T]. Initial conditions: u=I(x), u_t=V(x).\n",
"\n",
"T is the stop time for the simulation.\n",
"dt is the desired time step.\n",
"C is the Courant number (=c*dt/dx), which specifies dx.\n",
"f(x,t) is a function for the source term (can be 0 or None).\n",
"I and V are functions of x.\n",
"\n",
"user_action is a function of (u, x, t, n) where the calling\n",
"code can add visualization, error computations, etc.\n",
"\"\"\"\n",
"\n",
"import numpy as np\n",
"\n",
"def solver(I, V, f, c, L, dt, C, T, user_action=None):\n",
" \"\"\"Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\"\"\"\n",
" Nt = int(round(T/dt))\n",
" t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
" dx = dt*c/float(C)\n",
" Nx = int(round(L/dx))\n",
" x = np.linspace(0, L, Nx+1) # Mesh points in space\n",
" C2 = C**2 # Help variable in the scheme\n",
" # Make sure dx and dt are compatible with x and t\n",
" dx = x[1] - x[0]\n",
" dt = t[1] - t[0]\n",
"\n",
" if f is None or f == 0 :\n",
" f = lambda x, t: 0\n",
" if V is None or V == 0:\n",
" V = lambda x: 0\n",
"\n",
" u = np.zeros(Nx+1) # Solution array at new time level\n",
" u_1 = np.zeros(Nx+1) # Solution at 1 time level back\n",
" u_2 = np.zeros(Nx+1) # Solution at 2 time levels back\n",
"\n",
" import time; t0 = time.clock() # for measuring CPU time\n",
"\n",
" # Load initial condition into u_1\n",
" for i in range(0,Nx+1):\n",
" u_1[i] = I(x[i])\n",
"\n",
" if user_action is not None:\n",
" user_action(u_1, x, t, 0)\n",
"\n",
" # Special formula for first time step\n",
" n = 0\n",
" for i in range(1, Nx):\n",
" u[i] = u_1[i] + dt*V(x[i]) + \\\n",
" 0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" 0.5*dt**2*f(x[i], t[n])\n",
" u[0] = 0; u[Nx] = 0\n",
"\n",
" if user_action is not None:\n",
" user_action(u, x, t, 1)\n",
"\n",
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
" for n in range(1, Nt):\n",
" # Update all inner points at time t[n+1]\n",
" for i in range(1, Nx):\n",
" u[i] = - u_2[i] + 2*u_1[i] + \\\n",
" C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" dt**2*f(x[i], t[n])\n",
"\n",
" # Insert boundary conditions\n",
" u[0] = 0; u[Nx] = 0\n",
" if user_action is not None:\n",
" if user_action(u, x, t, n+1):\n",
" break\n",
"\n",
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
" cpu_time = t0 - time.clock()\n",
" return u, x, t, cpu_time\n",
"\n",
"def test_quadratic():\n",
" \"\"\"Check that u(x,t)=x(L-x)(1+t/2) is exactly reproduced.\"\"\"\n",
"\n",
" def u_exact(x, t):\n",
" return x*(L-x)*(1 + 0.5*t)\n",
"\n",
" def I(x):\n",
" return u_exact(x, 0)\n",
"\n",
" def V(x):\n",
" return 0.5*u_exact(x, 0)\n",
"\n",
" def f(x, t):\n",
" return 2*(1 + 0.5*t)*c**2\n",
"\n",
" L = 2.5\n",
" c = 1.5\n",
" C = 0.75\n",
" Nx = 6 # Very coarse mesh for this exact test\n",
" dt = C*(L/Nx)/c\n",
" T = 18\n",
"\n",
" def assert_no_error(u, x, t, n):\n",
" u_e = u_exact(x, t[n])\n",
" diff = np.abs(u - u_e).max()\n",
" tol = 1E-13\n",
" assert diff < tol\n",
"\n",
" solver(I, V, f, c, L, dt, C, T,\n",
" user_action=assert_no_error)\n",
"\n",
"def test_constant():\n",
" \"\"\"Check that u(x,t)=Q=0 is exactly reproduced.\"\"\"\n",
" u_const = 0 # Require 0 because of the boundary conditions\n",
" C = 0.75\n",
" dt = C # Very coarse mesh\n",
" u, x, t, cpu = solver(I=lambda x:\n",
" 0, V=0, f=0, c=1.5, L=2.5,\n",
" dt=dt, C=C, T=18)\n",
" tol = 1E-14\n",
" assert np.abs(u - u_const).max() < tol\n",
"\n",
"def viz(\n",
" I, V, f, c, L, dt, C, T, # PDE parameters\n",
" umin, umax, # Interval for u in plots\n",
" animate=True, # Simulation with animation?\n",
" tool='matplotlib', # 'matplotlib' or 'scitools'\n",
" solver_function=solver, # Function with numerical algorithm\n",
" ):\n",
" \"\"\"Run solver, store and visualize u at each time level.\"\"\"\n",
"\n",
" all_u = []\n",
" def plot_u_st(u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" plt.plot(x, u, 'r-',\n",
" xlabel='x', ylabel='u',\n",
" axis=[0, L, umin, umax],\n",
" title='t=%f' % t[n], show=True)\n",
" # Let the initial condition stay on the screen for 2\n",
" # seconds, else insert a pause of 0.2 s between each plot\n",
" time.sleep(2) if t[n] == 0 else time.sleep(0.2)\n",
" #plt.savefig('frame_%04d.png' % n) # for movie making\n",
" plt.savefig('tmp_%04d.png' % n) # for movie making \n",
" all_u.append(u.copy()) \n",
"\n",
" class PlotMatplotlib:\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" if n == 0:\n",
" plt.ion()\n",
" self.lines = plt.plot(x, u, 'r-')\n",
" plt.xlabel('x'); plt.ylabel('u')\n",
" plt.axis([0, L, umin, umax])\n",
" plt.legend(['t=%f' % t[n]], loc='lower left')\n",
" else:\n",
" self.lines[0].set_ydata(u)\n",
" plt.legend(['t=%f' % t[n]], loc='lower left')\n",
" plt.draw()\n",
" time.sleep(2) if t[n] == 0 else time.sleep(0.2)\n",
" plt.savefig('tmp_%04d.png' % n) # for movie making\n",
"\n",
" if tool == 'matplotlib':\n",
" import matplotlib.pyplot as plt\n",
" plot_u = PlotMatplotlib()\n",
" elif tool == 'scitools':\n",
" import scitools.std as plt # scitools.easyviz interface\n",
" plot_u = plot_u_st\n",
" import time, glob, os\n",
"\n",
" # Clean up old movie frames\n",
" for filename in glob.glob('tmp_*.png'):\n",
" os.remove(filename)\n",
"\n",
" # Call solver and do the simulaton\n",
" user_action = plot_u if animate else None\n",
" u, x, t, cpu = solver_function(\n",
" I, V, f, c, L, dt, C, T, user_action)\n",
"\n",
" # Make video files\n",
" fps = 4 # frames per second\n",
" codec2ext = dict(flv='flv', libx264='mp4', libvpx='webm',\n",
" libtheora='ogg') # video formats\n",
" filespec = 'tmp_%04d.png'\n",
" movie_program = 'ffmpeg'\n",
" for codec in codec2ext:\n",
" ext = codec2ext[codec]\n",
" cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\\\n",
" '-vcodec %(codec)s movie.%(ext)s' % vars()\n",
" os.system(cmd)\n",
"\n",
" if tool == 'scitools':\n",
" # Make an HTML play for showing the animation in a browser\n",
" plt.movie('tmp_*.png', encoder='html', fps=fps,\n",
" output_file='movie.html')\n",
" return cpu, np.array(all_u)\n",
"\n",
"def guitar(C):\n",
" \"\"\"Triangular wave (pulled guitar string).\"\"\"\n",
" L = 0.75\n",
" x0 = 0.8*L\n",
" a = 0.005\n",
" freq = 440\n",
" wavelength = 2*L\n",
" c = freq*wavelength\n",
" from math import pi\n",
" w = 2*pi*freq\n",
" num_periods = 1\n",
" T = 2*pi/w*num_periods\n",
" # Choose dt the same as the stability limit for Nx=50\n",
" dt = L/50./c\n",
"\n",
" def I(x):\n",
" return a*x/x0 if x < x0 else a/(L-x0)*(L-x)\n",
"\n",
" umin = -1.2*a; umax = -umin\n",
" cpu, all_u = viz(I, 0, 0, c, L, dt, C, T, umin, umax,\n",
" animate=True, tool='scitools')\n",
" # checking\n",
" #for e in all_u:\n",
" # print e[int(len(all_u[1])/2)]\n",
"\n",
"if __name__ == '__main__':\n",
" #test_quadratic()\n",
" import sys\n",
" try:\n",
" C = float(sys.argv[1])\n",
" print 'C=%g' % C\n",
" except IndexError:\n",
" C = 0.85\n",
" print 'Courant number: %.2f' % C\n",
" guitar(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Filename: `wave1D_u0_s_store`.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 3: Use a class for the user action function\n",
"\n",
"\n",
"Redo [Exercise 2: Add storage of solution in a user action function](#wave:exer:store:list) using a class for the user\n",
"action function. Let the `all_u` list be an attribute in this class\n",
"and implement the user action function as a method (the special method\n",
"`__call__` is a natural choice). The class versions avoid that the\n",
"user action function depends on parameters defined outside the\n",
"function (such as `all_u` in [Exercise 2: Add storage of solution in a user action function](#wave:exer:store:list)).\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Using a class, we get"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"#!/usr/bin/env python\n",
"\"\"\"\n",
"1D wave equation with u=0 at the boundary.\n",
"Simplest possible implementation.\n",
"\n",
"The key function is::\n",
"\n",
" u, x, t, cpu = solver(I, V, f, c, L, dt, C, T, user_action)\n",
"\n",
"which solves the wave equation u_tt = c**2*u_xx on (0,L) with u=0\n",
"on x=0,L, for t in (0,T]. Initial conditions: u=I(x), u_t=V(x).\n",
"\n",
"T is the stop time for the simulation.\n",
"dt is the desired time step.\n",
"C is the Courant number (=c*dt/dx), which specifies dx.\n",
"f(x,t) is a function for the source term (can be 0 or None).\n",
"I and V are functions of x.\n",
"\n",
"user_action is a function of (u, x, t, n) where the calling\n",
"code can add visualization, error computations, etc.\n",
"\"\"\"\n",
"\n",
"import numpy as np\n",
"\n",
"def solver(I, V, f, c, L, dt, C, T, user_action=None):\n",
" \"\"\"Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\"\"\"\n",
" Nt = int(round(T/dt))\n",
" t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
" dx = dt*c/float(C)\n",
" Nx = int(round(L/dx))\n",
" x = np.linspace(0, L, Nx+1) # Mesh points in space\n",
" C2 = C**2 # Help variable in the scheme\n",
" # Make sure dx and dt are compatible with x and t\n",
" dx = x[1] - x[0]\n",
" dt = t[1] - t[0]\n",
"\n",
" if f is None or f == 0 :\n",
" f = lambda x, t: 0\n",
" if V is None or V == 0:\n",
" V = lambda x: 0\n",
"\n",
" u = np.zeros(Nx+1) # Solution array at new time level\n",
" u_1 = np.zeros(Nx+1) # Solution at 1 time level back\n",
" u_2 = np.zeros(Nx+1) # Solution at 2 time levels back\n",
"\n",
" import time; t0 = time.clock() # for measuring CPU time\n",
"\n",
" # Load initial condition into u_1\n",
" for i in range(0,Nx+1):\n",
" u_1[i] = I(x[i])\n",
"\n",
" if user_action is not None:\n",
" user_action(u_1, x, t, 0)\n",
"\n",
" # Special formula for first time step\n",
" n = 0\n",
" for i in range(1, Nx):\n",
" u[i] = u_1[i] + dt*V(x[i]) + \\\n",
" 0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" 0.5*dt**2*f(x[i], t[n])\n",
" u[0] = 0; u[Nx] = 0\n",
"\n",
" if user_action is not None:\n",
" user_action(u, x, t, 1)\n",
"\n",
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
" for n in range(1, Nt):\n",
" # Update all inner points at time t[n+1]\n",
" for i in range(1, Nx):\n",
" u[i] = - u_2[i] + 2*u_1[i] + \\\n",
" C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" dt**2*f(x[i], t[n])\n",
"\n",
" # Insert boundary conditions\n",
" u[0] = 0; u[Nx] = 0\n",
" if user_action is not None:\n",
" if user_action(u, x, t, n+1):\n",
" break\n",
"\n",
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
" cpu_time = t0 - time.clock()\n",
" return u, x, t, cpu_time\n",
"\n",
"def test_quadratic():\n",
" \"\"\"Check that u(x,t)=x(L-x)(1+t/2) is exactly reproduced.\"\"\"\n",
"\n",
" def u_exact(x, t):\n",
" return x*(L-x)*(1 + 0.5*t)\n",
"\n",
" def I(x):\n",
" return u_exact(x, 0)\n",
"\n",
" def V(x):\n",
" return 0.5*u_exact(x, 0)\n",
"\n",
" def f(x, t):\n",
" return 2*(1 + 0.5*t)*c**2\n",
"\n",
" L = 2.5\n",
" c = 1.5\n",
" C = 0.75\n",
" Nx = 6 # Very coarse mesh for this exact test\n",
" dt = C*(L/Nx)/c\n",
" T = 18\n",
"\n",
" def assert_no_error(u, x, t, n):\n",
" u_e = u_exact(x, t[n])\n",
" diff = np.abs(u - u_e).max()\n",
" tol = 1E-13\n",
" assert diff < tol\n",
"\n",
" solver(I, V, f, c, L, dt, C, T,\n",
" user_action=assert_no_error)\n",
"\n",
"def test_constant():\n",
" \"\"\"Check that u(x,t)=Q=0 is exactly reproduced.\"\"\"\n",
" u_const = 0 # Require 0 because of the boundary conditions\n",
" C = 0.75\n",
" dt = C # Very coarse mesh\n",
" u, x, t, cpu = solver(I=lambda x:\n",
" 0, V=0, f=0, c=1.5, L=2.5,\n",
" dt=dt, C=C, T=18)\n",
" tol = 1E-14\n",
" assert np.abs(u - u_const).max() < tol\n",
"\n",
"def viz(\n",
" I, V, f, c, L, dt, C, T, # PDE parameters\n",
" umin, umax, # Interval for u in plots\n",
" animate=True, # Simulation with animation?\n",
" tool='matplotlib', # 'matplotlib' or 'scitools'\n",
" solver_function=solver, # Function with numerical algorithm\n",
" ):\n",
" \"\"\"Run solver, store and visualize u at each time level.\"\"\"\n",
"\n",
" class PlotUst:\n",
" def __init__(self):\n",
" self.all_u = []\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" plt.plot(x, u, 'r-',\n",
" xlabel='x', ylabel='u',\n",
" axis=[0, L, umin, umax],\n",
" title='t=%f' % t[n], show=True)\n",
" # Let the initial condition stay on the screen for 2\n",
" # seconds, else insert a pause of 0.2 s between each plot\n",
" time.sleep(2) if t[n] == 0 else time.sleep(0.2)\n",
" plt.savefig('tmp_%04d.png' % n) # for movie making \n",
" self.all_u.append(u.copy())\n",
"\n",
" class PlotMatplotlib:\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" if n == 0:\n",
" plt.ion()\n",
" self.lines = plt.plot(x, u, 'r-')\n",
" plt.xlabel('x'); plt.ylabel('u')\n",
" plt.axis([0, L, umin, umax])\n",
" plt.legend(['t=%f' % t[n]], loc='lower left')\n",
" else:\n",
" self.lines[0].set_ydata(u)\n",
" plt.legend(['t=%f' % t[n]], loc='lower left')\n",
" plt.draw()\n",
" time.sleep(2) if t[n] == 0 else time.sleep(0.2)\n",
" plt.savefig('tmp_%04d.png' % n) # for movie making\n",
"\n",
" if tool == 'matplotlib':\n",
" import matplotlib.pyplot as plt\n",
" plot_u = PlotMatplotlib()\n",
" elif tool == 'scitools':\n",
" import scitools.std as plt # scitools.easyviz interface\n",
" plot_u = PlotUst()\n",
" import time, glob, os\n",
"\n",
" # Clean up old movie frames\n",
" for filename in glob.glob('tmp_*.png'):\n",
" os.remove(filename)\n",
"\n",
" # Call solver and do the simulaton\n",
" user_action = plot_u if animate else None\n",
" u, x, t, cpu = solver_function(\n",
" I, V, f, c, L, dt, C, T, user_action)\n",
"\n",
" # Make video files\n",
" fps = 4 # frames per second\n",
" codec2ext = dict(flv='flv', libx264='mp4', libvpx='webm',\n",
" libtheora='ogg') # video formats\n",
" filespec = 'tmp_%04d.png'\n",
" movie_program = 'ffmpeg'\n",
" for codec in codec2ext:\n",
" ext = codec2ext[codec]\n",
" cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\\\n",
" '-vcodec %(codec)s movie.%(ext)s' % vars()\n",
" os.system(cmd)\n",
"\n",
" if tool == 'scitools':\n",
" # Make an HTML play for showing the animation in a browser\n",
" plt.movie('tmp_*.png', encoder='html', fps=fps,\n",
" output_file='movie.html')\n",
" return cpu, np.array(plot_u.all_u)\n",
"\n",
"def guitar(C):\n",
" \"\"\"Triangular wave (pulled guitar string).\"\"\"\n",
" L = 0.75\n",
" x0 = 0.8*L\n",
" a = 0.005\n",
" freq = 440\n",
" wavelength = 2*L\n",
" c = freq*wavelength\n",
" from math import pi\n",
" w = 2*pi*freq\n",
" num_periods = 1\n",
" T = 2*pi/w*num_periods\n",
" # Choose dt the same as the stability limit for Nx=50\n",
" dt = L/50./c\n",
"\n",
" def I(x):\n",
" return a*x/x0 if x < x0 else a/(L-x0)*(L-x)\n",
"\n",
" umin = -1.2*a; umax = -umin\n",
" cpu, all_u = viz(I, 0, 0, c, L, dt, C, T, umin, umax,\n",
" animate=True, tool='scitools')\n",
" # checking\n",
" #for e in all_u:\n",
" # print e[int(len(all_u[1])/2)]\n",
"\n",
"if __name__ == '__main__':\n",
" #test_quadratic()\n",
" import sys\n",
" try:\n",
" C = float(sys.argv[1])\n",
" print 'C=%g' % C\n",
" except IndexError:\n",
" C = 0.85\n",
" print 'Courant number: %.2f' % C\n",
" guitar(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Filename: `wave1D_u0_s2c`.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 4: Compare several Courant numbers in one movie\n",
"\n",
"\n",
"The goal of this exercise is to make movies where several curves,\n",
"corresponding to different Courant numbers, are visualized. Write a\n",
"program that resembles `wave1D_u0_s2c.py` in [Exercise 3: Use a class for the user action function](#wave:exer:store:list:class), but with a `viz` function that\n",
"can take a list of `C` values as argument and create a movie with\n",
"solutions corresponding to the given `C` values. The `plot_u` function\n",
"must be changed to store the solution in an array (see [Exercise 2: Add storage of solution in a user action function](#wave:exer:store:list) or [Exercise 3: Use a class for the user action function](#wave:exer:store:list:class) for\n",
"details), `solver` must be computed for each value of the Courant\n",
"number, and finally one must run through each time step and plot all\n",
"the spatial solution curves in one figure and store it in a file.\n",
"\n",
"The challenge in such a visualization is to ensure that the curves in\n",
"one plot correspond to the same time point. The easiest remedy is to\n",
"keep the time resolution constant and change the space resolution\n",
"to change the Courant number. Note that each spatial grid is needed for\n",
"the final plotting, so it is an option to store those grids too.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"Modifying the code to store all solutions for each $C$ value and also each corresponding spatial grid (needed for final plotting), we get"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"#!/usr/bin/env python\n",
"\"\"\"\n",
"1D wave equation with u=0 at the boundary.\n",
"Simplest possible implementation.\n",
"\n",
"The key function is::\n",
"\n",
" u, x, t, cpu = solver(I, V, f, c, L, dt, C, T, user_action)\n",
"\n",
"which solves the wave equation u_tt = c**2*u_xx on (0,L) with u=0\n",
"on x=0,L, for t in (0,T]. Initial conditions: u=I(x), u_t=V(x).\n",
"\n",
"T is the stop time for the simulation.\n",
"dt is the desired time step.\n",
"C is the Courant number (=c*dt/dx), which specifies dx.\n",
"f(x,t) is a function for the source term (can be 0 or None).\n",
"I and V are functions of x.\n",
"\n",
"user_action is a function of (u, x, t, n) where the calling\n",
"code can add visualization, error computations, etc.\n",
"\"\"\"\n",
"\n",
"import numpy as np\n",
"\n",
"def solver(I, V, f, c, L, dt, C, T, user_action=None):\n",
" \"\"\"Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\"\"\"\n",
" Nt = int(round(T/dt))\n",
" t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
" dx = c*dt/C\n",
" Nx = int(round(L/dx))\n",
" x = np.linspace(0, L, Nx+1) # Mesh points in space\n",
" C2 = C**2 # Help variable in the scheme\n",
" # Recompute to make sure dx and dt are compatible with x and t\n",
" dx = x[1] - x[0]\n",
" dt = t[1] - t[0]\n",
"\n",
" if f is None or f == 0 :\n",
" f = lambda x, t: 0\n",
" if V is None or V == 0:\n",
" V = lambda x: 0\n",
"\n",
" u = np.zeros(Nx+1) # Solution array at new time level\n",
" u_1 = np.zeros(Nx+1) # Solution at 1 time level back\n",
" u_2 = np.zeros(Nx+1) # Solution at 2 time levels back\n",
"\n",
" import time; t0 = time.clock() # for measuring CPU time\n",
"\n",
" # Load initial condition into u_1\n",
" for i in range(0,Nx+1):\n",
" u_1[i] = I(x[i])\n",
"\n",
" if user_action is not None:\n",
" user_action(u_1, x, t, 0)\n",
"\n",
" # Special formula for first time step\n",
" n = 0\n",
" for i in range(1, Nx):\n",
" u[i] = u_1[i] + dt*V(x[i]) + \\\n",
" 0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" 0.5*dt**2*f(x[i], t[n])\n",
" u[0] = 0; u[Nx] = 0\n",
"\n",
" if user_action is not None:\n",
" user_action(u, x, t, 1)\n",
"\n",
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
" for n in range(1, Nt):\n",
" # Update all inner points at time t[n+1]\n",
" for i in range(1, Nx):\n",
" u[i] = - u_2[i] + 2*u_1[i] + \\\n",
" C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n",
" dt**2*f(x[i], t[n])\n",
"\n",
" # Insert boundary conditions\n",
" u[0] = 0; u[Nx] = 0\n",
" if user_action is not None:\n",
" if user_action(u, x, t, n+1):\n",
" break\n",
"\n",
" # Switch variables before next step\n",
" u_2[:] = u_1; u_1[:] = u\n",
"\n",
" cpu_time = time.clock() - t0\n",
" return u, x, t, cpu_time\n",
"\n",
"def viz(\n",
" I, V, f, c, L, dt, C, T, # PDE parameters\n",
" umin, umax, # Interval for u in plots\n",
" animate=True, # Simulation with animation?\n",
" tool='matplotlib', # 'matplotlib' or 'scitools'\n",
" solver_function=solver, # Function with numerical algorithm\n",
" ):\n",
" \"\"\"\n",
" Run solver, store and viz. u at each time level with all C values.\n",
" \"\"\"\n",
" \n",
" class PlotUst:\n",
" def __init__(self):\n",
" self.all_u = []\n",
" self.all_u_for_all_C = []\n",
" self.x_mesh = [] # need each mesh for final plots\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" self.all_u.append(u.copy())\n",
" if t[n] == T: # i.e., whole time interv. done for this C\n",
" self.x_mesh.append(x)\n",
" self.all_u_for_all_C.append(self.all_u)\n",
" self.all_u = [] # reset to empty list\n",
" \n",
" if len(self.all_u_for_all_C) == len(C): # all C done\n",
" print 'Finished all C. Proceed with plots...'\n",
" # note: n will here be the last index in t[n]\n",
" for n_ in range(0, n+1): # for each tn\n",
" plt.plot(self.x_mesh[0], \n",
" self.all_u_for_all_C[0][n_],\n",
" axis=[0, L, umin, umax],\n",
" title='Solutions for all \\\n",
" C at t=%f' % t[n_])\n",
" plt.hold('on')\n",
" \n",
" for j in range(1, len(C)):\n",
" # build plot at this tn with each \n",
" # sol. from the different C values\n",
" plt.plot(self.x_mesh[j], \n",
" self.all_u_for_all_C[j][n_],\n",
" axis=[0, L, umin, umax])\n",
" plt.xlabel('x'); plt.ylabel('u')\n",
" plt.hold('off')\n",
" plt.show()\n",
" # Let the init. cond. stay on the screen for\n",
" # 2 sec, else insert a pause of 0.2 s \n",
" # between each plot \n",
" time.sleep(2) if t[n_] == 0 else \\\n",
" time.sleep(0.2) \n",
" plt.savefig('tmp_%04d.png' % n_) # for movie\n",
" \n",
"\n",
" class PlotMatplotlib:\n",
" def __init__(self):\n",
" self.all_u = []\n",
" self.all_u_for_all_C = []\n",
" def __call__(self, u, x, t, n):\n",
" \"\"\"user_action function for solver.\"\"\"\n",
" self.all_u.append(u.copy())\n",
" if t[n] == T: # i.e., whole time interv. done for this C\n",
" self.all_u_for_all_C.append(self.all_u)\n",
" self.all_u = [] # reset to empty list\n",
" \n",
" if len(self.all_u_for_all_C) == len(C): # all C done\n",
" print 'Finished all C. Proceed with plots...'\n",
" plt.ion()\n",
" # note: n will here be the last index in t[n]\n",
" for n_ in range(0, n+1): # for each tn\n",
" plt.plot(x, self.all_u_for_all_C[0][n_])\n",
" plt.axis([0, L, umin, umax])\n",
" plt.hold(True)\n",
" for j in range(1, len(C)):\n",
" # build plot at this tn with each \n",
" # sol. from the different C values\n",
" plt.plot(x, self.all_u_for_all_C[j][n_])\n",
" plt.axis([0, L, umin, umax])\n",
" plt.xlabel('x'); plt.ylabel('u')\n",
" plt.title('Solutions for all \\\n",
" C at t=%f' % t[n_])\n",
" plt.hold(False)\n",
" plt.draw()\n",
" # Let the init. cond. stay on the screen for\n",
" # 2 sec, else insert a pause of 0.2 s \n",
" # between each plot \n",
" time.sleep(2) if t[n_] == 0 else \\\n",
" time.sleep(0.2) \n",
" plt.savefig('tmp_%04d.png' % n_) # for movie\n",
"\n",
" if tool == 'matplotlib':\n",
" import matplotlib.pyplot as plt\n",
" plot_u = PlotMatplotlib()\n",
" elif tool == 'scitools':\n",
" import scitools.std as plt # scitools.easyviz interface\n",
" plot_u = PlotUst()\n",
" import time, glob, os\n",
"\n",
" # Clean up old movie frames\n",
" for filename in glob.glob('tmp_*.png'):\n",
" os.remove(filename)\n",
"\n",
" # Call solver and do the simulaton\n",
" user_action = plot_u if animate else None\n",
" for C_value in C:\n",
" print 'C_value --------------------------------- ', C_value\n",
" u, x, t, cpu = solver_function(\n",
" I, V, f, c, L, dt, C_value, T, user_action)\n",
"\n",
" # Make video files\n",
" fps = 4 # frames per second\n",
" codec2ext = dict(flv='flv', libx264='mp4', libvpx='webm',\n",
" libtheora='ogg') # video formats\n",
" filespec = 'tmp_%04d.png'\n",
" movie_program = 'ffmpeg'\n",
" for codec in codec2ext:\n",
" ext = codec2ext[codec]\n",
" cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\\\n",
" '-vcodec %(codec)s movie.%(ext)s' % vars()\n",
" os.system(cmd)\n",
"\n",
" if tool == 'scitools':\n",
" # Make an HTML play for showing the animation in a browser\n",
" plt.movie('tmp_*.png', encoder='html', fps=fps,\n",
" output_file='movie.html')\n",
" return cpu\n",
"\n",
"def guitar(C):\n",
" \"\"\"Triangular wave (pulled guitar string).\"\"\"\n",
" L = 0.75\n",
" x0 = 0.8*L\n",
" a = 0.005\n",
" freq = 440\n",
" wavelength = 2*L\n",
" c = freq*wavelength\n",
" from math import pi\n",
" w = 2*pi*freq\n",
" num_periods = 1\n",
" T = 2*pi/w*num_periods\n",
" # Choose dt the same as the stability limit for Nx=50\n",
" dt = L/50./c\n",
" dx = dt*c/float(C)\n",
" # Now dt is considered fixed and a list of C \n",
" # values is made by reducing increasing the dx value \n",
" # in steps of 10%. \n",
" all_C = [C]\n",
" all_C.append(c*dt/(1.1*dx))\n",
" all_C.append(c*dt/(1.2*dx))\n",
"\n",
" def I(x):\n",
" return a*x/x0 if x < x0 else a/(L-x0)*(L-x)\n",
"\n",
" umin = -1.2*a; umax = -umin\n",
" cpu = viz(I, 0, 0, c, L, dt, all_C, T, umin, umax,\n",
" animate=True, tool='scitools')\n",
" #cpu = viz(I, 0, 0, c, L, dt, all_C, T, umin, umax,\n",
" # animate=True, tool='matplotlib')\n",
" print 'cpu = ', cpu\n",
"\n",
"if __name__ == '__main__':\n",
" import sys\n",
" try:\n",
" C = float(sys.argv[1])\n",
" print 'C=%g' % C\n",
" except IndexError:\n",
" C = 0.85\n",
" print 'Courant number: %.2f' % C\n",
" # The list of C values will be generated from this C value\n",
" guitar(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Filename: `wave_numerics_comparison`.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Exercise 5: Implementing the solver function as a generator\n",
"\n",
"\n",
"The callback function `user_action(u, x, t, n)` is called from the\n",
"`solver` function (in, e.g., `wave1D_u0.py`) at every time level and lets\n",
"the user work perform desired actions with the solution, like plotting it\n",
"on the screen. We have implemented the callback function in the typical\n",
"way it would have been done in C and Fortran. Specifically, the code looks\n",
"like"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"if user_action is not None:\n",
" if user_action(u, x, t, n):\n",
" break"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Many Python programmers, however, may claim that `solver` is an iterative\n",
"process, and that iterative processes with callbacks to the user code is\n",
"more elegantly implemented as *generators*. The rest of the text has little\n",
"meaning unless you are familiar with Python generators and the `yield`\n",
"statement.\n",
"\n",
"Instead of calling `user_action`, the `solver` function\n",
"issues a `yield` statement, which is a kind of `return` statement:"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"yield u, x, t, n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The program control is directed back to the calling code:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"for u, x, t, n in solver(...):\n",
" # Do something with u at t[n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When the block is done, `solver` continues with the statement after `yield`.\n",
"Note that the functionality of terminating the solution process if\n",
"`user_action` returns a `True` value is not possible to implement in the\n",
"generator case.\n",
"\n",
"Implement the `solver` function as a generator, and plot the solution\n",
"at each time step.\n",
"Filename: `wave1D_u0_generator`.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Project 6: Calculus with 1D mesh functions\n",
"\n",
"\n",
"This project explores integration and differentiation of\n",
"mesh functions, both with scalar and vectorized implementations.\n",
"We are given a mesh function $f_i$ on a spatial one-dimensional\n",
"mesh $x_i=i\\Delta x$, $i=0,\\ldots,N_x$, over the interval $[a,b]$.\n",
"\n",
"\n",
"**a)**\n",
"Define the discrete derivative of $f_i$ by using centered\n",
"differences at internal mesh points and one-sided differences\n",
"at the end points. Implement a scalar version of\n",
"the computation in a Python function and write an associated unit test\n",
"for the linear case $f(x)=4x-2.5$ where the discrete derivative should\n",
"be exact.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"See code below.\n",
"\n",
"\n",
"\n",
"**b)**\n",
"Vectorize the implementation of the discrete derivative.\n",
"Extend the unit test to check the validity of the implementation.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"See code below.\n",
"\n",
"\n",
"\n",
"**c)**\n",
"To compute the discrete integral $F_i$ of $f_i$, we assume that\n",
"the mesh function $f_i$ varies linearly between the mesh points.\n",
"Let $f(x)$ be such a linear interpolant of $f_i$. We then\n",
"have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F_i = \\int_{x_0}^{x_i} f(x) dx\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The exact integral of a piecewise linear function $f(x)$ is\n",
"given by the Trapezoidal rule. Show\n",
"that if $F_{i}$ is already computed, we can find $F_{i+1}$\n",
"from"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F_{i+1} = F_i + \\frac{1}{2}(f_i + f_{i+1})\\Delta x\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make a function for the scalar implementation of the discrete integral\n",
"as a mesh function. That is, the function should return\n",
"$F_i$ for $i=0,\\ldots,N_x$.\n",
"For a unit test one can use the fact that the above defined\n",
"discrete integral of a linear\n",
"function (say $f(x)=4x-2.5$) is exact.\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"We know that the difference $F_{i+1} - F_i$ must amount to the area\n",
"of a trapezoid, which is exactly what $\\frac{1}{2}(f_i + f_{i+1})\\Delta x$ is.\n",
"To show the relation above, we may start with the Trapezoidal rule:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F_{i+1} = \\Delta x \\left[\\frac{1}{2}f(x_0) + \\sum_{j=1}^{n-1}f(x_j) + \\frac{1}{2}f(x_n) \\right] \\thinspace . \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since $n = i+1$, and since the final term in the sum may be separated out from the sum and split in two, this may be written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F_{i+1} = \\Delta x \\left[\\frac{1}{2}f(x_0) + \\sum_{j=1}^{i-1}f(x_j) + \\frac{1}{2}f(x_i) + \\frac{1}{2}f(x_i) + \\frac{1}{2}f(x_{i+1}) \\right] \\thinspace . \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This may further be written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F_{i+1} = \\Delta x \\left[\\frac{1}{2}f(x_0) + \\sum_{j=1}^{i-1}f(x_j) + \\frac{1}{2}f(x_i)\\right] + \\Delta x \\left[\\frac{1}{2}f(x_i) + \\frac{1}{2}f(x_{i+1}) \\right] \\thinspace . \\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, this gives"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F_{i+1} = F_i + \\frac{1}{2}(f_i + f_{i+1})\\Delta x\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See code below for implementation.\n",
"\n",
"\n",
"\n",
"**d)**\n",
"Vectorize the implementation of the discrete integral.\n",
"Extend the unit test to check the validity of the implementation.\n",
"\n",
"\n",
"\n",
"**Hint.**\n",
"Interpret the recursive formula for $F_{i+1}$ as a sum.\n",
"Make an array with each element of the sum and use the \"cumsum\"\n",
"(`numpy.cumsum`) operation to compute the accumulative sum:\n",
"`numpy.cumsum([1,3,5])` is `[1,4,9]`.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"See code below.\n",
"\n",
"\n",
"\n",
"**e)**\n",
"Create a class `MeshCalculus` that can integrate and differentiate\n",
"mesh functions. The class can just define some methods that call\n",
"the previously implemented Python functions. Here is an example\n",
"on the usage:"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"calc = MeshCalculus(vectorized=True)\n",
"x = np.linspace(0, 1, 11) # mesh\n",
"f = np.exp(x) # mesh function\n",
"df = calc.differentiate(f, x) # discrete derivative\n",
"F = calc.integrate(f, x) # discrete anti-derivative"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"**Solution.**\n",
"See code below.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"**Solution.**\n",
"The final version of the code reads"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"# -*- coding: utf-8 -*-\n",
"\"\"\"\n",
"Calculus with a 1D mesh function.\n",
"\"\"\"\n",
"import numpy as np\n",
"\n",
"class MeshCalculus:\n",
" def __init__(self, vectorized=True):\n",
" self.vectorized = vectorized\n",
"\n",
" def differentiate(self, f, x):\n",
" '''\n",
" Computes the derivative of f by centered differences, but \n",
" forw and back difference at the start and end, respectively.\n",
" '''\n",
" dx = x[1] - x[0]\n",
" Nx = len(x) - 1 # number of spatial steps\n",
" num_dfdx = np.zeros(Nx+1)\n",
" # Compute approximate derivatives at end-points first\n",
" num_dfdx[0] = (f(x[1]) - f(x[0]))/dx # FD approx.\n",
" num_dfdx[Nx] = (f(x[Nx]) - f(x[Nx-1]))/dx # BD approx.\n",
" # proceed with approximate derivatives for inner mesh points\n",
" if self.vectorized:\n",
" num_dfdx[1:-1] = (f(x[2:]) - f(x[:-2]))/(2*dx)\n",
" else: # scalar version\n",
" for i in range(1, Nx):\n",
" num_dfdx[i] = (f(x[i+1]) - f(x[i-1]))/(2*dx)\n",
" return num_dfdx\n",
" \n",
" def integrate(self, f, x):\n",
" '''\n",
" Computes the integral of f(x) over the interval \n",
" covered by x. \n",
" '''\n",
" dx = x[1] - x[0]\n",
" F = np.zeros(len(x)) \n",
" F[0] = 0 # starting value for iterative scheme\n",
" if self.vectorized:\n",
" all_trapezoids = np.zeros(len(x)-1) \n",
" all_trapezoids[:] = 0.5*(f(x[:-1]) + f(x[1:]))*dx \n",
" F[1:] = np.cumsum(all_trapezoids)\n",
" else: # scalar version\n",
" for i in range(0, len(x)-1):\n",
" F[i+1] = F[i] + 0.5*(f(x[i]) + f(x[i+1]))*dx \n",
" return F\n",
" \n",
"def test_differentiate():\n",
" def f(x):\n",
" return 4*x - 2.5\n",
" def dfdx(x):\n",
" derivatives = np.zeros(len(x))\n",
" derivatives[:] = 4\n",
" return derivatives\n",
" \n",
" a = 0; b = 1; Nx = 10\n",
" x = np.linspace(a, b, Nx+1) \n",
" exact_dfdx = dfdx(x) # compute exact derivatives\n",
" # test vectorized version\n",
" calc_v = MeshCalculus(vectorized=True)\n",
" num_dfdx = calc_v.differentiate(f, x)\n",
" print np.abs(num_dfdx - exact_dfdx)\n",
" diff = np.abs(num_dfdx - exact_dfdx).max()\n",
" tol = 1E-14\n",
" assert diff < tol\n",
" # test scalar version\n",
" calc = MeshCalculus(vectorized=False)\n",
" num_dfdx = calc.differentiate(f, x)\n",
" print np.abs(num_dfdx - exact_dfdx)\n",
" diff = np.abs(num_dfdx - exact_dfdx).max()\n",
" assert diff < tol\n",
" \n",
"def test_integrate():\n",
" def f(x):\n",
" return 4*x - 2.5\n",
" a = 0; b = 1; Nx = 10\n",
"# a = 2.5/4; b = 10; Nx = 2\n",
" x = np.linspace(a, b, Nx+1) \n",
" # The exact integral amounts to the total area of two triangles\n",
" I_exact = 0.5*abs(2.5/4 - a)*f(a) + 0.5*abs(b - 2.5/4)*f(b)\n",
" # test vectorized version \n",
" calc_v = MeshCalculus(vectorized=True) \n",
" F = calc_v.integrate(f, x)\n",
" print F, I_exact\n",
" diff = np.abs(F[-1] - I_exact)\n",
" print diff\n",
" tol = 1E-14\n",
" assert diff < tol \n",
" # test scalar version\n",
" calc = MeshCalculus(vectorized=False) \n",
" F = calc.integrate(f, x)\n",
" print F, I_exact\n",
" diff = np.abs(F[-1] - I_exact)\n",
" print diff\n",
" assert diff < tol\n",
" \n",
" \n",
"if __name__ == '__main__':\n",
" test_differentiate()\n",
" test_integrate()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"Filename: `mesh_calculus_1D`.\n",
"\n",
""
]
}
],
"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.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}