{
"cells": [
{
"cell_type": "markdown",
"id": "56889fcd-f08a-4475-9928-af817748f5c8",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "ac817f655f087b2a4acf370d42e32aa0",
"grade": false,
"grade_id": "cell-5e2b58ff19f0997c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Numerical Analysis - Fall semester 2025\n",
"# Serie 05 - Iterative methods for linear systems"
]
},
{
"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": "5afb622e-c357-4a77-9356-e75ad0400743",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "8d5f0a8dbed5c44de91251d3d6155f69",
"grade": false,
"grade_id": "cell-44ae0a875bc8bfaa",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import time"
]
},
{
"cell_type": "markdown",
"id": "c719caee-b4f2-4a4f-9037-f1faa22686c1",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "78051f63306b3a20b04fb6aa93bbbbd6",
"grade": false,
"grade_id": "cell-2a0ebd51555d5ee4",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Richardson's methods"
]
},
{
"cell_type": "markdown",
"id": "8a6c0438-a42a-45ee-82c5-96db83ec945c",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "ea6c15b1d6e1d9320ad50568047187b0",
"grade": false,
"grade_id": "cell-e7bec0601a76426c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 1:** Complete the below implementation of the Richardson's method with stopping criterion (cf. Algorithm 5.2 in the lecture notes).\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0af8fd7a-fdd3-4fce-a1f8-9dbccadf182f",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "66c996eed1b0f54fd136c4277b29e44f",
"grade": false,
"grade_id": "cell-1cdab4e3dd7d16d8",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def richardson(A, b, P, x_0, n_max, tol):\n",
" x = x_0\n",
" res = b - A @ x\n",
" niter = 0\n",
" while np.linalg.norm(res) > tol * np.linalg.norm(b) and niter < n_max:\n",
" z = np.linalg.solve(P, res)\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()\n",
" return x, res, niter"
]
},
{
"cell_type": "markdown",
"id": "53f7a680-78e7-4805-8d10-19182716ed7c",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "30fe66712c78012f539ef65d45f290ee",
"grade": false,
"grade_id": "cell-c6f33372b25653d9",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We now use\n",
"$$\n",
"A=\n",
"\\begin{pmatrix}\n",
" 1 & 0.5 & 0 \\\\\n",
" 0 & 2 & 1 \\\\\n",
" -1 & 1 & 1\n",
"\\end{pmatrix}, \\quad \\mathbf{b}=\n",
"\\begin{pmatrix}\n",
" 2 \\\\\n",
" 5 \\\\\n",
" 2\n",
"\\end{pmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "d1621ddb-0757-4713-954e-c561ba99dc6d",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "e18dd178b6305be7207641cc30e2d8df",
"grade": false,
"grade_id": "cell-dfc690a668a9cfce",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 2:** Implement the Jacobi method and iteratively solve the linear system $A\\mathbf{x}=\\mathbf{b}$ using the parameters $\\mathbf{x}^{(0)}=\\mathbf{0}$, $n_{\\max}=10^4$, and $tol=10^{-5}$.\n",
"\n",
"*Hint:* The Jacobi method is a Richardson's method for a specific choice of $P$, so you can use your implementation of the function `richardson` above.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f2b367d1-d9ed-4a2a-9d27-099a63e34a36",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "3e3f7e95775dcd35b797ff4c49f94ab6",
"grade": false,
"grade_id": "cell-1ca2775cef669156",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def jacobi(A, b, x_0, n_max, tol):\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()\n",
" return x, res, niter\n",
"\n",
"A = np.array([[1, 0.5, 0], [0, 2, 1], [-1, 1, 1]])\n",
"b = np.array([2, 5, 2])\n",
"x_0 = np.zeros(3)\n",
"n_max = 1e4\n",
"tol = 1e-5\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "fd5ff6d1-72cd-485c-801d-fa0e2fc99248",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "eb506589ec0bbc8110b9dd68996e2e5a",
"grade": false,
"grade_id": "cell-6c236e6268ad8289",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 3:** Repeat Exercise 2 using the Gauss-Seidel method.\n",
"\n",
"*Hint:* The NumPy function `np.tril` can be used to extract the lower triangular part of a matrix.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c1e9c75e-4b0c-4cc8-9396-f911091a318d",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "0b9f8fbb81bfb339c9af665ea40fd545",
"grade": false,
"grade_id": "cell-d16b29ce60202a28",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def gauss_seidel(A, b, x_0, n_max, tol):\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()\n",
" return x, res, niter\n",
"\n",
"\n",
"A = np.array([[1, 0.5, 0], [0, 2, 1], [-1, 1, 1]])\n",
"b = np.array([2, 5, 2])\n",
"x_0 = np.zeros(3)\n",
"n_max = 1e4\n",
"tol = 1e-5\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "06f277b2-f664-4809-8ffc-8980512eaeac",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "638eaa1f79615402d5ca0c0ee63e6ecf",
"grade": false,
"grade_id": "cell-1611b00377c077a3",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Let us now consider the new linear system $L\\mathbf{y}=\\mathbf{f}$ for the tridiagonal matrix $$\n",
" L=\\begin{pmatrix} 2.5 & -1 & & & \\\\ -1 & 2.5 & -1 & & \\\\ & \\ddots\n",
" & \\ddots & \\ddots & \\\\ & &-1 & 2.5 & -1 \\\\ &&&-1&2.5 \\end{pmatrix} \\in \\mathbb{R}^{n \\times n}, \\quad \\mathbf{f}=\n",
"\\begin{pmatrix}\n",
" 1.5 \\\\\n",
" 0.5 \\\\\n",
" 0.5 \\\\\n",
" \\vdots \\\\\n",
" 0.5 \\\\\n",
" 1.5\n",
"\\end{pmatrix} \\in \\mathbb{R}^{n}.\n",
"$$\n",
"which you can generate with the command `L = np.diag(2.5 * np.ones(n), 0) - np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)`."
]
},