{
"cells": [
{
"cell_type": "markdown",
"id": "56889fcd-f08a-4475-9928-af817748f5c8",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "6892ed0615afdbd623633895758f4b94",
"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": {
"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": "74151612c15f3ee6f7b216b837cc400a",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "9ec8884a89da712e2fe7af5c79223ab3",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "7a3b00c07622f675bdd489faad56cfa2",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "e515310eb1f2a84a2c6ba37a085e1493",
"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"
]
},
{
"cell_type": "markdown",
"id": "ef4bd76c-97bf-4537-bd82-b0e9965b7a34",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "7acb23d3a100aff97fc34e0294ffe9da",
"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"
]
},
{
"cell_type": "markdown",
"id": "3d2e5711-c73b-4d37-aa31-b7b48dff49da",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "df77dac680a57a84329f83c2baa69cc2",
"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": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "46bdef731547d6b01fe15c9167087722",
"grade": false,
"grade_id": "cell-29f17e8d808720cd",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def gauss_legendre_1(f):\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "43f58698-967e-4b32-99af-8869fc5c82aa",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "cdcd9b42ea58c7c4af9223a666c27bb1",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "f67b615fdcc276280ecc903f536132da",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "a7808beb67232966c0565b67de0a9160",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "5f79cb470433b5f9ea5d8ed6371f61be",
"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": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "e61b96071d1d553c58a9ef194e069e2f",
"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",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "3af54755-c66b-4346-a3a4-7651a552d8bb",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "355409b645344c0e16ed288657bea92c",
"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": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "635ae6e2903b846cad0a6de7a458c9ab",
"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",
" # YOUR CODE HERE\n",
" raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "73c280f1-90f7-46e6-bdac-ce3357d7dfe4",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "195d56ce8fdcca36ce2bd8c777a99707",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "6d98629932f1ab1787fb5e561e230fea",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "7a7e8fb54b75e3f6379b89a5da264b01",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "2c1a99217e8ab8744c3e68aff6326ca6",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "72c6ce2f8f1754c45e4bc28350d81c76",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "251bd95cbf46fcd886d1736864ced228",
"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": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "1e67ba4b6ce75636a121219405cb5246",
"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": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "d3ae9e6407f2cb1d092d50371f5cba4b",
"grade": false,
"grade_id": "cell-a47b7a912a8a2315",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "3b9a6aeb-924e-4006-924d-57223f097f1f",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "83365bdddd887ab1a0f76ea3a3a7344e",
"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": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "5ec908e28e2fdfa9b098d01363235839",
"grade": false,
"grade_id": "cell-0e38c4958fbdaafc",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "1ba42156-1c02-4a00-b8b1-fcc789dd6435",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "291a38aa6763607bcdd2f50d52d4f204",
"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
}