{
"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 05 - Iterative methods for linear systems"
]
},
{
"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 matplotlib.pyplot as plt\n",
"import time"
]
},
{
"cell_type": "markdown",
"id": "c719caee-b4f2-4a4f-9037-f1faa22686c1",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" x = x + z\n",
" res = b - A @ x\n",
" niter = niter + 1\n",
" ### END SOLUTION\n",
" return x, res, niter"
]
},
{
"cell_type": "markdown",
"id": "53f7a680-78e7-4805-8d10-19182716ed7c",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" P = np.diag(np.diag(A))\n",
" x, res, niter = richardson(A, b, P, x_0, n_max, tol)\n",
" ### END SOLUTION\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",
"### BEGIN SOLUTION\n",
"x, res, niter = jacobi(A, b, x_0, n_max, tol)\n",
"print(r\"obtained solution: x = {}\".format(x))\n",
"print(r\"residual of solution: res = {}\".format(res))\n",
"print(r\"number of iterations: niter = {}\".format(niter))\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "fd5ff6d1-72cd-485c-801d-fa0e2fc99248",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
" ### BEGIN SOLUTION\n",
" P = np.tril(A)\n",
" x, res, niter = richardson(A, b, P, x_0, n_max, tol)\n",
" ### END SOLUTION\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",
"### BEGIN SOLUTION\n",
"P_gauss_seidel = np.tril(A)\n",
"x, res, niter = gauss_seidel(A, b, x_0, n_max, tol)\n",
"print(r\"obtained solution: x = {}\".format(x))\n",
"print(r\"residual of solution: res = {}\".format(res))\n",
"print(r\"number of iterations: niter = {}\".format(niter))\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "06f277b2-f664-4809-8ffc-8980512eaeac",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"### BEGIN SOLUTION\n",
"y_0 = np.zeros(n)\n",
"tol = 1e-10\n",
"n_max = 1e8\n",
" \n",
"start_time = time.time()\n",
"y, res, niter = jacobi(L, f, y_0, n_max, tol)\n",
"end_time = time.time()\n",
"run_time = end_time - start_time\n",
"print(\"Jacobi:\\n\" + 7*\"-\")\n",
"print(r\"obtained solution: x = {}\".format(y))\n",
"print(r\"residual of solution: res = {}\".format(res))\n",
"print(r\"number of iterations: niter = {}\".format(niter))\n",
"print(r\"computing time: t = {}\".format(run_time))\n",
"\n",
"start_time = time.time()\n",
"y, res, niter = gauss_seidel(L, f, y_0, n_max, tol)\n",
"end_time = time.time()\n",
"run_time = end_time - start_time\n",
"print(\"\\nGauss-Seidel:\\n\" + 13*\"-\")\n",
"print(r\"obtained solution: x = {}\".format(y))\n",
"print(r\"residual of solution: res = {}\".format(res))\n",
"print(r\"number of iterations: niter = {}\".format(niter))\n",
"print(r\"computing time: t = {}\".format(run_time))\n",
"\n",
"# We notice here that both methods converge. The Jacobi method needs 53 iterations\n",
"# and approximately while the Gauss-Seidel method takes 27 iterations. The Gauss-Seidel\n",
"# method is approximately twice as fast, since we have ρ(B_GS) = ρ(B_J)² < ρ(B_J) for the\n",
"# iteration matrices B_J (Jacobi method) and B_GS (Gauss-Seidel method). See Theorem 5.1\n",
"# in the lecture notes. \n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "09e9a587-4ebb-456f-8cd8-7e289ab7ce5f",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "cell-f7debbcd055e2ef7",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"tol = 1e-10\n",
"n_max = 1e8\n",
"\n",
"err_list_jacobi = []\n",
"res_list_jacobi = []\n",
"niter_list_jacobi = []\n",
"\n",
"err_list_gaussseidel = []\n",
"res_list_gaussseidel = []\n",
"niter_list_gaussseidel = []\n",
"\n",
"cond_list = []\n",
"\n",
"n_list = 4 * np.arange(1, 51)\n",
"for n in n_list:\n",
"\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",
" y_0 = np.zeros(n)\n",
" y_ex = np.ones(n)\n",
" \n",
" y_c_jacobi, _, niter_jacobi = jacobi(L, f, y_0, n_max, tol)\n",
" y_c_gaussseidel, _, niter_gaussseidel = gauss_seidel(L, f, y_0, n_max, tol)\n",
" \n",
" cond_list.append(np.linalg.cond(L))\n",
" \n",
" res_list_jacobi.append(np.linalg.norm(f - L @ y_c_jacobi) / np.linalg.norm(f))\n",
" err_list_jacobi.append(np.linalg.norm(y_c_jacobi - y_ex) / np.linalg.norm(y_ex))\n",
" niter_list_jacobi.append(niter_jacobi)\n",
" \n",
" res_list_gaussseidel.append(np.linalg.norm(f - L @ y_c_gaussseidel) / np.linalg.norm(f))\n",
" err_list_gaussseidel.append(np.linalg.norm(y_c_gaussseidel - y_ex) / np.linalg.norm(y_ex))\n",
" niter_list_gaussseidel.append(niter_gaussseidel)\n",
"\n",
"plt.loglog(n_list, niter_list_jacobi, label=\"Jacobi\")\n",
"plt.loglog(n_list, niter_list_gaussseidel, label=\"Gauss-Seidel\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.ylabel(r\"number of iterations\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"# The number of iterations increases and then stabilizes to 103 for\n",
"# Jacobi's method and to 57 for Gauss-Seidel.\n",
"\n",
"plt.title(\"Jacobi method\")\n",
"plt.loglog(n_list, err_list_jacobi, label=r\"relative error $\\varepsilon$\")\n",
"plt.loglog(n_list, res_list_jacobi, label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"plt.title(\"Gauss-Seidel method\")\n",
"plt.loglog(n_list, err_list_gaussseidel, label=r\"relative error $\\varepsilon$\")\n",
"plt.loglog(n_list, res_list_gaussseidel, label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"plt.loglog(n_list, cond_list)\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.ylabel(r\"condition number $\\kappa(L)$\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"# We notice that for the matrix L, the residual is a good estimator of the error.\n",
"# We can explain this using the fact that the matrix L is well conditioned for\n",
"# all the values of n (see condition number graph). For example, for n=200,\n",
"# the condition number of the matrix is approximately 9. The error and the\n",
"# residual r_k at the k-th iteration are linked by the relation\n",
"# \n",
"# ||y - y_k|| / ||y|| <= κ(L) * ||r_k|| / ||f||\n",
"#\n",
"# Therefore, the residual is a good estimator of the error when the conditioning\n",
"# of the matrix is close to 1.\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "a695f1b2-1090-474d-b460-412f00debbf0",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"grade": false,
"grade_id": "cell-53da87053bc33e40",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"tol = 1e-10\n",
"n_max = 1e8\n",
"\n",
"err_list_jacobi = []\n",
"res_list_jacobi = []\n",
"niter_list_jacobi = []\n",
"\n",
"err_list_gaussseidel = []\n",
"res_list_gaussseidel = []\n",
"niter_list_gaussseidel = []\n",
"\n",
"cond_list = []\n",
"\n",
"n_list = 4 * np.arange(1, 13)\n",
"for n in n_list:\n",
"\n",
" L = np.diag(2 * np.ones(n), 0) - np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)\n",
" f = np.zeros(n)\n",
" f[[0, -1]] = 1\n",
" y_0 = np.zeros(n)\n",
" y_ex = np.ones(n)\n",
" \n",
" y_c_jacobi, _, niter_jacobi = jacobi(L, f, y_0, n_max, tol)\n",
" y_c_gaussseidel, _, niter_gaussseidel = gauss_seidel(L, f, y_0, n_max, tol)\n",
" \n",
" cond_list.append(np.linalg.cond(L))\n",
" \n",
" res_list_jacobi.append(np.linalg.norm(f - L @ y_c_jacobi) / np.linalg.norm(f))\n",
" err_list_jacobi.append(np.linalg.norm(y_c_jacobi - y_ex) / np.linalg.norm(y_ex))\n",
" niter_list_jacobi.append(niter_jacobi)\n",
" \n",
" res_list_gaussseidel.append(np.linalg.norm(f - L @ y_c_gaussseidel) / np.linalg.norm(f))\n",
" err_list_gaussseidel.append(np.linalg.norm(y_c_gaussseidel - y_ex) / np.linalg.norm(y_ex))\n",
" niter_list_gaussseidel.append(niter_gaussseidel)\n",
"\n",
"plt.loglog(n_list, niter_list_jacobi, label=\"Jacobi\")\n",
"plt.loglog(n_list, niter_list_gaussseidel, label=\"Gauss-Seidel\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.ylabel(r\"number of iterations\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"plt.title(\"Jacobi method\")\n",
"plt.loglog(n_list, err_list_jacobi, label=r\"relative error $\\varepsilon$\")\n",
"plt.loglog(n_list, res_list_jacobi, label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"plt.title(\"Gauss-Seidel method\")\n",
"plt.loglog(n_list, err_list_gaussseidel, label=r\"relative error $\\varepsilon$\")\n",
"plt.loglog(n_list, res_list_gaussseidel, label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"plt.loglog(n_list, cond_list)\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.ylabel(r\"condition number $\\kappa(A)$\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"# The necessary number of iterations increases when n increases.\n",
"# This is due to the fact that the norm of the iteration matrices B_J and B_GS\n",
"# gets closer and closer to 1 when n increases.\n",
"\n",
"# Moreover, we notice that in this case the residual is not at all a good\n",
"# estimator of the error. This is due to the fact that the conditioning of\n",
"# the matrix increases more and more when n increases, as we can see on the\n",
"# fourth plot.\n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "f928033e-eb2b-4bee-b367-d7a01c21b67e",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"The condition number for the $2$-norm of a matrix $L$, which is symmetric and positive definite, is given by\n",
"$\\kappa(L)=\\frac{\\lambda_{\\max}(L)}{\\lambda_{\\min}(L)}$. For the considered tridiagonal\n",
"matrix $L$, the eigenvalues are given by $\\lambda_i(L)=a+2|b|\\cos(\\pi i/(n+1))$ where $|\\cdot|$\n",
"is the absolute value. As the maximum and the minimum of $\\cos(\\pi i/(n+1))$ are reached in\n",
"$i=1$ and $i=n$, respectively, the maximal an minimal eigenvalues are given by\n",
"\n",
"\\begin{align*}\n",
"\\lambda_{\\max}(L)&=a+2|b|\\cos\\left(\\frac{\\pi}{n+1}\\right) \\\\\n",
"\\lambda_{\\min}(L)&=a+2|b|\\cos\\left(\\frac{\\pi n}{n+1}\\right)=a-2|b|\\cos\\left(\\frac{\\pi}{n+1}\\right),\n",
"\\end{align*}\n",
"\n",
"and the condition number for the $2$-norm of the matrix is therefore\n",
"\n",
"$$\n",
"\\kappa(L)=\\frac{a+2|b|\\cos\\left(\\frac{\\pi}{n+1}\\right)}{a-2|b|\\cos\\left(\\frac{\\pi}{n+1}\\right)}.\n",
"$$\n",
"\n",
"With $a=2.5$ and $b=-1$ we get for $n=200$ a condition number for the $2$-norm of $\\kappa(L)\\approx 8.995$ whereas for $a=2$ and $b=-1$, we have $\\kappa(L)\\approx 16373$ which is consistent with the numerical results obtained above.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "19b22028-cf39-4330-bf12-c431669f3ab8",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"Jacobi's method converges for any $\\mathbf{b}$ and $\\mathbf{x}^{(0)}$ if and only if $\\rho(B_J)<1$, where $B_J= I - \\operatorname{diag}(A)^{-1}A$. We have\n",
"$$\n",
"B_J=\n",
"\\begin{pmatrix}\n",
" 1 & 0 \\\\\n",
" 0 & 1 \\\\\n",
"\\end{pmatrix}\n",
"-\n",
"\\begin{pmatrix}\n",
" 1 & 0 \\\\\n",
" 0 & 1 \\\\\n",
"\\end{pmatrix}^{-1}\n",
"\\begin{pmatrix}\n",
" 1 & \\beta \\\\\n",
" -2 & 1 \n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
" 0 & -\\beta \\\\\n",
" 2 & 0 \n",
"\\end{pmatrix}\n",
"$$\n",
"and therefore the eigenvalues of $B_J$ are the roots of\n",
"$$\n",
"0\n",
"= \\det(\\lambda I -B_J) \n",
"= \\det \n",
"\\begin{pmatrix}\n",
" \\lambda & \\beta \\\\\n",
" -2 & \\lambda\n",
"\\end{pmatrix}\n",
"= \\lambda^2 + 2 \\beta\n",
"\\Rightarrow\n",
"\\lambda = \\pm \\sqrt{-2\\beta} %= \\pm i \\sqrt{2 \\beta}\n",
"$$\n",
"and we have $|\\lambda_{max}|<1 \\Leftrightarrow \\big\\lvert \\sqrt{-2\\beta} \\, \\big\\lvert < 1$. \n",
"Therefore, the right answer is the fourth, namely $\\beta \\in (-1/2,1/2)$ (let us remind that the eigenvalues can be complex).\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "95444730-2c43-4d9c-b2ee-b25c263e0372",
"metadata": {
"nbgrader": {
"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": {
"nbgrader": {
"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": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"$$\n",
"\\left\\{\n",
"\\begin{array}{l}\n",
"x_1^{(k+1)}=\\frac{1}{a_{11}}\\, \\left( b_1 - a_{12}\\, x_2^{(k)} - a_{13}\\, x_3^{(k)} \\right)= 1 - \\alpha\\, x_3^{(k)} \\\\[5mm]\n",
"x_2^{(k+1)}=\\frac{1}{a_{22}}\\, \\left( b_2 - a_{21}\\, x_1^{(k)} - a_{23}\\, x_3^{(k)} \\right)= 3 - \\beta\\, x_1^{(k)} \\\\[5mm]\n",
"x_3^{(k+1)}=\\frac{1}{a_{33}}\\, \\left( b_3 - a_{31}\\, x_1^{(k)} - a_{32}\\, x_2^{(k)} \\right)= 2 - \\gamma\\, x_2^{(k)}.\n",
"\\end{array}\n",
"\\right.\n",
"$$\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "03dc7df6-8a81-4dd6-b7d3-efa674e7dafa",
"metadata": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"We just have to choose in the previous formulas $k=0$ and\n",
"$x_1^{(0)}=x_2^{(0)}=x_3^{(0)}=1$. We get\n",
"\n",
"$$\n",
"\\left\\{\n",
"\\begin{array}{l}\n",
"x_1^{(1)}= 1 - \\alpha\\, x_3^{(0)} = 1 - \\alpha \\\\[5mm]\n",
"x_2^{(1)}= 3 - \\beta\\, x_1^{(0)}=3 -\\beta \\\\[5mm]\n",
"x_3^{(1)}= 2 - \\gamma\\, x_2^{(0)} =2-\\gamma.\n",
"\\end{array}\n",
"\\right.\n",
"$$\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "e01ab5cc-0803-4266-8ed2-8a9aba466ec8",
"metadata": {
"nbgrader": {
"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",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"The iteration matrix for the Jacobi method is given by\n",
"\n",
"$$\n",
" B_J \n",
" = \\Big( \\operatorname{diag}(A) \\Big)^{-1}\\Big( \\operatorname{diag}(A)-A\\Big) \n",
" =I-\\operatorname{diag}(A)^{-1}A =\n",
" \\begin{pmatrix}\n",
" 0 & 0 & -\\alpha \\\\\n",
" -\\beta & 0 & 0 \\\\\n",
" 0 & -\\gamma & 0\n",
" \\end{pmatrix}.\n",
"$$\n",
"\n",
"The eigenvalues of $B_J$ are given by\n",
"\n",
"$$\n",
" \\det \\begin{pmatrix}\n",
" \\lambda & 0 & \\alpha \\\\\n",
" \\beta & \\lambda & 0 \\\\\n",
" 0 & \\gamma & \\lambda\n",
" \\end{pmatrix} =\n",
" \\lambda^3 + \\alpha\\beta\\gamma = 0 \\quad \\Longrightarrow \\quad \n",
" \\lambda = (-\\alpha\\beta\\gamma)^{\\frac{1}{3}}.\n",
"$$\n",
"\n",
"The method converges if and only if $\\rho(B_J) < 1$, so we require $|\\alpha\\beta\\gamma| < 1$.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"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 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
}