
{
"cell_type": "markdown",
"id": "f040cfbe-b7b6-4480-8c6e-6b366664285e",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "6b09378d780ac7d5406df9762150d092",
"grade": false,
"grade_id": "cell-c9bda5c4eb1f25c2",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 4:** For $n=4$, solve the linear system with the Jacobi and Gauss-Seidel methods you've implemented in Exercises 2 and 3. Use the initial data $\\mathbf{y}^{(0)}=\\mathbf{0}$, a tolerance of $10^{-10}$, and a maximum number of iterations $10^8$. Note down the computing time as well as the number of iterations. What can we say about the convergence of the two methods? Which one converges faster?\n",
"*Hint:* You can measure the time (in seconds) an operation takes with\n",
"```python\n",
"start_time = time.time()\n",
"# Some operation ...\n",
"end_time = time.time()\n",
"run_time = end_time - start_time\n",
"```\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "71dec30c-7c3d-4766-a1fb-b78811a0b68f",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "c99f8be58953828e4892a1750aac4813",
"grade": false,
"grade_id": "cell-d79263a4e3edcad8",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# Assemble the matrix and right-hand side\n",
"n = 4\n",
"L = np.diag(2.5 * np.ones(n), 0) - np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)\n",
"f = 0.5 * np.ones(n)\n",
"f[[0, -1]] = 1.5\n",
"\n",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "09e9a587-4ebb-456f-8cd8-7e289ab7ce5f",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "1391aa01d77561dc9ee63ea16477f764",
"grade": false,
"grade_id": "cell-2634a67813bc0524",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 5:** Solve the linear system $L \\mathbf{y} = \\mathbf{f}$ from above for $n=4,8,12,...,200$ with initial iterate $\\mathbf{y}^{(0)}=\\mathbf{0}$, tolerance $10^{-10}$, and maximum number of iterations $10^8$. How does the number of iterations in terms of $n$ behave? Knowing that the exact solution of the system is $\\mathbf{y}=(1,...,1)^\\top$ and denoting the approximate solution as $\\mathbf{y}_c$, plot the relative error $\\|\\mathbf{y}_c-\\mathbf{y}\\|/\\|\\mathbf{y}\\|$ and the normalized residual $\\|\\mathbf{f} - L \\mathbf{y}_c\\|/\\|\\mathbf{f}\\|$ in terms of $n$ in a graph in logarithmic scale. Is the residual a good estimator of the error? How does the condition number of $L$ change in terms of $n$?\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "16bdb9c3-c365-4b0e-a52d-642aa948670d",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "a9280471021bb32bfe1faae30b1bf4e7",
"grade": false,
"grade_id": "cell-f7debbcd055e2ef7",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "a695f1b2-1090-474d-b460-412f00debbf0",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "5718f147b52cca6532566c8033b4a8ba",
"grade": false,
"grade_id": "cell-4e07eebdd69727c0",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 6:** For $n=4, 8, 12, \\dots, 48$, repeat the previous exercise on the matrix\n",
"$$\n",
" L=\\begin{pmatrix} 2 & -1 & & & \\\\ -1 & 2 & -1 & & \\\\ & \\ddots\n",
" & \\ddots & \\ddots & \\\\ & &-1 & 2 & -1 \\\\ &&&-1&2 \\end{pmatrix} \\in \\mathbb{R}^{n \\times n}\n",
"$$\n",
"and the right-hand side $\\mathbf{f} = (1, 0, \\dots, 0, 1)^{\\top}$, such that the exact solution will still be the same.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "60c6de6c-3bec-4fd9-ac37-eb147019956c",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "3f26517588ad27f8d93d5a967c3a8fa2",
"grade": false,
"grade_id": "cell-53da87053bc33e40",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "f928033e-eb2b-4bee-b367-d7a01c21b67e",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "dcf281364ebd3047a2ef8eec0df23602",
"grade": false,
"grade_id": "cell-d5b27901d0dd7ae9",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"The eigenvalues of a tridiagonal matrix of the form\n",
"\n",
"$$\n",
"L = \\begin{bmatrix} a & b & & & \\\\ c & a & b & & \\\\ & \\ddots & \\ddots & \\ddots & \\\\ & & c & a & b \\\\ &&&c&a \\end{bmatrix} \\in\\mathbb{R}^{n\\times n}\n",
"$$\n",
"are $\\lambda_i(L)=a+2\\sqrt{bc}\\cos\\left(\\frac{\\pi i}{n+1}\\right)$, $i=1,...,n$. "
]
},