{ "cells": [ { "cell_type": "markdown", "id": "56889fcd-f08a-4475-9928-af817748f5c8", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1c1f4137ca4dc1719d0a9183f952507f", "grade": false, "grade_id": "cell-5e2b58ff19f0997c", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "## Numerical Analysis - Fall semester 2025\n", "# Serie 09 - Splines, numerical differentiation and integration" ] }, { "cell_type": "markdown", "id": "5d0c899b-58da-49e8-ace0-828c8d5b60ca", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "10edd7f8d5f42ef44d458651f4e462e8", "grade": false, "grade_id": "cell-e56f5ef08590e711", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "First, we will need to import some of the usual packages. You will have to run this cell every time you restart your notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "7d640d3e-d688-47ad-95b6-96eed7152684", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "c421cd6444adfd593791bcaa50c264d3", "grade": false, "grade_id": "cell-6bbe893666203fcf", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy as sp" ] }, { "cell_type": "markdown", "id": "7199b23b-ac3a-4087-9c30-83d8249316d2", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "31f22b9e563d37eef98227b45a848d62", "grade": false, "grade_id": "cell-3141cab342a74bd8", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Interpolation of the Runge function" ] }, { "cell_type": "markdown", "id": "041c94fb-2730-4dfc-951c-7a333a16fcff", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1c8664b11da13180d29f135705200651", "grade": false, "grade_id": "cell-cf07ec01b52c57c1", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "As in the previous exercise sheet, we are given the data $(x_i,y_i)$, $i=0,1,\\ldots,n$, where the $x_i$ are $n+1$ nodes in the interval $[a,b]=[-5,5]$ and\n", "the values $y_i$ come from evaluation of the *Runge* function\n", "$$\n", " g(x)= \\frac{1}{1+x^2}\n", "$$\n", "without any measurement error, namely $y_i = g(x_i)$, and want to analyze the accuracy of different polynomial interpolation strategies of $g$ for the data $(x_i,y_i)$. To measure the accuracy, we compute the maximum error the interpolant makes on a set of points `x_eval` specified by the user. The below two functions do this for the interpolants with uniform and Chebyshev nodes based on last week's implementations. " ] }, { "cell_type": "code", "execution_count": null, "id": "be0932b8-4fa7-4edf-8ed7-12aee1ad56db", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "055fa810421d791ab24d0c0325950253", "grade": false, "grade_id": "cell-0ccc0945dcf855d4", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "def uniform_interpolation_error(g, a, b, n, x_eval):\n", " x_unif = a + np.arange(n + 1) / n * (b - a)\n", " y_unif = g(x_unif)\n", " a_unif = np.polyfit(x_unif, y_unif, n)\n", " y_eval = np.polyval(a_unif, x_eval)\n", " error = np.max(np.abs(y_eval - g(x_eval)))\n", " return error\n", "\n", "def chebyshev_interpolation_error(g, a, b, n, x_eval):\n", " x_cheb = (a + b) / 2 - (b - a) / 2 * np.cos(np.pi * (np.arange(n + 1) + 0.5) / (n + 1)) \n", " y_cheb = g(x_cheb)\n", " a_cheb = np.polyfit(x_cheb, y_cheb, n)\n", " y_eval = np.polyval(a_cheb, x_eval)\n", " error = np.max(np.abs(y_eval - g(x_eval)))\n", " return error" ] }, { "cell_type": "markdown", "id": "1807a8b7-bc33-4a8c-9ede-a33e99073484", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "6a409c87cce7a4b86c05dc3d0c7c5fea", "grade": false, "grade_id": "cell-5b489d6cc8e222e5", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "**Exercise 1:** Create an equivalent function for cubic spline interpolation.\n", "\n", "*Hint:* As last week, you can use the SciPy function `spline = sp.interpolate.CubicSpline(x, y)` to interpolate some data `x` and `y` with a cubic spline, and return a function `spline`, which can be called to evaluate the spline at some points. \n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "530ada3b-f8cb-472e-a212-72994545bc70", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "7381a371bc862493ad2685c48642b102", "grade": false, "grade_id": "cell-966694aa6c140d06", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def spline_interpolation_error(g, a, b, n, x_eval):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", " return error" ] }, { "cell_type": "markdown", "id": "0ad574c7-4a26-4453-968b-94a2cb0a9817", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "ecd5855b25cc0a478a90a49a3181a105", "grade": false, "grade_id": "cell-3713f2f8c72dc60b", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "**Exercise 2:** For $n=1, 3, 5, ..., 25$, compute the errors of the three methods (uniform, Chebyshev, spline) and visualize them in the same plot with a logarithmic $y$-axis. Choose `x_eval` as 1'000 uniformly spaced values in $[-5, 5]$. What do you observe?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "a4403392-5586-4756-ba1f-310ebd20a19d", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "e31fcb4d52c170b304d7d113de1ebbc8", "grade": false, "grade_id": "cell-7ef3e42056d8af80", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def g(x):\n", " return 1 / (1 + x ** 2)\n", "\n", "a = -5\n", "b = 5\n", "x_eval = np.linspace(-5, 5, 1000)\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "6483f48d-323a-42cb-8ee9-83da5c2ef338", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "7560692374465eda8393303a237e11ae", "grade": false, "grade_id": "cell-4ed297a599d137b1", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Forward and central finite differences" ] }, { "cell_type": "markdown", "id": "86016eb9-2439-4dd4-bd48-5810584540b2", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "aa42ef050b9d554192345ed828e15d18", "grade": false, "grade_id": "cell-36f727029cd3481f", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 3:** Implement the functions `delta_h_backward` and `delta_h_central` which, given a function $f$, a point $\\overline{x}$, and a value $h$, compute the backward finite differences approximation \n", "\\begin{equation}\n", "\\delta_h^{-}f(\\overline{x}) = \\frac{f(\\overline{x}) - f(\\overline{x} - h)}{h}\n", "\\end{equation}\n", "and the central finite difference approximation \n", "\\begin{equation}\n", "\\delta_h^{c}f(\\overline{x}) = \\frac{f(\\overline{x} + h) - f(\\overline{x} - h)}{2h}\n", "\\end{equation}\n", "respectively.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "36db09e5-d835-4c12-abd2-9c0da7cae88a", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "a33b6e6ceb352efc583043445defdcd9", "grade": false, "grade_id": "cell-3bae7b0b57da42c1", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def delta_h_backward(f, x_bar, h):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "def delta_h_central(f, x_bar, h):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "5839fe26-5f85-4d0c-8af4-bdb796c7ad2f", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "c67844a1f4b3de28258f5179a24a7284", "grade": false, "grade_id": "cell-9d11fe8c9f9f17c4", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "Let us run a few simple checks on your implementation, to see if there are no issues." ] }, { "cell_type": "code", "execution_count": null, "id": "d394e023-0442-494a-970e-556f5195d1d3", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "9de14ef563f29a1ecbfc981874a0b836", "grade": false, "grade_id": "cell-3f4849577575e280", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "assert np.isclose(D := delta_h_backward(lambda x: x, 0, 0.1), 1), f\"'delta_h_backward(f, 0, 0.1)' for 'f(x) = x' should return 'f'(0) ≈ 1', but got {D}\";assert np.isclose(D := delta_h_central(lambda x: x, 0, 0.1), 1), f\"'delta_h_central(f, 0, 0.1)' for 'f(x) = x' should return 'f'(0) ≈ 1', but got {D}\";assert np.isclose(D := delta_h_backward(lambda x: x**2, 0.5, 1e-10), 1), f\"'delta_h_backward(f, 0, 0.1)' for 'f(x) = x^2' should return 'f'(1/2) ≈ 1', but got {D}\";assert np.isclose(D := delta_h_central(lambda x: x**2, 0.5, 1e-10), 1), f\"'delta_h_central(f, 0, 0.1)' for 'f(x) = x^2' should return 'f'(1/2) ≈ 1', but got {D}\";print(\"Nice! Your function worked well on the simple examples.\")" ] }, { "cell_type": "markdown", "id": "0613dd97-3b02-4b78-9e67-4ad7e3fba93e", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "b2b0913a3e69b35b8e617b68fe6d1366", "grade": false, "grade_id": "cell-d3b2430fda49eada", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "We now consider the function\n", "\n", "\\begin{equation}\n", "f(x)=x\\log(x)-\\sin^2(x)\n", "\\end{equation}\n", "\n", "for $x>0$ for which we want to approximate the first derivative in\n", "$\\bar{x}=1.$" ] }, { "cell_type": "markdown", "id": "76b80441-c2fd-4aa6-b171-e644a68fba3d", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "91185cbc95e7ca4eb4ee707e059c13a7", "grade": false, "grade_id": "cell-2a8bcca3c6359c15", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "**Exercise 4:** Approximate the derivative of $f$ using the two methods you have implemented in the previous exercises for $h = 0.1 \\cdot (\\frac{1}{2})^{i}$ for $i=0,1,2,3,4,5$, and compute for each $h$ the error $\\varepsilon_h=\\vert f^{'}(\\bar{x})-\\delta_h (\\bar{x})\\vert$, where $\\delta_h$ is either the forward finite differences or the central finite differences operator. Plot the errors $\\varepsilon_h$ in terms of $h$ on a graph in logarithmic scale. Comment on the obtained results.\n", "\n", "*Hint:* Use `plt.loglog` and compare the graphs of $\\varepsilon_h$ with the graphs defined by the curves $(h,h)$ and $(h,h^2)$.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "6d9ae130-d571-42ed-a051-48fd44e41070", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "0996a46bad93790a7aae9f9cb413b62e", "grade": false, "grade_id": "cell-4025f88f67eb6a6a", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def f(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "def df(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "c15b4d35-60b1-4245-bfad-184aec6ce9bc", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "96219d9f77104cccd5cebf4574fefb32", "grade": false, "grade_id": "cell-73274d942c17aede", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "**Exercise 5:** Repeat the previous exercise for $i = 0, 1,\\ldots, 30$.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "a974f3f0-0f76-4bf2-b43e-93392e7fbecf", "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "1d77e1b4a64fea5681a312b9e66baf4b", "grade": false, "grade_id": "cell-149a97ac6376015a", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "8aa303c9-d65e-49f9-a569-bf54a6883f7a", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "95871a8915cd988212968ec054468fa0", "grade": false, "grade_id": "cell-384ca392df10ebfa", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": true } }, "source": [ "
\n", " \n", "**Exercise 6 (Theoretical):** Similarly like you have seen in the lecture, justify why in the previous exercise the error is minimized for a particular value of $h$. \n", "
\n" ] }, { "cell_type": "markdown", "id": "e06c1508-431f-417a-a736-482d984265e2", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "45de3e5d5f6c5f0b1f82e774a7931a85", "grade": false, "grade_id": "cell-9abeb0ed6a9a88e3", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Mysterious finite differences formula" ] }, { "cell_type": "markdown", "id": "a443a391-6fbb-4869-972b-8c683202e88c", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1e2fa83ee4a67f552b43ae66873f3967", "grade": false, "grade_id": "cell-743b7a3e7bc0a719", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": true } }, "source": [ "
\n", " \n", "**Exercise 7 (Theoretical):** Consider the general finite differences formula\n", "\n", "$$\n", "D_hf({x})=\\frac{\\frac{1}{3} f({x}+h)+\\beta\\, f({x})- (\\frac{1}{2}+\\beta)\\,f(x-h)+\n", "\\frac{1}{6}f(x-2h)}{h}\n", "$$\n", "\n", "to approximate $f'(x)$ with $\\beta\\, \\in \\mathbb{R}$ independent of $h$. If we suppose that the function $f$ is regular enough, which of the following statements is true? [Only one answer is possible]\n", "\n", "1. If $\\beta = \\frac{1}{2}$, the formula $D_h f$ is of second order but not of third order.\n", "\n", "2. If $\\beta \\neq \\frac{1}{2}$, the formula does not approximate $f'(x)$ when $h\\rightarrow 0$.\n", "\n", "3. If $\\beta = -\\frac{5}{6}$, the formula $D_h f$ is of first order.\n", "\n", "4. The formula $D_h f$ converges to $f'(x)$ for all $\\beta \\in \\mathbb{R}$ when $h\\rightarrow 0$.\n", "\n", "5. If $\\beta = -\\frac{1}{6}$, the formula $D_h f$ is of first order.\n", "\n", "*Hint:* Insert the Taylor series expansions of $f(x + h)$, $f(x - h)$, and $f(x - 2h)$ into the formula and group all terms of the same derivative order ($f(x),f'(x),f''(x), \\dots$) together to see how it behaves depending on $\\beta$ and $h$.\n", "
\n" ] }, { "cell_type": "markdown", "id": "d80f1e35-4574-4fd1-a540-c41347806321", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "604ae2eea0e4f4517c6f5eaab013f05e", "grade": false, "grade_id": "cell-6d7f999d0bbb6be0", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Mysterious finite differences formula" ] }, { "cell_type": "markdown", "id": "70c95a59-1acb-4249-b216-c837da0ab4a5", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "bedaaec5485f9af1d9bf3ce957b68315", "grade": false, "grade_id": "cell-aa5d4fd1d4f799a0", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": true } }, "source": [ "
\n", " \n", "**Exercise 8 (Theoretical):** Compute the approximation errors for the trapezoidal rule and Simpson rule applied to the integrals:\n", "$$\n", " \\int_0^1 x^4\\,\\mathrm{d}x \\quad\\text{and}\\quad \\int_{0}^1 x^5\\,\\mathrm{d}x\n", "$$\n", "
\n" ] }, { "cell_type": "markdown", "id": "9697b5ac-be57-41a6-b36e-08a66e112a68", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "b6162c5e76d8c6eea21b505d95935b8b", "grade": false, "grade_id": "cell-169300f28af196fc", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": true } }, "source": [ "
\n", " \n", "**Exercise 9 (Theoretical):** Find $C \\in \\mathbb{R}$ such that the trapezoidal rule gives the exact result for the integral\n", "$$\n", " \\int_0^1 x^5 - C x^4\\,\\mathrm{d}x .\n", "$$\n", "
\n" ] }, { "cell_type": "markdown", "id": "1ba42156-1c02-4a00-b8b1-fcc789dd6435", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "0bceed9ea07b0657952e95ed45e2537b", "grade": false, "grade_id": "cell-5b39b5b957674ceb", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "## The end\n", "\n", "Splendid! You have reached the end of the ninth exercise notebook." ] } ], "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }