{
"cells": [
{
"cell_type": "markdown",
"id": "56889fcd-f08a-4475-9928-af817748f5c8",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-5e2b58ff19f0997c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Numerical Analysis - Fall semester 2025\n",
"# Serie 10 - Numerical integration"
]
},
{
"cell_type": "markdown",
"id": "5d0c899b-58da-49e8-ace0-828c8d5b60ca",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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"
]
},
{
"cell_type": "markdown",
"id": "f63dc15a-a707-43e2-9722-59d5ba8eb67f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-70742cef19695e2b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Gauss--Legendre quadrature"
]
},
{
"cell_type": "markdown",
"id": "c8a56847-1bb0-4299-93a3-b272b6a049d8",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ccd629c1b2daf0be",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"To estimate the integral of a continous function $f : [-1, 1] \\to \\mathbb{R}$, consider the formula\n",
"$$\n",
"2 \\sum_{i=0}^n w_i f(x_i),\n",
"$$\n",
"where $n \\in \\mathbb{N} \\setminus \\{0\\}$, $-1 \\le x_0 < \\cdots < x_n \\le 1$, and $w_0, \\dots, w_n \\in \\mathbb{R}$ are called the weights. The Newton--Cotes formula uses equispaced abscissae in $[-1, 1]$: $x_i \\coloneq -1+2i/n$ for all $i \\in \\{0, \\dots, n\\}$. Instead, the Gauss--Legendre rule uses abscissae and weights such that all polynomial functions of degree at most $2n+1$ are correctly integrated by the formula."
]
},