
{
"cell_type": "markdown",
"id": "4943a41f-3fb6-4d70-a04d-5279e663634d",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "7eabf1c3100086f560a9cfd6cc3931fb",
"grade": false,
"grade_id": "cell-3d53ad275e551be4",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 7 (Theoretical):** Give the condition number for the $2$-norm $\\kappa(L)$ of the tridiagonal matrix $L$ explicitly in the case where $b=c$ and $|b|\\leq a/2$, i.e. for positive semi-definite $L$. Can this explain the numerical results obtained in the two previous exercises?\n",
"\n",
"*Hint:* For a symmetric matrix $L$, the condition number for the $2$-norm is given by $\\kappa(L) = \\lambda_{\\max}(L) / \\lambda_{\\min}(L)$.\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "19b22028-cf39-4330-bf12-c431669f3ab8",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "16ebcb9a6d2ae5f0ede08fa35be6efbe",
"grade": false,
"grade_id": "cell-14c1e0214f42f407",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Convergence of Jacobi's method (previous exam question)"
]
},
{
"cell_type": "markdown",
"id": "73d158fc-9e90-4f85-8612-f96f29c79546",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "16ed77c35e79d87f16f8d38c5511407b",
"grade": false,
"grade_id": "cell-79c845784ea50757",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Given\n",
"$$\n",
"A =\n",
"\\begin{pmatrix}\n",
" 1 & \\beta \\\\\n",
" -2 & 1\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"and $\\mathbf{b} \\in \\mathbb{R}^{2}$ which form the linear system $A\\mathbf{x}=\\mathbf{b}$."
]
},