diff --git a/Week 2/5 Discrete and Fast Fourier Transforms.ipynb b/Week 2/5 Discrete and Fast Fourier Transforms.ipynb new file mode 100644 index 0000000..289d694 --- /dev/null +++ b/Week 2/5 Discrete and Fast Fourier Transforms.ipynb @@ -0,0 +1,441 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "4ec40081b048ce2f34f3f4fedbb0be10", + "grade": false, + "grade_id": "cell-98f724ece1aacb67", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "# CDS: Numerical Methods Assignments\n", + "\n", + "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n", + "\n", + "- Solutions must be submitted via the Jupyter Hub.\n", + "\n", + "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n", + "\n", + "## Submission\n", + "\n", + "1. Name all team members in the the cell below\n", + "2. make sure everything runs as expected\n", + "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n", + "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n", + "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n", + "6. submit your solutions **in time (before the deadline)**" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "team_members = \"\"" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "2d40509e944f45e5c23cc3e732547aad", + "grade": false, + "grade_id": "cell-2a0ab25436430f05", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "## Discrete and Fast Fourier Transforms (DFT and FFT)\n", + "\n", + "In the following we will implement a DFT algorithm and, based on that, a FFT algorithm. Our aim is to experience the drastic improvement of computational time in the FFT case." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "90044d2a3233721361765db06afba03b", + "grade": true, + "grade_id": "cell-abee6fbaf30772f2", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Import packages here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "bbe05b6dd0bc5b479cf66199114d7e4d", + "grade": false, + "grade_id": "cell-a1c3327dc1cad0db", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 1\n", + "\n", + "Implement a Python function $\\text{DFT(yk)}$ which returns the Fourier transform defined by\n", + "\n", + "\\begin{equation}\n", + "\\beta_j = \\sum^{N-1}_{k=0} f(x_k) e^{-ij x_k}\n", + "\\end{equation}\n", + "\n", + "with $x_k = \\frac{2\\pi k}{N}$ and $j = 0, 1, ..., N-1$. The $\\text{yk}$ should represent the array corresponding to $y_k = f(x_k)$. Please note that this definition is slightly different to the one we introduced in the lecture. Here we follow the notation of Numpy and Scipy.\n", + "\n", + "Hint: try to write the sum as a matrix-vector product and use $\\text{numpy.dot()}$ to evaluate it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "2fded00d23ce12e99ba32c09a6370b3c", + "grade": true, + "grade_id": "cell-5f6638846212c9d1", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def DFT(yk):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "15164465870c44b9b4e6328f56eaed20", + "grade": false, + "grade_id": "cell-74e9ce917ff9d690", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 2 \n", + "\n", + "Make sure your function $\\text{DFT(yk)}$ and Numpy’s FFT function $\\text{numpy.fft.fft(yk)}$ return\n", + "the same data by plotting $|\\beta_j|$ vs. $j$ for\n", + "\n", + "\\begin{equation}\n", + " y_k = f(x_k) = e^{20i x_k} + e^{40 i x_k}\n", + "\\end{equation}\n", + "and\n", + "\\begin{equation}\n", + " y_k = f(x_k) = e^{i 5 x_k^2}\n", + "\\end{equation}\n", + "\n", + "using $N = 128$ for both routines." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "14cef25098dec15c058916f4fc58c133", + "grade": true, + "grade_id": "cell-7cc28776346ee714", + "locked": false, + "points": 1, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "1674c10abab84f4d7c15bc7e5fea53b2", + "grade": false, + "grade_id": "cell-e04e5bc7ed412f64", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 3\n", + "\n", + "Analyze the evaluation-time scaling of your $\\text{DFT(yk)}$ function with the help of the timeit\n", + "module. Base your code on the following example:\n", + "\n", + "```python\n", + "import timeit\n", + "\n", + "tOut = timeit.repeat(stmt=lambda: DFT(yk), number=10, repeat=5)\n", + "tMean = np.mean(tOut)\n", + "```\n", + "This example evaluates $\\text{DFT(yk)}$ 5 × 10 times and stores the resulting 5 evaluation times in tOut. Afterwards we calculate the mean value of these 5 repetitions. \n", + "Use this example to calculate and plot the evaluation time of your $\\text{DFT(yk)}$ function for $N = 2^2, 2^3, ..., 2^M$. Depending on your implementation you might be able to go up to $M = 10$. Be careful and increase M just step by step!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "0e92c1c2d548a1c9b1476dd295415c2e", + "grade": true, + "grade_id": "cell-0ab81532ab86e322", + "locked": false, + "points": 4, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "3063d5b408e133b1930c56fa45fed54e", + "grade": false, + "grade_id": "cell-73ca4de972164356", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 4\n", + "\n", + "A very simple FFT algorithm can be derived by the following separation of the sum from\n", + "above:\n", + "\n", + "\\begin{align}\n", + " \\beta_j = \\sum^{N-1}_{k=0} f(x_k) e^{-ij \\frac{2\\pi k}{N}} \n", + " &= \\sum^{N/2 - 1}_{k=0} f(x_{2k}) e^{-ij \\frac{2\\pi 2k}{N}} \n", + " + \\sum^{N/2 - 1}_{k=0} f(x_{2k+1}) e^{-ij \\frac{2\\pi (2k+1)}{N}}\\\\ \n", + " &= \\sum^{N/2 - 1}_{k=0} f(x_{2k}) e^{-ij \\frac{2\\pi k}{N/2}}\n", + " + \\sum^{N/2 - 1}_{k=0} f(x_{2k+1}) e^{-ij \\frac{2\\pi k}{N/2}} e^{-ij \\frac{2\\pi}{N}}\\\\\n", + " &= \\beta^{\\text{even}}_j + \\beta^{\\text{odd}}_j e^{-ij \\frac{2\\pi}{N}}\n", + "\\end{align}\n", + "\n", + "where $\\beta^{\\text{even}}_j$ is the Fourier transform based on only even $k$ (or $x_k$) and $\\beta^{\\text{odd}}_j$ the Fourier transform based on only odd $k$. In case $N = 2^M$ this even/odd separation can be done again and again in a recursive way. \n", + "\n", + "Use the template below to implement a $\\text{FFT(yk)}$ function, making use of your $\\text{DFT(yk)}$ function from above. Make sure that you get the same results as before by comparing the results from $\\text{DFT(yk)}$\n", + "and $\\text{FFT(yk)}$ for both functions defined in task 2.\n", + "\n", + "```python\n", + "def FFT(yk):\n", + " \"\"\"Don't forget to write a docstring ...\n", + " \"\"\"\n", + " N = # ... get the length of yk\n", + " \n", + " assert # ... check if N is a power of 2. Hint: use the % (modulo) operator\n", + " \n", + " if(N <= 2):\n", + " return # ... call DFT with all yk points\n", + " \n", + " else:\n", + " betaEven = # ... call FFT but using just even yk points\n", + " betaOdd = # ... call FFT but using just odd yk points\n", + " \n", + " expTerms = np.exp(-1j * 2.0 * np.pi * np.arange(N) / N)\n", + " \n", + " # Remember : beta_j is periodic in j !\n", + " betaEvenFull = np.concatenate([betaEven, betaEven])\n", + " betaOddFull = np.concatenate([betaOdd, betaOdd])\n", + " \n", + " return betaEvenFull + expTerms * betaOddFull\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "fdad5c990a077e2bc2c58d4c4da46973", + "grade": true, + "grade_id": "cell-ce8233802d8ccb83", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def FFT(yk):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "33b415b478341202e57888621063fc35", + "grade": true, + "grade_id": "cell-a27cf0cb147b31ae", + "locked": false, + "points": 1, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Do your plotting here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "c1d2dd0c2543e228d48d84e94720a3a8", + "grade": false, + "grade_id": "cell-c61537808f7cce68", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 5\n", + "\n", + "Analyze the evaluation-time scaling of your $\\text{FFT(yk)}$ function with the help of the timeit module and compare it to the scaling of the $\\text{DFT(yk)}$ function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "bb4f1eee977aad469245f60fca596938", + "grade": true, + "grade_id": "cell-aaf90559928426bf", + "locked": false, + "points": 1, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb b/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb new file mode 100644 index 0000000..a83c9cd --- /dev/null +++ b/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb @@ -0,0 +1,313 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "4ec40081b048ce2f34f3f4fedbb0be10", + "grade": false, + "grade_id": "cell-98f724ece1aacb67", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "# CDS: Numerical Methods Assignments\n", + "\n", + "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n", + "\n", + "- Solutions must be submitted via the Jupyter Hub.\n", + "\n", + "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n", + "\n", + "## Submission\n", + "\n", + "1. Name all team members in the the cell below\n", + "2. make sure everything runs as expected\n", + "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n", + "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n", + "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n", + "6. submit your solutions **in time (before the deadline)**" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "team_members = \"\"" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "177d489a4104e3c95a4de1a4c7768c01", + "grade": false, + "grade_id": "cell-1e89a94d71771bb6", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "## Composite Numerical Integration: Trapezoid and Simpson Rules\n", + "\n", + "In the following we will implement the composite trapezoid and Simpson rules to calculate definite integrals. These rules are defined by\n", + "\n", + "\\begin{align}\n", + "\t\\int_a^b \\, f(x)\\, dx &\\approx \\frac{h}{2} \\left[ f(a) + 2 \\sum_{j=1}^{n-1} f(x_j) + f(b) \\right] \n", + " &\\text{trapezoid} \\\\\n", + " &\\approx \\frac{h}{3} \\left[ f(a) + 2 \\sum_{j=1}^{n/2-1} f(x_{2j}) + 4 \\sum_{j=1}^{n/2} f(x_{2j-1}) + f(b) \\right]\t \n", + " &\\text{Simpson}\n", + "\\end{align}\n", + " \n", + "with $a = x_0 < x_1 < \\dots < x_{n-1} < x_n = b$ and $x_k = a + kh$. Here $k = 0, \\dots, n$ and $h = (b-a) / n$ is the step size." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "a60a63e0450a3157dd421b394288f18a", + "grade": true, + "grade_id": "cell-44d29c12deac7ed7", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Import packages here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "e11bb7c15d840e7a9397f209769ebb66", + "grade": false, + "grade_id": "cell-ce9a56067e726f36", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 1\n", + "\n", + "Implement both integration schemes as Python functions $\\text{trapz(yk, dx)}$ and $\\text{simps(yk, dx)}$. The argument $\\text{yk}$ is an array of length $n+1$ representing $y_k = f(x_k)$ and $\\text{dx}$ is the step size $h$. Compare your results with Scipy's functions $\\text{scipy.integrate.trapz(yk, xk)}$ and $\\text{scipy.integrate.simps(yk, xk)}$ for a $f(x_k)$ of your choice.\n", + "\n", + "Try both even and odd $n$. What do you see? Why?\n", + "\n", + "Hint: go to the Scipy documentation. How are even and odd $n$ handled there?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "2bc6bd3985c2b7ab4ab051ebe94496f9", + "grade": true, + "grade_id": "cell-59f0de06f77dce3e", + "locked": false, + "points": 6, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def trapz(yk, dx):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "def simps(yk, dx):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "9599217f233689affb19148157e62b41", + "grade": true, + "grade_id": "cell-ff04b1d785ea895f", + "locked": false, + "points": 1, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Compare your results here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "f3a2f1f2b9ba3ffeb8646c346797d95a", + "grade": false, + "grade_id": "cell-1a7e33464e3be83b", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 2\n", + "\n", + "Implement at least one test function for each of your integration functions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "c3b7ee0d3ced97054590e89bba97e031", + "grade": true, + "grade_id": "cell-d8f2e0aa55860e08", + "locked": false, + "points": 6, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def test_trapz():\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "def test_simps():\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "test_trapz()\n", + "test_simps()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "ead1d68798b82e5c9c5dba354a255abb", + "grade": false, + "grade_id": "cell-71d20f6b94c6ed05", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 3\n", + "\n", + "Study the accuracy of these integration routines by calculating the following integrals for a variety of step sizes $h$:\n", + "\n", + "- $\\int_0^1 \\, x\\, dx$\n", + "- $\\int_0^1 \\, x^2\\, dx$\n", + "- $\\int_0^1 \\, x^\\frac{1}{2}\\, dx$\n", + "\n", + "The integration error is defined as the difference (not the absolute difference) between your numerical results and the exact results. Plot the integration error as a function of $h$ for both integration routines and all listed functions. Comment on the comparison between both integration routines. Does the sign of the error match your expectations? Why?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "90eaf3d9beb2347589518aba4e8ad3c4", + "grade": true, + "grade_id": "cell-b0bb51b7eae7769b", + "locked": false, + "points": 4, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}