
{
"cell_type": "markdown",
"id": "50482f9b-e62b-4065-8f03-6b0d422ca775",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-52e76849a5800363",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 1 (Theoretical):** Write a system of equations expressing this condition. How many equations and unknowns does it have?\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"Let $p : \\mathbb{R} \\to \\mathbb{R}$ be a polynomial function of degree at most $2n+1$. Then, there exist $a_0, \\dots, a_{2n+1} \\in \\mathbb{R}$ such that, for all $x \\in \\mathbb{R}$,\n",
"$$\n",
"p(x)\n",
"= \\sum_{j=0}^{2n+1} a_j x^j\n",
"= \\sum_{j=0}^n \\left(a_{2j} x^{2j} + a_{2j+1} x^{2j+1}\\right).\n",
"$$\n",
"Thus,\n",
"$$\n",
"2 \\sum_{i=0}^n w_i p(x_i)\n",
"= 2 \\sum_{j=0}^n a_{2j} \\sum_{i=0}^n w_i x_i^{2j} + 2 \\sum_{j=0}^n a_{2j+1} \\sum_{i=0}^n w_i x_i^{2j+1}\n",
"$$\n",
"and, since, for every $j \\in \\mathbb{N}$, $\\int_{-1}^1 x^j \\mathrm{d}x$ equals $\\frac{2}{j+1}$ if $j$ is even and $0$ if $j$ is odd,\n",
"$$\n",
"\\int_{-1}^1 p\n",
"= \\sum_{j=0}^n \\left(a_{2j} \\int_{-1}^1 x^{2j} \\mathrm{d}x + a_{2j+1} \\int_{-1}^1 x^{2j+1} \\mathrm{d}x\\right)\n",
"= 2 \\sum_{j=0}^n a_{2j} \\frac{1}{2j+1}.\n",
"$$\n",
"Equating both right-hand sides yields\n",
"$$\n",
"\\sum_{j=0}^n a_{2j} \\left(\\sum_{i=0}^n w_i x_i^{2j} - \\frac{1}{2j+1}\\right) + \\sum_{j=0}^n a_{2j+1} \\sum_{i=0}^n w_i x_i^{2j+1} = 0.\n",
"$$\n",
"This equality must hold for all $a_0, \\dots, a_{2n+1} \\in \\mathbb{R}$, which yields\n",
"$$\n",
"\\left\\{\\begin{array}{l}\n",
"{\\displaystyle\\sum_{i=0}^n w_i x_i^{2j} = \\frac{1}{2j+1}}\\\\[5mm]\n",
"{\\displaystyle\\sum_{i=0}^n w_i x_i^{2j+1} = 0}\n",
"\\end{array}\\right. \\text{ for all } j \\in \\{0, \\dots, n\\}.\n",
"$$\n",
"This system has $2n+2$ equations and $2n+2$ unknowns.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "ef4bd76c-97bf-4537-bd82-b0e9965b7a34",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-949cf8e190ef453a",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 2 (Theoretical):** Compute the Gauss--Legendre abscissae and weights for all $n \\in \\{1, 2\\}$.\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"For $n=1$, $w_0 = 1/2 = w_1$ and $x_1 = 1/\\sqrt{3} = -x_0$. For $n=2$, $w_0 = 5/18 = w_2$, $w_1 = 4/9$, $x_2 = \\sqrt{3/5} = -x_0$, and $x_1 = 0$.\n",
" \n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "3d2e5711-c73b-4d37-aa31-b7b48dff49da",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-00a1fcd66ddb47cd",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 3:** Complete the below function `gauss_legendre_1` which approximates the integral $\\int_{-1}^{1} f(x)~\\mathrm{d}x$ of any function $f$ using the Gauss--Legendre abscissae and weights you have derived above for $n=1$.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "51be4861-632d-426a-805c-c8fb8fcc32bf",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-29f17e8d808720cd",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def gauss_legendre_1(f):\n",
" ### BEGIN SOLUTION\n",
" return f(- 1 / np.sqrt(3)) + f(1 / np.sqrt(3))\n",
" ### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "43f58698-967e-4b32-99af-8869fc5c82aa",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-736edc0c0248dbdb",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Let us quickly verify that your implementation works properly on a few simple integrals."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f4f96042-1639-400e-b507-70a5655e3755",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-29a72364be21b449",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert not (gauss_legendre_1(lambda x: x) is None), f\"your function 'gauss_legendre_1' returned nothing, make sure to 'return' your result\";assert isinstance(G := gauss_legendre_1(lambda x: x), float), f\"expected 'gauss_legendre_1' to return a floating-point number, but got {type(G)}\";assert np.isclose(G := gauss_legendre_1(lambda x: 1), 2), f\"'gauss_legendre_1(f)' for 'f(x) = 1' should return 'I(f) = 2', but got {G}\";assert np.isclose(G := gauss_legendre_1(lambda x: x), 0), f\"'gauss_legendre_1(f)' for 'f(x) = x' should return 'I(f) = 0', but got {G}\";assert np.isclose(G := gauss_legendre_1(lambda x: x**2), 2/3), f\"'gauss_legendre_1(f)' for 'f(x) = x^2' should return 'I(f) = 2/3', but got {G}\";assert np.isclose(G := gauss_legendre_1(lambda x: x**3), 0), f\"'gauss_legendre_1(f)' for 'f(x) = x^3' should return 'I(f) = 0', but got {G}\";assert np.isclose(G := gauss_legendre_1(lambda x: x**4), 2/9), f\"'gauss_legendre_1(f)' for 'f(x) = x^4' should return 'I(f) = 2/9', but got {G}\";print(\"Nice! Your function correctly integrated our simple examples.\")"
]
},
{
"cell_type": "markdown",
"id": "b79b22f0-17aa-4d29-a138-6ab793f66c9a",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-db2df4acbeb03035",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Below, we provide you with a function which computes the composite Gauss--Legendre formula for $n = 1$ based on your implementation of the simple Gauss--Legendre formula from the previous exercise."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d696bf61-dbb8-4acb-ab02-ca357ce08f27",
"metadata": {},
"outputs": [],
"source": [
"def gauss_legendre_1_composite(n, f):\n",
" subinterval_borders = np.linspace(-1, 1, n + 1)\n",
" Q_f = 0\n",
" for i in range(n):\n",
" a = subinterval_borders[i]\n",
" b = subinterval_borders[i + 1]\n",
" f_transformed = lambda x: f((a + b) / 2 + (b - a) / 2 * x)\n",
" Q_f += gauss_legendre_1(f_transformed) * (b - a) / 2\n",
" return Q_f"
]
},