From dd8187bba2f01d406ee2c4fd1965b1c431d0e6db Mon Sep 17 00:00:00 2001
From: Nicolas Aspert <nicolas.aspert@epfl.ch>
Date: Tue, 18 Mar 2025 11:33:12 +0100
Subject: [PATCH] add week 5
---
.../Linear systems - Poisson equation.ipynb | 673 ++++++++++++++++++
1 file changed, 673 insertions(+)
create mode 100644 exercises/week05/Linear systems - Poisson equation.ipynb
diff --git a/exercises/week05/Linear systems - Poisson equation.ipynb b/exercises/week05/Linear systems - Poisson equation.ipynb
new file mode 100644
index 0000000..523959f
--- /dev/null
+++ b/exercises/week05/Linear systems - Poisson equation.ipynb
@@ -0,0 +1,673 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "044b8034",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "outputs": [],
+ "source": [
+ "# Initialize Otter\n",
+ "import otter\n",
+ "grader = otter.Notebook(\"Linear systems - Poisson equation.ipynb\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d133bce9-d159-4b1d-a25f-69bed8da493a",
+ "metadata": {},
+ "source": [
+ "# Matrix Analysis 2025 - EE312\n",
+ "\n",
+ "## Week 5 - Discrete Poisson equation \n",
+ "[LTS2](https://lts2.epfl.ch)\n",
+ "\n",
+ "### Objectives\n",
+ "Apply your knowledge about linear systems to Poisson equation resolution."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "654ad4ee-c060-4bc3-b218-3a6709952750",
+ "metadata": {
+ "tags": []
+ },
+ "source": [
+ "## Poisson equation\n",
+ "Let $u,v \\in \\mathbb{R}^n$ represent a physical quantity $f$ and $g: \\mathbb{R} \\mapsto \\mathbb{R}$ sampled at $n$ equi-spaced locations $x_i$, i.e $u_i = f(x_i)$, $v_i = g(x_i)$. \n",
+ "\n",
+ "Let us assume that the underlying continuous object satisfies the Poisson equation: $\\frac{d^2}{dx^2} f (x)= g(x)$ with constraints $f(x_j) = y_j$ for a subset of $m$ locations $j \\in \\{i_1, \\ldots i_m \\}$. \n",
+ "\n",
+ "We will assume that **all the $i_k$ are distinct**.\n",
+ "\n",
+ "This equation governs a number of physical principles, e.g. gravity, fluid mechanics and electrostatics. In the latter we have $\\Delta \\varphi = -\\frac{\\rho}{\\varepsilon}$, where $\\Delta$ is the Laplacian operator $(\\frac{d^2}{dx^2} + \\frac{d^2}{dy^2} + \\frac{d^2}{dz^2})$, $\\varphi$ is the electric potential, $\\rho$ the density of electric charges and $\\varepsilon$ the electric permittivity.\n",
+ "\n",
+ "---"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "06ac5cfa-62ef-4477-a957-e240ddeb85d8",
+ "metadata": {},
+ "source": [
+ "### 1D case"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7f3cb74d-72f0-460b-a978-a47f87de8557",
+ "metadata": {
+ "tags": []
+ },
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9fd89d0b-50df-4711-a32b-ec1026e6b122",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "1. Write down a matrix equation for the discrete version of $\\frac{d^2}{dx^2} f (x)= g(x)$, using the finite-difference approximation of the derivative $\\frac{d^2}{dx^2} f = f(x_{k+1}) - 2f(x_k) +f(x_{k-1})$. For the boundary conditions we will make the assumption that $x_{-1}=x_0$ and $x_{n-1}=x_n$ (also referred to as Dirichlet boundary conditions)\n",
+ "\n",
+ "What is the rank of the corresponding matrix $D$ ? "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d1f8bd46",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "aed76fe2-8cb5-4b28-bd3b-fa0bea113a46",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<Your answers here>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cceac118-c43d-4782-9f7e-4392e5ecc451",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "2. Implement a function that creates the $D$ matrix (also called Laplacian). The [diag](https://numpy.org/doc/stable/reference/generated/numpy.diag.html) function in numpy might be useful."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d41105d8-6933-4ef8-b2f0-c6ac053e1f55",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "def lapl_matrix(N):\n",
+ " ..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "756d9370",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "outputs": [],
+ "source": [
+ "grader.check(\"q2\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4cea41dc-6f12-42f2-8644-b51e4918504d",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "3. Write down a matrix equation for the discrete version of $f(x_j) = y_j$. What is the rank of the corresponding matrix $B$ ?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fc67bfcb",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "418bc83f-7d46-41c2-8a7e-df892857320c",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "4. Implement a function that creates matrix $B$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1d294eda-7edc-42c2-9640-13bec779bd48",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "def b_matrix(N, idx_list): # idx_list is a list of valid indices, e.g. N=5, idx_list=[1,3]\n",
+ " ..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d592616d",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "outputs": [],
+ "source": [
+ "grader.check(\"q4\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b77f215e-e4d0-4976-9ed6-7ee210c689a6",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "3. Write down a matrix equation for the full problem (Hint: Consider the matrix $C=\\begin{pmatrix}D\\\\B\\end{pmatrix}$). Discuss the rank of the matrix and deduce a way to numerically solve the linear system. Implement a function that returns matrix $C$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "70f5637a",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2eaf5a8c-b83a-4f53-9c4e-6192bf4e6fa4",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "def c_matrix(D, B):\n",
+ " ..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d77116d4",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "outputs": [],
+ "source": [
+ "grader.check(\"q4\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "10a3d136-f5f3-4c1b-9c0a-ece58f891dfe",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "5. What explicit formula can you use to compute the pseudo-inverse of $C$ (justify)?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ce8aa19a",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3726eb03-0c82-4aed-b179-77ff947583c3",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "6. Implement a function that return the solution of the problem. You can either use the formula above (you may then use the [linalg.inv](https://numpy.org/doc/stable/reference/generated/numpy.linalg.inv.html) or compute the pseudo-inverse using [linalg.pinv](https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "79f8c6d1-7b92-48b6-bb90-e888aaa41243",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "# v is a vector of size N\n",
+ "# u is a vector of size len(idx_list)\n",
+ "def solve_poisson(N, v, idx_list, u): \n",
+ " ..."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ff6184e1-e308-47d1-8b72-ab4559b48adb",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "7. Let us now consider the solutions of the particular case $g(x) = 0, \\forall x\\in\\mathbb{R}$. What are the analytical solutions that would solve Poisson equation ?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "53e63141",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4a33db39-077d-4775-a55d-ae5bd7400b23",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "Let us verify that our implementation of the solution works with a small example. Fill the values for $u$ s.t. they belong to the solution space."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "96381701-e737-41c1-8642-fa1a5c13924c",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "N = 50\n",
+ "v = np.zeros(N)\n",
+ "idx_list = [10, 20, 30, 40]\n",
+ "u = ...\n",
+ "sol = solve_poisson(N, v, idx_list, u)\n",
+ "plt.plot(sol) # plot the result"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "04e6fc55-6c8a-4dfa-8ca0-9b45ed2f157c",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "9. Are the results conform to what we expected ? What happens near the boundaries ?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e1c9c95f",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ecf46b87-accd-4043-8ab6-06835830699b",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "10. Let us use a step function for $f$ and try to solve the system"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d81b4501-b2a5-4412-a266-d04ffa8dd2c7",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "outputs": [],
+ "source": [
+ "idx_list = [10, 20, 30, 40]\n",
+ "u2 = [1, 1, -1, -1]\n",
+ "sol2 = solve_poisson(N, v, idx_list, u2)\n",
+ "plt.scatter(idx_list, u2, marker='+', color='r')\n",
+ "plt.plot(sol2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4918d31c-40a2-4df3-88ec-d18a80e879e7",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "What do you observe and what is the reason for that ? "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1482c041",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b520b2c6-7c84-4c3e-be16-8eb866648a97",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n",
+ "<!-- BEGIN QUESTION -->\n",
+ "\n",
+ "11. Let us now consider the application of the Poisson equation to electrostatic. Knowing a distribution of electric charges, we can solve it to compute the electric potential. The distribution of electric charges will be very simple, consisting of two punctual charges"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "96932aa1-74aa-4f4e-bbf5-df0f0d85c47c",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "outputs": [],
+ "source": [
+ "N = 50\n",
+ "v3 = np.zeros(N)\n",
+ "v3[10] = 1\n",
+ "v3[40] = 1\n",
+ "plt.plot(v3)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "42861911-53e1-4cb3-8dc7-aa617c01a7e8",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "- What are the analytical solutions for this problem ? (remember the discrete derivative of a step function is a Dirac function, i.e. 0 everywhere except at one point where it is 1). \n",
+ "- Plot the analytical solution\n",
+ "- Compute and plot the numerical solution\n",
+ "- What do you observe ? "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1b19e6d2",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "source": [
+ "_Type your answer here, replacing this text._"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f42fc23b-c412-4791-b5b6-bef54fa3bab3",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f2f9bd66-067e-4590-8088-87b7969621bb",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "# Plot analytical solution\n",
+ "..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "22eae7ab-4c91-4000-b123-7105cc4c023a",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "idx_list = [10, 20, 30, 40]\n",
+ "u3 = [1, 1, 1, 1]\n",
+ "sol3 = solve_poisson(N, v3, idx_list, u3)\n",
+ "plt.scatter(idx_list, u3, marker='+', color='r')\n",
+ "plt.plot(sol3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b0310f39-9bbf-4518-be7d-af0e3f6724a7",
+ "metadata": {
+ "tags": [
+ "otter_answer_cell"
+ ]
+ },
+ "outputs": [],
+ "source": [
+ "# check the 2nd derivative of the solution computed\n",
+ "..."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c0f0755f",
+ "metadata": {
+ "deletable": false,
+ "editable": false
+ },
+ "source": [
+ "<!-- END QUESTION -->\n",
+ "\n"
+ ]
+ }
+ ],
+ "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.11.10"
+ },
+ "otter": {
+ "OK_FORMAT": true,
+ "tests": {
+ "q2": {
+ "name": "q2",
+ "points": null,
+ "suites": [
+ {
+ "cases": [
+ {
+ "code": ">>> def testD(N):\n... D = lapl_matrix(N)\n... np.testing.assert_array_almost_equal(np.sum(D, axis=1), np.zeros(N))\n... np.testing.assert_array_almost_equal(np.diagonal(D, 1), np.ones(N - 1))\n... np.testing.assert_array_almost_equal(np.diagonal(D, -1), np.ones(N - 1))\n... v = -2 * np.ones(N)\n... v[0] = -1\n... v[N - 1] = -1\n... np.testing.assert_array_almost_equal(np.diagonal(D), v)\n>>> testD(5)\n>>> testD(12)\n",
+ "failure_message": "Check your implementation",
+ "hidden": false,
+ "locked": false,
+ "success_message": "Good, results seem correct"
+ }
+ ],
+ "scored": true,
+ "setup": "",
+ "teardown": "",
+ "type": "doctest"
+ }
+ ]
+ },
+ "q4": {
+ "name": "q4",
+ "points": null,
+ "suites": [
+ {
+ "cases": [
+ {
+ "code": ">>> B = b_matrix(5, [1, 3])\n>>> Btest = np.array([[0, 1, 0, 0, 0], [0, 0, 0, 1, 0]])\n>>> np.testing.assert_array_almost_equal(B, Btest)\n",
+ "failure_message": "Check your implementation",
+ "hidden": false,
+ "locked": false,
+ "success_message": "Good, results seem correct"
+ }
+ ],
+ "scored": true,
+ "setup": "",
+ "teardown": "",
+ "type": "doctest"
+ }
+ ]
+ }
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
--
GitLab