{ "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 }