
{
"cell_type": "markdown",
"id": "cd36ad68-2945-4c3b-9b8c-3884a2707126",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "d4d1eed19c05049723b5203590e5e4dc",
"grade": false,
"grade_id": "cell-dd03d79b0a985dd3",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 8 (Theoretical):** For which values of $\\beta$ does Jacobi's method converge for any $\\mathbf{b}$ and any initial vector $\\mathbf{x}^{(0)}$? \n",
"\n",
"1. For $\\beta \\in (-1, 0]$\n",
"\n",
"2. For $\\beta = -1/2$\n",
"\n",
"3. For $\\beta \\in [0,1)$\n",
"\n",
"4. For $\\beta \\in (-1/2,1/2)$\n",
"\n",
"5. For $\\beta \\in (-1/2,\\infty)$\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "95444730-2c43-4d9c-b2ee-b25c263e0372",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "03a7b0a46d6e1a0a7c3ba302595efa52",
"grade": false,
"grade_id": "cell-a0ffcfeef32899e8",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Convergence of Jacobi's method"
]
},
{
"cell_type": "markdown",
"id": "feb91ee6-062c-4d90-9b51-bb490bc797f7",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "734a95d036274f4efc6f71fe4700e804",
"grade": false,
"grade_id": "cell-7642445980a17a6f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We consider the linear system $A \\mathbf{x}=\\mathbf{b}$ with\n",
"\n",
"$$\n",
" A = \\begin{pmatrix}\n",
" 1 & 0 & \\alpha \\\\\n",
" \\beta & 1 & 0 \\\\\n",
" 0 & \\gamma & 1 \n",
" \\end{pmatrix}\n",
" , \\qquad \\mathbf{b}= \\begin{pmatrix} 1\\\\ 3\\\\ 2\\end{pmatrix}\n",
"$$\n",
"where $\\alpha, \\beta, \\gamma \\in \\mathbb{R}$ are given."
]
},
{
"cell_type": "markdown",
"id": "e8099b42-b810-446a-9439-6dbb53688f0b",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "db487e64b116976b785be1d5d01cf053",
"grade": false,
"grade_id": "cell-62e05fed462f74a2",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 9 (Theoretical):** Write Jacobi's method to solve the linear system $A\\bf{x}=\\bf{b}$. That is, express the new iterates $x_1^{(k+1)}$, $x_2^{(k+1)}$, and $x_3^{(k+1)}$ in terms of the previous iterates $x_1^{(k)}$, $x_2^{(k)}$, and $x_3^{(k)}$.\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "03dc7df6-8a81-4dd6-b7d3-efa674e7dafa",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "b6363e478ea0a254e81ae724035ad129",
"grade": false,
"grade_id": "cell-63470da636d4d851",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 10 (Theoretical):** Compute the vector ${\\bf{x}}^{(1)}$ obtained after the first iteration, from the initial vector ${\\bf{x}}^{(0)}=(1,1,1)^\\top$.\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "e01ab5cc-0803-4266-8ed2-8a9aba466ec8",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "fb15101ab8cec9693fb3c82058c16984",
"grade": false,
"grade_id": "cell-3d76369d26e6aa3e",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 11 (Theoretical):** Find a condition on $\\alpha$, $\\beta$, $\\gamma$ which ensures that the Jacobi method converges.\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "2443aa45-0fca-4113-995c-70e9adc40d89",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "50fc7ce8acde27cc2e6ad5e34a4f40da",
"grade": false,
"grade_id": "cell-a26fae9678ff662c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"## The end\n",
"\n",
"Congratulations! You have made it to the end of the fifth 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
}