{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise: Gauss-Seidel Iteration\n",
"\n",
"This exercise focuses on solving the Poisson equation in one dimension with homogeneous Dirichlet conditions:\n",
"$$\n",
"- u''(x) = b(x), \\quad u(0) = u(1) = 0, \\quad b(x) := \\exp(x).\n",
"$$\n",
"A discretization of this equation on a uniform grid using the finite difference method leads to the following linear system:\n",
"$$\n",
"\\frac{1}{h^2}\n",
"\\begin{pmatrix}\n",
" 2 & -1 \\\\\n",
" -1 & 2 & -1 \\\\\n",
" & -1 & 2 & -1 \\\\\n",
" & & \\ddots & \\ddots & \\ddots & \\\\\n",
" & & & -1 & 2 & -1 \\\\\n",
" & & & & -1 & 2 \\\\\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
" u_1 \\\\\n",
" u_2 \\\\\n",
" u_3 \\\\\n",
" \\vdots \\\\\n",
" u_{n-2} \\\\\n",
" u_{n-1}\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
" b(x_1) \\\\\n",
" b(x_2) \\\\\n",
" b(x_3) \\\\\n",
" \\vdots \\\\\n",
" b(x_{n-2}) \\\\\n",
" b(x_{n-1})\n",
"\\end{pmatrix},\n",
"\\quad\n",
"h := \\frac{1}{n},\n",
"\\quad\n",
"x_i := ih.\n",
"$$\n",
"where $h$ is the spacing between the discretization points, and $(u_1, u_2, \\cdots, u_{n-1})$ are the sought values of the unknown function $u$ at the interior points of the domain $[0, 1]$.\n",
"\n",
"1. Compute the solution of the system for $n = 50$ using Julia's `\\` method.\n",
"\n",
" > **Hint**: To build the matrix of the linear system, one can use `LinearAlgebra.SymTridiagonal` and the `fill` function.\n",
"\n",
"2. Implement the Gauss-Seidel method (or its generalization, the relaxation method) to solve the linear system for $n = 50$, this time without using Julia's `\\` and `inv` functions or any software library. Initialize the iteration with $\\mathbf x_0 = \\mathbf 0$.\n",
"\n",
"3. Verify on a graph that the obtained approximate solution is close to the exact solution, which is given by $$u(x) = \\exp(x) - 1 - (\\exp(1) - 1)x.$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise: steepest descent algorithm\n",
"\n",
"Consider the linear system:\n",
"\n",
"$$\n",
"\\mathbf{A}\\mathbf{x} :=\n",
"\\begin{pmatrix}\n",
"3 & 1 \\\\\n",
"1 & 3\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"x_1 \\\\\n",
"x_2\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"1 \\\\\n",
"1\n",
"\\end{pmatrix} =: \\mathbf{b}.\n",
"$$\n",
"\n",
"1. Show that $\\mathbf{A}$ is positive definite.\n",
" \n",
"2. Draw the contour lines of the function:\n",
"\n",
"$$\n",
"f(\\mathbf{x}) = \\frac{1}{2} \\mathbf{x}^T \\mathbf{A} \\mathbf{x} - \\mathbf{b}^T \\mathbf{x}.\n",
"$$\n",
"\n",
"3. Plot the contour lines of $f$ in Julia using the function `contourf` from the package `Plots`.\n",
"\n",
"4. Estimate the number $K$ of iterations of the steepest descent algorithm required to guarantee that $E_K \\leq 10^{-8}$, starting from the vector $\\mathbf{x}^{(0)} = (2~3)^T$.\n",
"\n",
"5. Implement the steepest descent method to find the solution to the given linear system, and plot the iterates as linked dots over the filled contour of $f$.\n",
"\n",
"6. Plot the error $E_k$ as a function of the iteration index, using a linear scale for the $x$ axis and a logarithmic scale for the $y$ axis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise: Conjugate Gradient method\n",
"\n",
"Consider a linear system of the form\n",
"$$\n",
"\\tag{*} \\mathsf A \\mathbf x = \\mathbf b,\n",
"$$\n",
"where $\\mathsf A \\in \\mathbb R^{n\\times n}$ is a positive definite symmetric matrix. The goal of this exercise is to theoretically derive and then implement the conjugate gradient method. Let $\\mathbf x_*$ be the exact solution of the system.\n",
"\n",
"1. **Theoretical Derivation of the Method**.\n",
" Assuming a vector $\\mathbf x_0$ is given, define, for $k \\in \\mathbb N$,\n",
" the Krylov space\n",
" $$\n",
" \\mathcal K_k := \\text{Span} \\Bigl\\{ \\mathbf r_0, \\mathsf A \\mathbf r_0, \\mathsf A^2 \\mathbf r_0, \\dotsc, \\mathsf A^{k-1} \\mathbf r_0 \\Bigr\\},\n",
" \\qquad \\mathbf r_0 := \\mathsf A \\mathbf x_0 - \\mathbf b.\n",
" $$\n",
" Let $\\mathbf x_k$ be the minimizer of $f$ in the affine space $\\mathbf x_0 + \\mathcal K_k$,\n",
" where $f$ is the function\n",
" $$\n",
" f(\\mathbf x) := \\frac{1}{2} \\mathbf x^T \\mathsf A \\mathbf x - \\mathbf b^T \\mathbf x.\n",
" $$\n",
" The goal of this exercise is to arrive at a simple recurrence formula to construct $\\mathbf x_k$ iteratively.\n",
"\n",
" - Show that there exists $k_* \\in \\{0, \\dotsc, n\\}$ such that $\\mathbf x_k = \\mathbf x_*$ for all $k \\geq k_*$.\n",
"\n",
" *Hint 1*: Note that $\\mathsf A \\mathbf x = \\mathbf b$ if and only if $\\mathsf A (\\mathbf x - \\mathbf x_0) = - \\mathbf r_0$.\n",
"\n",
" *Hint 2*: Show that there exists $k_* \\in \\{0, \\dotsc, n\\}$ such that $\\mathsf A \\mathcal K_{k_*} = \\mathcal K_{k_*}$,\n",
" and in particular, $\\mathbf r_0 \\in \\mathsf A \\mathcal K_{k_*}$.\n",
"\n",
" - Show that\n",
" $f(\\mathbf x) = \\frac{1}{2} \\lVert \\mathbf x - \\mathbf x_* \\rVert_{\\mathsf A}^2 + f(\\mathbf x_*)$,\n",
" where $\\langle \\mathbf x, \\mathbf y \\rangle_{\\mathsf A} := \\mathbf x^T \\mathsf A \\mathbf y$\n",
" and $\\lVert \\cdot \\rVert_{\\mathsf A}$ is the associated norm.\n",
"\n",
" *Note*: Notice that this point implies that $\\mathbf x_k - \\mathbf x_0$ is the $\\mathsf A$-orthogonal projection of $\\mathbf x_* - \\mathbf x_0$ onto $\\mathcal K_k$.\n",
"\n",
" - Let $\\mathbf r_k := \\mathsf A \\mathbf x_k - \\mathbf b$.\n",
" Show that $\\mathbf r_k = \\mathbf r_0 + \\mathsf A (\\mathbf x_k - \\mathbf x_0)$,\n",
" and deduce that $\\mathbf r_k \\in \\mathcal K_{k+1}$.\n",
"\n",
" - Show that $\\langle \\mathbf r_{k}, \\mathbf v \\rangle = 0$ for all $\\mathbf v \\in \\mathcal K_{k}$.\n",
"\n",
" *Hint*: Note that $\\mathbf r_k = \\nabla f(\\mathbf x_k)$.\n",
"\n",
" - Deduce from the two previous points that if $\\mathbf r_k \\neq \\mathbf 0$,\n",
" then $\\mathbf r_k \\in \\mathcal K_{k+1} \\setminus \\mathcal K_k$.\n",
"\n",
" - Also deduce that $\\langle \\mathbf r_{k}, \\mathbf v \\rangle_{\\mathsf A} = 0$ for all $\\mathbf v \\in \\mathcal K_{k-1}$,\n",
" for $k \\geq 1$.\n",
"\n",
" *Hint*: Use the definition of the Krylov spaces $\\mathcal K_k$.\n",
"\n",
" - Conclude that if $\\mathbf r_k \\neq \\mathbf 0$,\n",
" then the vector $\\mathbf d_k$ defined by the recurrence relation\n",
" $$\n",
" \\tag{CG1}\n",
" \\mathbf d_{k} = \\mathbf r_{k} - \\frac{\\mathbf r_{k}^T \\mathsf A \\mathbf d_{k-1}}{\\mathbf d_{k-1}^T \\mathsf A \\mathbf d_{k-1}} \\mathbf d_{k-1},\n",
" \\qquad \\mathbf d_0 = \\mathbf r_0\n",
" $$\n",
" \n",
" is the unique (up to a multiplicative constant) element of $\\mathcal K_{k+1}$\n",
" that is $\\mathsf A$-orthogonal to $\\mathcal K_k$.\n",
"\n",
" - Show that $\\mathbf x_{k+1} - \\mathbf x_{k}$ is $\\mathsf A$-orthogonal to $\\mathcal K_k$,\n",
" hence $\\mathbf x_{k+1} = \\mathbf x_{k} - \\omega_k \\mathbf d_k$ for some $\\omega_k \\in \\mathbb R$ when $\\mathbf r_k \\neq 0$.\n",
"\n",
" *Hint*: Note that $\\langle \\mathbf x_{k+1} - \\mathbf x_k, \\mathbf v \\rangle_{\\mathsf A} = \\langle \\mathbf r_{k+1} - \\mathbf r_k, \\mathbf v \\rangle$,\n",
" and use one of the results proven above.\n",
"\n",
" - Finally, show that\n",
" $$\n",
" \\tag{CG2}\n",
" \\mathbf x_{k+1} = \\mathbf x_k - \\frac{\\mathbf r_k^T \\mathbf d_k}{\\mathbf d_k^T \\mathsf A \\mathbf d_k} \\mathbf d_k.\n",
" $$\n",
" \n",
" *Hint*: Find the parameter $\\omega_k$ from the previous point by minimizing $f(\\mathbf x_k - \\omega_k \\mathbf d_k)$.\n",
"\n",
"1. **Implementation of the Method**.\n",
" The relations (CG1) and (CG2) allow for constructing vectors $(\\mathbf x_k)$ recursively,\n",
" with a cost per iteration independent of $k$; this is the **conjugate gradient method**.\n",
"\n",
" - Based on these equations,\n",
" implement a Julia function to solve the system (*).\n",
"\n",
" - Test your implementation on the matrices $\\mathsf A$ and $\\mathbf b$ in the code below,\n",
" which correspond to a finite difference discretization of the Poisson equation on the domain $\\Omega = (0, 2) \\times (0, 1)$,\n",
" with homogeneous Dirichlet boundary condition:\n",
" $$\n",
" \\begin{aligned}\n",
" - \\Delta u(x, y) &= b(x, y), \\qquad (x, y) \\in \\Omega, \\\\\n",
" u(x, y) &= 0, \\quad \\qquad \\quad (x, y) \\in \\partial \\Omega.\n",
" \\end{aligned}\n",
" $$\n",
" The right-hand side considered is\n",
" $$\n",
" b(x, y) = \\sin(4\\pi x) + \\sin(2\\pi y).\n",
" $$\n",
" Use a vector of 1 as the initial vector $\\mathbf x_0$ which can be built using the `ones` function.\n",
"\n",
" - Create a graph illustrating the evolution of the residual $\\lVert \\mathbf r_k \\rVert$, the error $\\lVert \\mathbf x_k - \\mathbf x_* \\rVert$,\n",
" and the error $\\lVert \\mathbf x_k - \\mathbf x_* \\rVert_{\\mathsf A}$ as a function of $k$.\n",
" Choose appropriate scales for the graph (linear-linear, log-log, linear-log, or log-linear)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"vscode": {
"languageId": "julia"
}
},
"outputs": [],
"source": [
"using LaTeXStrings\n",
"using LinearAlgebra\n",
"using Plots\n",
"import SparseArrays\n",
"\n",
"# Domain size\n",
"Lx, Ly = 2, 1\n",
"\n",
"# Number of discretization points along x and y, including the boundary points\n",
"nx, ny = 101, 101\n",
"\n",
"function discretize(nx, ny)\n",
" hx, hy = Lx/(nx - 1), Ly/(ny - 1)\n",
" Dxx = (1/hx^2) * SymTridiagonal(2ones(nx-2), -ones(nx-3))\n",
" Dyy = (1/hy^2) * SymTridiagonal(2ones(ny-2), -ones(ny-3))\n",
" A = kron(Dxx, I(ny-2)) + kron(I(nx-2), Dyy)\n",
" xgrid = Lx/(nx-1) * (1:nx-2)\n",
" ygrid = Ly/(ny-1) * (1:ny-2)\n",
" x_2d = reshape([x for y in ygrid, x in xgrid], (nx-2)*(ny-2))\n",
" y_2d = reshape([y for y in ygrid, x in xgrid], (nx-2)*(ny-2))\n",
" b = sin.(4π*x_2d) + sin.(2π*y_2d)\n",
" return SparseArrays.SparseMatrixCSC(A), b\n",
"end\n",
"\n",
"function plot_solution(f)\n",
" f = reshape(f, ny-2, nx-2)\n",
" z = [zeros(nx)'; zeros(ny-2) f zeros(ny-2); zeros(nx)'] # Add boundary\n",
" xgrid = Lx/(nx-1) * (0:nx-1)\n",
" ygrid = Ly/(ny-1) * (0:ny-1)\n",
" heatmap(xgrid, ygrid, z, c=:viridis, levels=50)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"vscode": {
"languageId": "julia"
}
},
"outputs": [],
"source": [
"# Your function should return a list of iterates x_0, x_1, …\n",
"function conjugate_gradients(A, b)\n",
" # Relpace the line below with your code\n",
" [A\\b]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"vscode": {
"languageId": "julia"
}
},
"outputs": [],
"source": [
"A, b = discretize(nx, ny)\n",
"xs = conjugate_gradients(A, b)\n",
"plot_solution(xs[end])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"vscode": {
"languageId": "julia"
}
},
"outputs": [],
"source": [
"# Produce plots of the residual and errors here"
]
}
],
"metadata": {
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 2
}