{
"cells": [
{
"cell_type": "markdown",
"id": "56889fcd-f08a-4475-9928-af817748f5c8",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "6c9d4b0987e070d74134bdadde3233b7",
"grade": false,
"grade_id": "cell-5e2b58ff19f0997c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"## Numerical Analysis - Fall semester 2025\n",
"# Serie 06 - Iterative methods for linear systems (continued)"
]
},
{
"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": "6cd575064c6d1b921dc2102233fad9c6",
"grade": false,
"grade_id": "cell-44ae0a875bc8bfaa",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"id": "9b66147c-34be-4f25-a102-09d342197a12",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "29724a59515dddf953c251b80f938e2f",
"grade": false,
"grade_id": "cell-2ba7c9bcdf5b69b6",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Convergence of Richarson's methods"
]
},
{
"cell_type": "markdown",
"id": "ef554769-b473-4595-8d73-b5e72b533764",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "5186a9effa303b633001167d333d807c",
"grade": false,
"grade_id": "cell-ba3e3d0fe7e9fb50",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Consider the linear system $A\\mathbf{x}=\\mathbf{b}$ where\n",
"\\begin{equation*}\n",
"\tA= \\begin{pmatrix} 2&-1&0\\\\-1&2&-1\\\\0&-1&2 \\end{pmatrix} \\text{ and } \\mathbf{b}=\\begin{pmatrix} 1\\\\0\\\\1\\end{pmatrix}\n",
"\\end{equation*}"
]
},
{
"cell_type": "markdown",
"id": "492ef7ad-b8ab-4040-8f1e-a345b657b8fa",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "8d6b51582e11488df4ee595a975a2dd1",
"grade": false,
"grade_id": "cell-43a5eb3bdd9c7f8e",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 1 (Theoretical):** Determine if Jacobi's method is guaranteed to converge. \n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "568f339d-843b-4ed4-8e0b-657c70d0fe81",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "a10168e2f69d1d9c93f427d71afbb23a",
"grade": false,
"grade_id": "cell-47917b531846e05f",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
" Now consider the following iterative method\n",
" \\begin{equation}\n",
"\t L\\mathbf{x}^{(k+1)}=L\\mathbf{x}^{(k)}+\\delta(\\mathbf{b}-A\\mathbf{x}^{(k)}) \\quad k\\geq 0 \\qquad \\qquad (\\ast)\n",
" \\end{equation}\n",
"where $\\delta>0$ is a parameter and\n",
"\\begin{equation*}\n",
"\tL:=\\begin{pmatrix}2&0&0\\\\-1&2&0\\\\0&-1&2\\end{pmatrix}\n",
"\\end{equation*}\n"
]
},
{
"cell_type": "markdown",
"id": "44d2bd4a-d692-40a7-b59b-bc8d072537cc",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "a6ada9f1935133020f17ebc58a423367",
"grade": false,
"grade_id": "cell-458b1a1c4a66c6fb",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 2 (Theoretical):** Rewrite the method $(\\ast)$ in the form $\\mathbf{x}^{(k+1)}= B_\\delta\\mathbf{x}^{(k)}+\\mathbf{z}_\\delta, k\\geq0$, for a suitable matrix $B_\\delta$ and vector $\\mathbf{z}_{\\delta}$ which is to be determined.\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "27fdf655-0917-4bfd-a906-eccc83f19da7",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "4269f4c607ab07dd185e55b4af645eda",
"grade": false,
"grade_id": "cell-141db8e0f524fdf8",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 3 (Theoretical):** Establish for which values of the parameter $\\delta> 0 $ the method $(\\ast)$ converges.\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "cbef3a28-04e0-4e51-82a0-7bbc9e4c2fdc",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "d3743ebeb6de399f350776f8bff683c0",
"grade": false,
"grade_id": "cell-f398bd87fc88c149",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Gradient method and generalizations"
]
},
{
"cell_type": "markdown",
"id": "86399130-c4bf-4ea1-9756-89ad8e025860",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "62ecdae3e720d0ef768c46d3cb37e622",
"grade": false,
"grade_id": "cell-f4169769f3e0965a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 4:** Complete the function `gradient_method` which implements the gradient method (Algorithm 5.3 in lecture notes) for iteratively solving the linear system $A \\mathbf{x} = \\mathbf{b}$.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c63f6964-c104-4331-9672-c963d0483768",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "94ab1825a86d18ada3233bafda0e6507",
"grade": false,
"grade_id": "cell-9a493deabc3677ef",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def gradient_method(A, b, x_0, n_max, tol):\n",
" x = x_0\n",
" r = b - A @ x\n",
" k = 0\n",
" while np.linalg.norm(r) > tol * np.linalg.norm(b) and k < n_max:\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()\n",
" return x, r, k"
]
},
{
"cell_type": "markdown",
"id": "1ab06a98-aa30-4b0b-bad0-3baebf9f20d2",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "dcfbfed30c671b6ebbe597067ab814ef",
"grade": false,
"grade_id": "cell-62c391c01606cb0b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We now want to compare the gradient method with the Jacobi method from last week. Therefore, we provide you with an implementation of the Jacobi method."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "01f7d0dc-712f-4d3f-b163-8b4f15bf5b91",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "be77a1b1d2a434797bf647c382cee06e",
"grade": false,
"grade_id": "cell-74a586f28587eaa6",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"def jacobi(A, b, x_0, n_max, tol):\n",
" P = np.diag(np.diag(A))\n",
" x = x_0\n",
" r = b - A @ x\n",
" k = 0\n",
" while np.linalg.norm(r) > tol * np.linalg.norm(b) and k < n_max:\n",
" z = np.linalg.solve(P, r)\n",
" x = x + z\n",
" r = b - A @ x\n",
" k = k + 1\n",
" return x, r, k"
]
},
{
"cell_type": "markdown",
"id": "aea1ef05-c4a0-480a-a06b-1e061f361191",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "11663510c855317a81fdb3724ca7d2e7",
"grade": false,
"grade_id": "cell-d5c5e63cef0e1806",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"As a test case, we consider the symmetric positive definite Laplacian matrix $A \\in \\mathbb{R}^{m \\times m}$. You can generate it with the following function:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "deb82fa8-f451-414e-8a44-9337bab0655d",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "24ad07f29c52b365c49fa02b2fda0c05",
"grade": false,
"grade_id": "cell-dbd3e282f0e9f9f8",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"def laplacian_matrix(m):\n",
" n = int(np.sqrt(m))\n",
" d = np.ones(n ** 2)\n",
" mat = sp.sparse.spdiags([- d, 2 * d, - d], [-1, 0, 1], n, n)\n",
" I = sp.sparse.eye(n)\n",
" A = sp.sparse.kron(I, mat) + sp.sparse.kron(mat, I) \n",
" return A.toarray()\n",
"\n",
"A = laplacian_matrix(100)\n",
"plt.title(r\"$100 \\times 100$ Laplacian matrix\")\n",
"plt.spy(A)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "4449e7ab-df68-4bdb-be48-8aa68d14518a",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "b41d98fe41bac15b2d39cb479f99a4d5",
"grade": false,
"grade_id": "cell-8f79118cb2cd1479",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
"\n",
"**Exercise 5:** Solve the system $A \\mathbf{x} = \\mathbf{b}$ for the Laplacian matrix $A \\in \\mathbb{R}^{100 \\times 100}$ from Exercise 4 for number of iterations $n_{\\max} = 10, 20, 30, \\dots, 200$ with the Jacobi and the gradient method. We \"handpick\" the solution $\\mathbf{x} = (\\frac{1}{m}, \\frac{2}{m}, \\dots, \\frac{m}{m})^{\\top}$, and set $\\mathbf{b} = A \\mathbf{x}$. Use $\\mathbf{x}^{(0)} = \\boldsymbol{0}$ and take `tol = 0` such that the stopping criterion won't be satisfied. Plot the errors $\\lVert \\mathbf{x}_c - \\mathbf{x} \\rVert$ of the iterative solution $\\mathbf{x}_c$ for both methods against the number of iterations $n_{\\max}$ for a logarithmic $y$-axis (`plt.semilogy`). Explain what you observe.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49812b10-55fd-4698-807d-1c4b79ec7697",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "add69605c519f263fdb5aac35eeacb62",
"grade": false,
"grade_id": "cell-d58e8c1f9b517430",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "88e4b840-2b24-482a-8498-b2780b459fc9",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "be36d64a1ebf4903d9612fcd9f5f9c9b",
"grade": false,
"grade_id": "cell-63dcdba40019bec4",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
"\n",
"**Exercise 6:** Repeat the previous exercise to compare the gradient method with the conjugate gradient method. Explain the difference by analyzing the condition number for the $2$-norm of $A$ using `np.linalg.cond(A)`.\n",
"\n",
"*Hint:* You can use `x, _ = sp.sparse.linalg.cg(A, b, x_0, maxiter=n_max, rtol=tol)` to obtain the conjugate gradient solution `x` of the linear system with matrix `A` and right-hand side `b`, with starting vector `x_0`, maximum number of iterations `n_max`, and tolerance `tol`.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5126b589-9a14-440d-b643-b12043b7b8c5",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "ddf69aabd028591f317fc833ca53713d",
"grade": false,
"grade_id": "cell-23d8f61d4013aafd",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
]
},
{
"cell_type": "markdown",
"id": "2443aa45-0fca-4113-995c-70e9adc40d89",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "429a879875accc24ff544dcadfcf2e1c",
"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 sixth 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
}