"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"
"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",
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)$.
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 \}$.
We will assume that **all the $i_k$ are distinct**.
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.
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)
What is the rank of the corresponding matrix $D$ ?
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.
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$.
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).
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 ?
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
- 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).
