{
"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 06 - Iterative methods for linear systems (continued)"
]
},
{
"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": "5afb622e-c357-4a77-9356-e75ad0400743",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"We will investigate the spectral radius of the iteration matrix $B_J$:\n",
" \\begin{align*}\n",
" B_J=I-D^{-1}A&=\\begin{pmatrix} 1&0&0\\\\0&1&0\\\\0&0&1\\end{pmatrix}-\\begin{pmatrix} 1/2&0&0\\\\0&1/2&0\\\\0&0&1/2\\end{pmatrix}\\begin{pmatrix} 2&-1&0\\\\-1&2&-1\\\\0&-1&2\\end{pmatrix}\\\\&=\\begin{pmatrix} 0&1/2&0\\\\1/2&0&1/2\\\\0&1/2&0 \\end{pmatrix}\n",
" \\end{align*}\n",
" which has eigenvalues $\\lambda_1=0,\\lambda_{2,3}=\\pm1/\\sqrt{2}$ Hence, $\\rho(B_J)=1/\\sqrt{2}<1$. We therefore conclude that the Jacobi method will converge. \n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "568f339d-843b-4ed4-8e0b-657c70d0fe81",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"We have\n",
"\\begin{align*}\n",
" L\\mathbf{x}^{(k+1)}&=L\\mathbf{x}^{(k)}+\\delta\\mathbf{b}-A\\mathbf{x}^{(k)}\\\\\n",
" \\Rightarrow L\\mathbf{x}^{(k+1)}&=(L-\\delta A)\\mathbf{x}^{(k)}+\\delta\\mathbf{b}\n",
"\\end{align*}\n",
"and observe that $L$ is invertible. Hence, we write\n",
"\\begin{align*}\n",
" \\mathbf{x}^{(k+1)}=L^{-1}(L-\\delta A)\\mathbf{x}^{(k)}+\\delta L^{-1}\\mathbf{b} = (I-\\delta L^{-1}A)\\mathbf{x}^{(k)}+\\delta L^{-1}\\mathbf{b}. \n",
"\\end{align*}\n",
"We note that\n",
"\\begin{equation*}\n",
" L^{-1}=\\begin{pmatrix}1/2&0&0\\\\1/4&1/2&0\\\\1/8&1/4&1/2\\end{pmatrix},\n",
"\\end{equation*}\n",
"Therefore, \n",
"\\begin{equation*}\n",
" B_{\\delta}=\\begin{pmatrix}1-\\delta&\\delta/2&0\\\\0&1-3\\delta/4&\\delta/2\\\\0&\\delta/8&1-3\\delta/4\\end{pmatrix}.\n",
"\\end{equation*}\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "27fdf655-0917-4bfd-a906-eccc83f19da7",
"metadata": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"The eigenvalues of $B_\\delta$ are the zeros of its characteristic polynomial\n",
"\\begin{equation*}\n",
" (1-\\delta-\\lambda)\\left(\\lambda^2-2\\lambda\\left(1-\\frac{3}{4}\\delta\\right)+\\left(1-\\frac{3}{4}\\delta\\right)^2-\\frac{\\delta^2}{16}\\right)=0,\n",
"\\end{equation*}\n",
"which, after some algebraic manipulations, are\n",
"\\begin{equation*}\n",
" \\lambda_1=1-\\delta,\\hspace{0.3cm}\\lambda_2=1-\\delta,\\hspace{0.3cm}\\lambda_3=1-\\delta/2.\n",
"\\end{equation*}\n",
"Hence,\n",
"\\begin{equation*}\n",
" \\rho(B_\\delta)=\\max\\left\\{\\left|1-\\delta\\right|,\\left|1-\\delta/2\\right|\\right\\},\\hspace{0.3cm}\\delta>0.\n",
"\\end{equation*}\n",
"from which we conclude that $\\rho(B_{\\delta}) < 1$ whenever $0<\\delta<2$.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "cbef3a28-04e0-4e51-82a0-7bbc9e4c2fdc",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" w = A @ r\n",
" alpha = r.T @ r / (r.T @ w)\n",
" x = x + alpha * r\n",
" r = r - alpha * w\n",
" k = k + 1\n",
" ### END SOLUTION\n",
" return x, r, k"
]
},
{
"cell_type": "markdown",
"id": "1ab06a98-aa30-4b0b-bad0-3baebf9f20d2",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "cell-d58e8c1f9b517430",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"m = 100\n",
"A = laplacian_matrix(m)\n",
"x_ex = np.arange(1, m + 1) / m\n",
"b = A @ x_ex\n",
"x_0 = np.zeros(m)\n",
"tol = 0\n",
"\n",
"err_jacobi_list = []\n",
"err_gradient_list = []\n",
"\n",
"n_max_list = 10 * np.arange(1, 21)\n",
"for n_max in n_max_list:\n",
" x_jacobi, _, _ = jacobi(A, b, x_0, n_max, tol)\n",
" x_gradient, _, _ = gradient_method(A, b, x_0, n_max, tol)\n",
" err_jacobi_list.append(np.linalg.norm(x_jacobi - x_ex))\n",
" err_gradient_list.append(np.linalg.norm(x_gradient - x_ex))\n",
"\n",
"plt.semilogy(n_max_list, err_jacobi_list, label=r\"Jacobi method\")\n",
"plt.semilogy(n_max_list, err_gradient_list, label=r\"gradient method\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"number of iterations $n_{\\max}$\")\n",
"plt.ylabel(r\"absolute error $|| \\mathbf{x}_c - \\mathbf{x} ||$\")\n",
"plt.show()\n",
"\n",
"# The absolute error decays exponentially for both methods, as is ensured by\n",
"# Theorem 5.1 (Jacobi method) and Theorem 5.4 (gradient method).\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "88e4b840-2b24-482a-8498-b2780b459fc9",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "cell-23d8f61d4013aafd",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"m = 100\n",
"A = laplacian_matrix(m)\n",
"x_ex = np.arange(1, m + 1) / m\n",
"b = A @ x_ex\n",
"x_0 = np.zeros(m)\n",
"tol = 0\n",
"\n",
"err_gradient_list = []\n",
"err_conjugate_list = []\n",
"\n",
"n_max_list = 10 * np.arange(1, 21)\n",
"for n_max in n_max_list:\n",
" x_gradient, _, _ = gradient_method(A, b, x_0, n_max, tol)\n",
" x_conjugate, _ = sp.sparse.linalg.cg(A, b, x_0, maxiter=n_max, rtol=tol)\n",
" err_gradient_list.append(np.linalg.norm(x_gradient - x_ex))\n",
" err_conjugate_list.append(np.linalg.norm(x_conjugate - x_ex))\n",
"\n",
"plt.semilogy(n_max_list, err_gradient_list, label=r\"gradient method\")\n",
"plt.semilogy(n_max_list, err_conjugate_list, label=r\"conjugate gradient method\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"number of iterations $n_{\\max}$\")\n",
"plt.ylabel(r\"absolute error $|| \\mathbf{x}_c - \\mathbf{x} ||$\")\n",
"plt.show()\n",
"\n",
"cond_A = np.linalg.cond(A)\n",
"base_gradient = (cond_A - 1) / (cond_A + 1)\n",
"base_conjugate_gradient = (np.sqrt(cond_A) - 1) / (np.sqrt(cond_A) + 1)\n",
"print(\"gradient method: (κ(A) - 1) / (κ(A) + 1) = {:.3f}\".format(base_gradient))\n",
"print(\"conjugate gradient method: (√κ(A) - 1) / (√κ(A) + 1) = {:.3f}\".format(base_conjugate_gradient))\n",
"\n",
"# We notice that the conjugate gradient method converges significantly faster \n",
"# than the gradient method as the reduction factor of the error is controlled\n",
"# by √κ(A) instead of κ(A) (Theorem 5.4 vs. Theorem 5.5). Indeed, for this\n",
"# specific A, the reduction factor for the gradient method 0.959, while it is\n",
"# 0.749 for the conjugate gradient\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "2443aa45-0fca-4113-995c-70e9adc40d89",
"metadata": {
"nbgrader": {
"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
}