
{
"cell_type": "markdown",
"id": "ff68b4d3-0aa1-4a0a-9ebc-76b1a7e2c00e",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-826834571814288b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 4:** Consider the function $f(x)= x^5 \\operatorname{sign}(x)$ where $\\operatorname{sign}(x) = 1$ if $x \\geq 0$ and $-1$ otherwise. For an increasing number of sub-intervals $n = 10^0, 10^1, \\dots, 10^3$, approximate the integral $I(f)$ with the composite Gauss--Legendre quadrature. For each $n$, compute the error $|I(f)-Q_h(f)|$ and plot it against $h = 2/n$ on a logarithmic scale. What order of convergence do you observe?\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4cad856a-95e0-490c-93a0-692588335883",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-aba5b5e2d521f90d",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def f(x):\n",
" return x ** 5 * (1 - 2 * (x < 0))\n",
"\n",
"### BEGIN SOLUTION\n",
"n_list = 10 ** np.arange(0, 4)\n",
"h_list = 2 / n_list\n",
"\n",
"errors = []\n",
"\n",
"# The integral is easy to compute by hand by separately considering the intervals \n",
"# [-1, 0] and [0, 1]\n",
"I_f = 1 / 3\n",
"\n",
"for n in n_list:\n",
" Q_h = gauss_legendre_1_composite(n, f)\n",
" errors.append(abs(I_f - Q_h))\n",
"\n",
"plt.loglog(h_list, errors, label=\"Gauss-Legendre\")\n",
"plt.loglog(h_list, h_list ** 4, linestyle=\"--\", color=\"k\", label=r\"$(h, h^4)$\")\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"error\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# We observe a convergence with order 4. Since the function f(x) = x⁵ sign(x)\n",
"# is in C⁴ and the Gauss--Legendre quadrature has degree of exactness 3.\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "3af54755-c66b-4346-a3a4-7651a552d8bb",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0fa0cb98faf8646f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 5:** The quadrature weights and nodes to approximate $\\int_{-1}^{1} f(x)~\\mathrm{d}x$ can be obtained with the function `np.polynomial.legendre.leggauss` for any degree $d \\in \\mathbb{N}$. Use them to implement the function `gauss_legendre`, which approximates the integral $\\int_{a}^{b} f(x)~\\mathrm{d}x$ with the Gauss--Legendre quadrature of arbitrary degrees $d$ and intervals $[a, b]$.\n",
"\n",
"*Hint:* Since the quadrature nodes $\\tilde{x}_1, \\tilde{x}_2, \\dots, \\tilde{x}_d$ and weights $\\tilde{w}_1, \\tilde{w}_2, \\dots, \\tilde{w}_d$ from `np.polynomial.legendre.leggauss` are only valid for integrals on the interval $[-1, 1]$, you will have to transform these weights to approximate integrals defined on the interval $[a, b]$. You can do this by transforming the quadrature nodes with\n",
"\n",
"$$\n",
"x_i = \\frac{a + b}{2} + \\frac{b - a}{2} \\tilde{x}_i, i=1, 2, \\dots, d,\n",
"$$\n",
"\n",
"and the quadrature weights with\n",
"\n",
"$$\n",
"w_i = \\frac{b - a}{2} \\tilde{w}_i, i=1, 2, \\dots, d.\n",
"$$\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c4c61760-89f7-49bb-bf76-2eb726609d49",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-646c0d5f55e9f7d7",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def gauss_legendre(a, b, d, f):\n",
" nodes, weights = np.polynomial.legendre.leggauss(d)\n",
" ### BEGIN SOLUTION\n",
" nodes_transformed = (a + b) / 2 + (b - a) / 2 * nodes\n",
" weights_transformed = (b - a) / 2 * weights\n",
" return np.sum(weights_transformed * f(nodes_transformed))\n",
" ### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "73c280f1-90f7-46e6-bdac-ce3357d7dfe4",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-df55ee99deaf9198",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We now run some tests for you to see if your function does what it is supposed to do."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ed1b8162-9596-4e82-9ae4-e36f40940a33",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-2acb93fed6bc8136",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert not (gauss_legendre(-1, 1, 2, lambda x: x) is None), f\"your function 'gauss_legendre' returned nothing, make sure to 'return' your result\";assert isinstance(G := gauss_legendre(-1, 1, 2, lambda x: x), float), f\"expected 'gauss_legendre' to return a floating-point number, but got {type(G)}\";assert np.isclose(G := gauss_legendre(-1, 1, 2, lambda x: 1), 2), f\"'gauss_legendre(-1, 1, 2, f)' for 'f(x) = 1' should return 'I(f) = 2', but got {G}\";assert np.isclose(G := gauss_legendre(-1, 1, 2, lambda x: x), 0), f\"'gauss_legendre(-1, 1, 2, f)' for 'f(x) = x' should return 'I(f) = 0', but got {G}\";assert np.isclose(G := gauss_legendre(-1, 1, 2, lambda x: x**2), 2/3), f\"'gauss_legendre(-1, 1, 2, f)' for 'f(x) = x^2' should return 'I(f) = 2/3', but got {G}\";assert np.isclose(G := gauss_legendre(-1, 1, 2, lambda x: x**3), 0), f\"'gauss_legendre(-1, 1, 2, f)' for 'f(x) = x^3' should return 'I(f) = 0', but got {G}\";assert np.isclose(G := gauss_legendre(-1, 1, 2, lambda x: x**4), 2/9), f\"'gauss_legendre(-1, 1, 2, f)' for 'f(x) = x^4' should return 'I(f) = 2/9', but got {G}\";assert np.isclose(G := gauss_legendre(0, 2, 2, lambda x: x**2), 8 / 3), f\"'gauss_legendre(0, 2, 2, f)' for 'f(x) = x^2' should return 'I(f) = 8/3', but got {G}\";assert np.isclose(G := gauss_legendre(-2, 3, 2, lambda x: x**2), 35 / 3), f\"'gauss_legendre(-2, 3, 2, f)' for 'f(x) = x^2' should return 'I(f) = 35/3', but got {G}\";print(\"Looks good!\")"
]
},
{
"cell_type": "markdown",
"id": "e1669f37-7d4f-49c5-8cd9-ba12c6324499",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-f5b97135a9488f35",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Composite Newton-Cotes formulas"
]
},
{
"cell_type": "markdown",
"id": "f76a0d23-ab0a-4956-b8ab-0eeaa754ca14",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-47b7ebc90041be6d",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We want to compute the integral\n",
"\n",
"$$\n",
"I(f)=\\int_{a}^b f(x) \\, dx,\n",
"$$\n",
"\n",
"for different functions $f$. Below, we provide you with the function implementations of two quadrature rules (`trapezium` and `simpson`) which all take the same arguments: `a` and `b` define the interval,`n` is the number of subintervals of the composite formula and `f` is the function to integrate."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8080c1a2-a3d9-40e7-842a-d58e2daf1dea",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0da9ae129de27946",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"def trapezium(a, b, n, f):\n",
" # Composite trapezium formula\n",
" # - a,b: boundaries of the integration interval\n",
" # - n: number of sub-intervals\n",
" # - f: function to integrate\n",
" h = (b - a) / n\n",
" xi = np.linspace(a, b, n + 1) # quadrature nodes\n",
" alphai = np.hstack((h / 2, np.full(n - 1, h), h / 2)) # weights\n",
" Qh_trap = np.dot(alphai, f(xi)) # quadrature formula\n",
" return Qh_trap"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b1e11187-4a5d-4c1c-b154-90dfc853cc47",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-7ea713f03d4cf201",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"def simpson(a, b, n, f):\n",
" # Composite Simpson formula\n",
" # - a,b: boundaries of the integration interval\n",
" # - n: number of sub-intervals\n",
" # - f: function to integrate\n",
" h = (b - a) / n\n",
" xi = np.linspace(a, b, n + 1)\n",
" # sub-interval boundaries\n",
" alphai = (h / 3) * np.hstack((0.5, np.ones(n - 1), 0.5))\n",
" # weights at x_i\n",
" ci = np.linspace(a + h / 2, b - h / 2, n)\n",
" # sub-interval mid-points\n",
" betai = (2 * h / 3) * np.ones(n)\n",
" # weights at c_i\n",
" Qh_simp = np.dot(alphai, f(xi)) + np.dot(betai, f(ci))\n",
" return Qh_simp"
]
},
{
"cell_type": "markdown",
"id": "d5e6ec4d-53e9-4e7f-bb12-b8b8086fb72f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-61667314b16bc0c7",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 6:** Consider the function $f_1(x)=e^{-x}\\sin(x)$ on the interval $[a,b]=[0,2]$. For an increasing number of sub-intervals $n = 10^1, 10^2, \\dots, 10^5$, approximate the integral $I(f_1)$ with the composite trapezium and Simpson formulas (use the functions `trapezium` and `simpson`). For each $n$, compute the error $|I(f_1)-Q_h(f_1)|$ and plot it against $h = (b - a)/n$ on a logarithmic scale. What order of convergence do you observe?\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6f38744c-9977-4933-94f5-a83932f30358",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a47b7a912a8a2315",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"def g_1(x):\n",
" return np.exp(-x) * np.sin(x)\n",
"\n",
"a = 0\n",
"b = 2\n",
"n_list = 10 ** np.arange(1, 6)\n",
"h_list = (b - a) / n_list\n",
"\n",
"errors_trap = []\n",
"errors_simpson = []\n",
"\n",
"# Exact integral:\n",
"I_f_1 = 1 / 2 - (np.cos(2) + np.sin(2)) / (2 * np.exp(2))\n",
"\n",
"for n in n_list:\n",
" Q_h_trap = trapezium(a, b, n, g_1)\n",
" errors_trap.append(abs(I_f_1 - Q_h_trap))\n",
" Q_h_simpson = simpson(a, b, n, g_1)\n",
" errors_simpson.append(abs(I_f_1 - Q_h_simpson))\n",
"\n",
"plt.loglog(h_list, errors_trap, label=\"trapezium\")\n",
"plt.loglog(h_list, errors_simpson, label=\"simpson\")\n",
"plt.loglog(h_list, h_list ** 2, linestyle=\"--\", color=\"k\", label=r\"$(h, h^2)$\")\n",
"plt.loglog(h_list, h_list ** 4, linestyle=\":\", color=\"k\", label=r\"$(h, h^4)$\")\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"error\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# The error of the trapezium formula goes down with the same\n",
"# slope as the graph (h, h²), i.e. every time the length of the\n",
"# sub-intervals is divided by 10, the error is reduced by a factor 100.\n",
"# These two methods are therefore of order 2.\n",
"# We further notice that the error with the trapezium formula\n",
"# is bigger (precisely twice as big) than the midpoint method.\n",
"# Simpson's method converges with order 4 (until we reach the precision of the\n",
"# machine 10⁻¹⁶).\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "3b9a6aeb-924e-4006-924d-57223f097f1f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-66a491c3bb918666",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 7:** Repeat the previous exercise for the function $f_2(x) = \\sqrt{|x|^3}$ on the interval $[a, b] = [-2, 2]$. How do the orders of convergence compare to what you have found for the function $f_1$?\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "015f9784-c234-45a5-81ce-a2f3f814e5d2",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0e38c4958fbdaafc",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"def g_2(x):\n",
" return np.sqrt(np.abs(x) ** 3)\n",
"\n",
"a = -2\n",
"b = 2\n",
"n_list = 10 ** np.arange(1, 6)\n",
"h_list = (b - a) / n_list\n",
"\n",
"errors_trap = []\n",
"errors_simpson = []\n",
"\n",
"# Exact integral:\n",
"I_f_2 = 16 * np.sqrt(2) / 5 \n",
"\n",
"for n in n_list:\n",
" Q_h_trap = trapezium(a, b, n, g_2)\n",
" errors_trap.append(abs(I_f_2 - Q_h_trap))\n",
" Q_h_simpson = simpson(a, b, n, g_2)\n",
" errors_simpson.append(abs(I_f_2 - Q_h_simpson))\n",
"\n",
"plt.loglog(h_list, errors_trap, label=\"trapezium\")\n",
"plt.loglog(h_list, errors_simpson, label=\"simpson\")\n",
"plt.loglog(h_list, h_list ** 2, linestyle=\"--\", color=\"k\", label=r\"$(h, h^2)$\")\n",
"plt.loglog(h_list, h_list ** 4, linestyle=\":\", color=\"k\", label=r\"$(h, h^4)$\")\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"error\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# We notice that the trapezium methods still converge \n",
"# with order 2, whereas Simpson's method converges with an order smaller\n",
"# than 4. This is due to the fact that the function √|x|³ is not regular\n",
"# enough. It is in fact only C¹ and not C² as its second derivative approaches\n",
"# infinity in x=0.\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "1ba42156-1c02-4a00-b8b1-fcc789dd6435",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-5b39b5b957674ceb",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"## The end\n",
"\n",
"Great! You have reached the end of the tenth 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
}