{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical Methods\n",
"\n",
"## Lecture 7: Numerical Linear Algebra III\n",
"\n",
"### Extra exercises - solutions"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import numpy as np\n",
"import scipy.linalg as sl"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 1: diagonally dominant matrices ($*$)\n",
"\n",
"A matrix $A$ is said to be diagonally dominant if for each row $i$ the absolute value of the diagonal element is larger than the sum of the absolute values of all the other terms of the row.\n",
"\n",
"\n",
"- write this definition in a mathematical form.\n",
"\n",
"\n",
"- write a code that checks if a matrix is diagonally dominant.\n",
"\n",
"\n",
"- test it with well chosen 2x2 and 3x3 examples.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Def*: $A = (A_{i,j})_{i=1..n, j=1..n}$ is diagonally dominant if: \n",
"$$\\forall i \\in [|1,n|]\\,,\\quad \\sum_{\\substack{j=1\\\\j\\neq i}}^n |A_{i,j}| \\leq |A_{i,i}|$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A1 is diagonally dominant\n",
"A2 is diagonally dominant\n",
"A3 is diagonally dominant\n",
"!! Row 0 is not diagonally dominant\n"
]
}
],
"source": [
"def is_diag_dom(A):\n",
" \n",
" (n,m) = A.shape\n",
" if n != m: \n",
" print(\"Error: matrix must be square, not %dx%d\" % (n,m))\n",
" exit(1)\n",
" \n",
" for i in range(n):\n",
" abs_row = np.sum(np.abs(A[i,0:i])) + np.sum(np.abs(A[i,i+1:n]))\n",
" if abs_row > abs(A[i,i]):\n",
" print(\"!! Row %d is not diagonally dominant\" %i)\n",
" return False\n",
" \n",
" return True\n",
"\n",
"\n",
"A1 = np.array([[2,1],[2,3]])\n",
"if is_diag_dom(A1): \n",
" print(\"A1 is diagonally dominant\")\n",
" \n",
"A2 = np.array([[-2,1],[2,3]])\n",
"if is_diag_dom(A2): \n",
" print(\"A2 is diagonally dominant\")\n",
" \n",
"A3 = np.array([[-2,2],[2,3]])\n",
"if is_diag_dom(A3): \n",
" print(\"A3 is diagonally dominant\")\n",
" \n",
"A4 = np.array([[-2,4],[2,3]])\n",
"if is_diag_dom(A4): \n",
" print(\"A4 is diagonally dominant\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 2: singular matrices and ill-conditioning ($*$)\n",
"\n",
"For the following matrixes, compute the determinant and the condition number, and classify them as singular, ill conditioned or well conditioned:\n",
"$$ (i)\\quad A = \n",
" \\begin{pmatrix}\n",
" 1 & 2 & 3 \\\\\n",
" 2 & 3 & 4 \\\\\n",
" 3 & 4 & 5 \\\\\n",
" \\end{pmatrix}\n",
"\\quad\\quad\\quad\\quad\n",
"(ii)\\quad A = \n",
" \\begin{pmatrix}\n",
" 2.11 & -0.80 & 1.72 \\\\\n",
" -1.84 & 3.03 & 1.29 \\\\\n",
" -1.57 & 5.25 & 4.30 \\\\\n",
" \\end{pmatrix}\n",
"$$\n",
"$$ (iii)\\quad A = \n",
" \\begin{pmatrix}\n",
" 2 & -1 & 0 \\\\\n",
" -1 & 2 & -1 \\\\\n",
" 0 & -1 & 2 \\\\\n",
" \\end{pmatrix}\n",
"\\quad\\quad\\quad\\quad\n",
"(iv)\\quad A = \n",
" \\begin{pmatrix}\n",
" 4 & 3 & -1 \\\\\n",
" 7 & -2 & 3 \\\\\n",
" 5 & -18 & 13 \\\\\n",
" \\end{pmatrix}\\,.\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"determinant: -7.401486830834414e-17 condition number: 4.065294633547468e+16\n",
"determinant: 0.05886699999999722 condition number: 3218.3325414808673\n",
"determinant: 4.0 condition number: 5.828427124746193\n",
"determinant: 289.99999999999994 condition number: 10.295133664507079\n"
]
}
],
"source": [
"A = np.array([[1.,2.,3.],[2.,3.,4.],[3.,4.,5.]])\n",
"print(\"determinant: \",np.linalg.det(A), \" condition number: \", np.linalg.cond(A))\n",
"# => singular\n",
"\n",
"A = np.array([[2.11,-0.8,1.72],[-1.84,3.03,1.29],[-1.57,5.25,4.30]])\n",
"print(\"determinant: \",np.linalg.det(A), \" condition number: \", np.linalg.cond(A))\n",
"# => ill conditioned\n",
"\n",
"A = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])\n",
"print(\"determinant: \",np.linalg.det(A), \" condition number: \", np.linalg.cond(A))\n",
"# => well conditioned\n",
"\n",
"A = np.array([[4,3,-1],[7,-2,3],[5,-18,3]])\n",
"print(\"determinant: \",np.linalg.det(A), \" condition number: \", np.linalg.cond(A))\n",
"# => well conditioned\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 3: Hilbert matrices ($**$) \n",
"\n",
"The *Hilbert matrix* is a classic example of ill-conditioned matrix:\n",
"\n",
"$$\n",
"A = \n",
" \\begin{pmatrix}\n",
" 1 & 1/2 & 1/3 & \\cdots \\\\\n",
" 1/2 & 1/3 & 1/4 & \\cdots \\\\\n",
" 1/3 & 1/4 & 1/5 & \\cdots \\\\\n",
" \\vdots & \\vdots & \\vdots & \\ddots \\\\\n",
"\\end{pmatrix}\\,.\n",
"$$\n",
"\n",
"Let's consider the linear system $A\\pmb{x}=\\pmb{b}$ where \n",
"$$ b_i = \\sum_{j=1}^n A_{ij},\\quad \\textrm{for}\\quad i=1,2,\\ldots, n.$$\n",
"\n",
"\n",
"- How can you write entry $A_{ij}$ for any $i$ and $j$ ?\n",
"\n",
"\n",
"- Convince yourself by pen and paper that $ \\pmb{x} = \\left[ 1, 1, \\cdots 1\\right]^T$ is the solution of the system.\n",
"\n",
"\n",
"- Write a function that returns $A$ and $b$ for a given $n$.\n",
"\n",
"\n",
"- For a range of $n$, compute the condition number of $A$, solve the linear system and compute the error ($err = \\sum_{i=1}^n \\left|x_{computed, i}-x_{exact, i}\\right|$). What do you observe ?"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n",
"1.1102230246251565e-15\n",
"2.098321516541546e-14\n",
"1.1090017792980689e-12\n",
"9.708789328044531e-12\n",
"1.1745374672855746e-09\n",
"1.864989607192058e-08\n",
"1.3061090481381044e-06\n",
"4.266230945881855e-05\n",
"0.0013671789195044415\n",
"0.024070190939258218\n",
"1.6852308180873334\n",
"32.9242317402154\n",
"53.06593350973159\n",
"49.70691401254196\n",
"35.03838051744423\n",
"102.34414169551886\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\Matt\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:12: LinAlgWarning: Ill-conditioned matrix (rcond=2.55091e-17): result may not be accurate.\n",
" if sys.path[0] == '':\n",
"C:\\Users\\Matt\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:12: LinAlgWarning: Ill-conditioned matrix (rcond=1.1721e-18): result may not be accurate.\n",
" if sys.path[0] == '':\n",
"C:\\Users\\Matt\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:12: LinAlgWarning: Ill-conditioned matrix (rcond=5.97124e-19): result may not be accurate.\n",
" if sys.path[0] == '':\n",
"C:\\Users\\Matt\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:12: LinAlgWarning: Ill-conditioned matrix (rcond=9.14407e-19): result may not be accurate.\n",
" if sys.path[0] == '':\n",
"C:\\Users\\Matt\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:12: LinAlgWarning: Ill-conditioned matrix (rcond=7.09847e-19): result may not be accurate.\n",
" if sys.path[0] == '':\n",
"C:\\Users\\Matt\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:12: LinAlgWarning: Ill-conditioned matrix (rcond=5.32825e-19): result may not be accurate.\n",
" if sys.path[0] == '':\n"
]
}
],
"source": [
"\n",
"def hilbert(n):\n",
" A = np.zeros((n,n))\n",
" b = np.zeros(n)\n",
" for i in range(n):\n",
" for j in range(n):\n",
" A[i,j] = 1./(i+j+1)\n",
" b[i] = np.sum(A[i,:])\n",
" return A,b\n",
"\n",
"for n in range(1,18):\n",
" A, b = hilbert(n)\n",
" x = sl.solve(A,b)\n",
" x_exact = np.ones(n)\n",
" error = np.sum(abs(x-x_exact))\n",
" print(error)\n",
" \n",
"# Error between numerical solution and exact solution is large and varies a lot for n > 14 => ill conditioned"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 4: Vandermonde matrices ($**$) \n",
"\n",
"A *Vandermonde matrix* is defined as follows, for any $\\alpha_1, \\dots, \\alpha_n$ real numbers:\n",
"$$V=\\begin{pmatrix}\n",
"1 & \\alpha_1 & {\\alpha_1}^2 & \\dots & {\\alpha_1}^{n-1}\\\\\n",
"1 & \\alpha_2 & {\\alpha_2}^2 & \\dots & {\\alpha_2}^{n-1}\\\\\n",
"1 & \\alpha_3 & {\\alpha_3}^2 & \\dots & {\\alpha_3}^{n-1}\\\\\n",
"\\vdots & \\vdots & \\vdots & &\\vdots \\\\\n",
"1 & \\alpha_n & {\\alpha_n}^2 & \\dots & {\\alpha_n}^{n-1}\\\\\n",
"\\end{pmatrix}$$\n",
"\n",
"\n",
"- Write a function that takes a real number $\\alpha$ and an integer $n$ as input, and returns a **vector** $v = \\left(1, \\alpha, \\alpha^2, \\dots, \\alpha^{n-1}\\right)$\n",
"\n",
"\n",
"- Using this function, write a function that takes a vector $a = \\left(\\alpha_1, \\alpha_2, \\dots, \\alpha_n\\right)$ as input and returns the corresponding Vandermonde matrix.\n",
"\n",
"\n",
"- For different sets of randomly chosen $(\\alpha_i)$, compute the determinant of the corresponding Vandermonde matrix. What does it tell us regarding the matrix conditioning ?\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"287.99999999999676"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from pprint import pprint\n",
"\n",
"def vdm_row(alpha,n):\n",
" row = np.zeros(n)\n",
" cur_alpha = 1\n",
" for i in range(n):\n",
" row[i] = cur_alpha\n",
" cur_alpha *= alpha\n",
" return row\n",
"\n",
"def vdm(alpha_vec):\n",
" n = alpha_vec.size\n",
" A = np.zeros((n,n))\n",
" for i in range(n):\n",
" A[i,:] = vdm_row(alpha_vec[i],n)\n",
" return A\n",
"\n",
"alphas = np.array([1,2,3, 4, 5])\n",
"V = vdm(alphas)\n",
"\n",
"sl.det(V)\n",
"\n",
"# det is large compared to matruix entries = not very well conditioned"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5: LU solve ($**$) \n",
"\n",
"Write a function that solves a linear system $A\\pmb{x}=\\pmb{b}$ using the LU decomposition method.\n",
"\n",
"Hint: you can re-use the function you have written in lecture 6, or use the built-in function `scipy.linalg.lu` to compute the LU decomposition. Write code that performs the forward substitution and backward substitution. Compare your result to the one given by `scipy.linalg.solve`.\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def LU_dec(A):\n",
" # upper triangular matrix contains gaussian elimination result\n",
" # we won't change A in-place but create a local copy\n",
" A = A.copy()\n",
" m, n = A.shape\n",
" # lower triangular matrix has identity diagonal and scaling factors\n",
" L = np.identity(n)\n",
" # Loop over each pivot row.\n",
" for k in range(n):\n",
" # Loop over each equation below the pivot row.\n",
" for i in range(k+1, n):\n",
" # Define the scaling factor outside the innermost\n",
" # loop otherwise its value gets changed.\n",
" s = (A[i,k]/A[k,k])\n",
" for j in range(k, n):\n",
" A[i,j] = A[i,j] - s*A[k,j]\n",
" # scaling factors make up the lower matrix \n",
" L[i,k] = s\n",
" # A now is the upper triangular matrix U\n",
" return L, A\n",
"\n",
"# This function assumes that A is already an upper triangular matrix.\n",
"def backward_substitution(A, b):\n",
" n = np.size(b)\n",
" \n",
" x = np.zeros(n)\n",
" for k in range(n-1, -1, -1):\n",
" s = 0.\n",
" for j in range(k+1, n):\n",
" s = s + A[k, j]*x[j]\n",
" x[k] = (b[k] - s)/A[k, k]\n",
" \n",
" return x\n",
"\n",
"\n",
"# This function assumes that A is already a lower triangular matrix.\n",
"def forward_substitution(A, b):\n",
" n = np.size(b)\n",
" \n",
" x = np.zeros(n)\n",
" for k in range(n):\n",
" s = 0.\n",
" for j in range(k):\n",
" s = s + A[k, j]*x[j]\n",
" x[k] = (b[k] - s)/A[k, k]\n",
" \n",
" return x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Reminder*: if $A = LU$ then:\n",
"\n",
"$$ A {\\bf x} = {\\bf b} \\Leftrightarrow LU{\\bf x} = {\\bf b} \\Leftrightarrow L(U{\\bf x}) = {\\bf b} \\Leftrightarrow \\left\\{\\begin{align*}L{\\bf y} = {\\bf b} \\quad (1) \\\\ U{\\bf x} = {\\bf y} \\quad (2) \\end{align*} \\right. $$\n",
"\n",
"\n",
"Hence the algorithm to solve $A {\\bf x} = {\\bf b}$ involves the following steps: \n",
"\n",
"- find the LU decomposition of $A$\n",
"\n",
"\n",
"- solve (1) with forward substitution\n",
"\n",
"\n",
"- then solve (2) with backward substitution"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"True\n"
]
}
],
"source": [
"def LU_solve(A,b):\n",
" \n",
" L,U = LU_dec(A)\n",
" y = forward_substitution(L,b)\n",
" x = backward_substitution(U,y)\n",
" \n",
" return x\n",
"\n",
"\n",
"A = np.array([[1.,0.,3.,7.],[2.,1.,0.,4.],[5.,4.,1.,-2.],[4.,1.,6.,2.]])\n",
"b = np.array([1.,2.,-3.,2.])\n",
"\n",
"# Always check that our implemtation produces the right result\n",
"print(np.allclose(sl.solve(A,b), LU_solve(A,b)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6: Gauss-Seidel relaxation ($***$)\n",
"\n",
"Convergence of the Gauss-Seidel method can be improved by a technique known\n",
"as relaxation. The idea is to take the new value of x i as a weighted average of its previous value and the value predicted by the regular Gauss-Seidel iteration. \n",
"\n",
"The corresponding formula for the $k^{th}$ iteration of the algorithm and the $i^{th}$ row is:\n",
"\n",
"$$x_i^{(k)} = \\frac{\\omega}{A_{ii}}\\left(b_i- \\sum_{\\substack{j=1}}^{i-1} A_{ij}x_j^{(k)} - \\sum_{\\substack{j=i+1}}^n A_{ij}x_j^{(k-1)}\\right) + (1-\\omega)x_i^{(k-1)},\\quad i=1,2,\\ldots, n.$$\n",
"\n",
"where the weight $\\omega$ is called the relaxation factor and is usually positive.\n",
"\n",
"\n",
"- What does the algorithm give for $\\omega = 0$; for $\\omega = 1$; for $0 < \\omega < 1$? NB. When $\\omega > 1$, the method is called \"over-relaxation\".\n",
"\n",
"\n",
"- Write a function that solves a system with the relaxed Gauss-Seidel's algorithm, for a given $\\omega$.\n",
"\n",
"\n",
"- Use this function to solve the system from Lecture 7, for different values of $\\omega$. How many iterations are necessary to reach a tolerance of $10^{-6}$ for each value of $\\omega$ ?\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"For omega=0.50 48 relaxed GS iterations are needed to reach a tolerance of 1e-6\n",
"For omega=0.90 18 relaxed GS iterations are needed to reach a tolerance of 1e-6\n",
"For omega=0.95 16 relaxed GS iterations are needed to reach a tolerance of 1e-6\n",
"For omega=1.00 14 relaxed GS iterations are needed to reach a tolerance of 1e-6\n",
"For omega=1.10 12 relaxed GS iterations are needed to reach a tolerance of 1e-6\n",
"For omega=1.50 29 relaxed GS iterations are needed to reach a tolerance of 1e-6\n"
]
}
],
"source": [
"A = np.array([[10., 2., 3., 5.],\n",
" [1., 14., 6., 2.],\n",
" [-1., 4., 16.,-4],\n",
" [5. ,4. ,3. ,11.]])\n",
"b = np.array([1., 2., 3., 4.])\n",
"\n",
"def gauss_seidel_rel_V1(A, b, omega, maxit=500, tol=1.e-6):\n",
" m, n = A.shape\n",
" x = np.zeros_like(b)\n",
" for k in range(maxit):\n",
" for i in range(m):\n",
" x[i] = omega/A[i,i] * (b[i] - np.dot(A[i,:i], x[:i]) - np.dot(A[i,i+1:], x[i+1:])) + (1-omega)*x[i]\n",
" residual = sl.norm(A@x - b)\n",
"# print(\"iteration: %d residual: %e\" %(k,residual))\n",
" if (residual < tol): break\n",
" \n",
" return x,k\n",
"\n",
"# for omega = 0: nothing happens and the system is not solved\n",
"# for omega = 0: standard Gauss-Seidel\n",
"# for 0 < omega < 1: the new x is taken as an average of the old x and the new one. \n",
"# Instead of moving directly from x^(k) to the x^(k+1) given by GS, you stop somewhere in the middle\n",
"# ==> usually slows-down the convergence\n",
"\n",
"for omega in [0.5, 0.9, 0.95, 1, 1.1, 1.5]:\n",
" x,i = gauss_seidel_rel_V1(A, b, omega)\n",
" print(\"For omega=%1.2f %d relaxed GS iterations are needed to reach a tolerance of 1e-6\" % (omega,i))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\omega$ cannot be determined beforehand for an arbitrary system. \n",
"However, an estimate can be computed during run time. \n",
"\n",
"Let $\\Delta_x^{(k)} = | x^{(k)} - x^{(k-1)} |$ be the magnitude of the change in x during the $k^{th}$ iteration. \n",
"If $k$ is sufficiently large (say $k \\geq 5$), it can be shown that an approximation of the optimal value of \\omega is:\n",
"$$\n",
"\\omega_{opt} \\approx \\frac{2}{1+\\sqrt{1-\\Delta x^{(k+1)} / \\Delta x^{(k)}}} \\,.\n",
"$$\n",
"\n",
"The relaxed Gauss-Seidel algorithm can be summarised as follows: \n",
"Carry out $k$ iterations with $\\omega = 1$ (usually $k=10$ for big systems) \n",
"Record \t$\\Delta x^{(k)}$ \n",
"Perform an additional iteration \n",
"Record \t$\\Delta x^{(k+1)}$ \n",
"Compute $\\omega_{opt}$ \n",
"Perform all subsequent iterations with $\\omega = \\omega_{opt}$\n",
"\n",
"\n",
" - Modify previous function to compute automatically the relaxation parameter $\\omega$. Compute $\\omega_{opt}$ after $k=6$ iterations as the system is small.\n",
" - Solve the previous system with this new function. What is the value of $\\omega$ ? How many iterations are necessary to reach a tolerance of $10^{-6}$ ?\n",
" \n",
" \n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Omega given by the formula after 6 iterations: 1.100. 12 relaxed GS iterations are needed to reach a tolerance of 1e-6\n"
]
}
],
"source": [
"from math import sqrt\n",
"\n",
"def gauss_seidel_rel_V2(A, b, k_relax=10, omega_ini=1., maxit=500, tol=1.e-6):\n",
" m, n = A.shape\n",
" x = np.zeros_like(b)\n",
" omega = omega_ini\n",
" for k in range(maxit):\n",
" for i in range(m):\n",
" x[i] = omega/A[i,i] * (b[i] - np.dot(A[i,:i], x[:i]) - np.dot(A[i,i+1:], x[i+1:])) + (1-omega)*x[i]\n",
" residual = sl.norm(A@x - b)\n",
"# print(\"iteration: %d residual: %e\" %(k,residual))\n",
" if (residual < tol): break\n",
" \n",
" if (k == k_relax): \n",
" res_prev = residual\n",
" if (k == k_relax+1): \n",
" res = residual\n",
" omega = 2/(1+sqrt(1-res/res_prev))\n",
" return x,k,omega\n",
"\n",
"k = 6\n",
"x,i,omega = gauss_seidel_rel_V2(A, b, k_relax=k)\n",
"print(\"Omega given by the formula after %d iterations: %1.3f. %d relaxed GS iterations are needed to reach a tolerance of 1e-6\" % (k,omega,i))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### A bigger example\n",
"\n",
"Let's consider $A\\pmb{x}=\\pmb{b}$ where:\n",
"\n",
"$$\n",
"A = \\begin{pmatrix}\n",
"5 & -2 & 0 & 0 & \\cdots & 0 \\\\\n",
"-2 & 5 & -2 & 0 & \\cdots & 0 \\\\\n",
"0 & -2 & 5 & -2 & \\cdots & 0 \\\\\n",
"\\vdots & & & \\ddots & & \\vdots \\\\ \n",
" & & & & 5 & -2 \\\\\n",
"0 & \\cdots & & & -2 & 5 \\\\ \n",
"\\end{pmatrix}\n",
"$$\n",
"and\n",
"$$\n",
"b = \\left(0, 0, \\cdots 0, 1000 \\right)^T\n",
"$$\n",
"\n",
" - Solve $A\\pmb{x}=\\pmb{b}$ using the relaxed Gauss-Seidel algorithm for $n=3000$. Compare the number of iterations with the algorithm without relaxation."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Without relaxation (omega=1.000) 44 iterations are required to reach a tolerance of 1e-6\n",
"With relaxation (omega=1.250) 31 iterations are required to reach a tolerance of 1e-6\n"
]
}
],
"source": [
"n = 3000\n",
"A_big = 5*np.eye(n)\n",
"for i in range(n-1):\n",
" A_big[i,i+1] = -2.\n",
" A_big[i+1,i] = -2.\n",
"b_big = np.zeros(n)\n",
"b_big[n-1] = 1000.\n",
"\n",
"x,i,omega = gauss_seidel_rel_V2(A_big, b_big, k_relax=-2)\n",
"print(\"Without relaxation (omega=%1.3f) %d iterations are required to reach a tolerance of 1e-6\" % (omega,i))\n",
"\n",
"x,i,omega = gauss_seidel_rel_V2(A_big, b_big, k_relax=10)\n",
"print(\"With relaxation (omega=%1.3f) %d iterations are required to reach a tolerance of 1e-6\" % (omega,i))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"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.7.3"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": false,
"sideBar": false,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": false,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}