environment.yml

@@ -11,4 +11,5 @@ dependencies:
   - scikit-image=0.24.0
   - scikit-learn=1.5.1
   - scipy=1.14.1
-  
\ No newline at end of file
+  - pip:
+    - otter-grader==6.1.0

exercises/week05/Linear systems - Poisson equation.ipynb

+%% Cell type:code id:044b8034 tags:
+
+``` python
+# Initialize Otter
+import otter
+grader = otter.Notebook("Linear systems - Poisson equation.ipynb")
+```
+
+%% Cell type:markdown id:d133bce9-d159-4b1d-a25f-69bed8da493a tags:
+
+# Matrix Analysis 2025 - EE312
+
+## Week 5 - Discrete Poisson equation
+[LTS2](https://lts2.epfl.ch)
+
+### Objectives
+Apply your knowledge about linear systems to Poisson equation resolution.
+
+%% Cell type:markdown id:654ad4ee-c060-4bc3-b218-3a6709952750 tags:
+
+## Poisson equation
+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.
+
+---
+
+%% Cell type:markdown id:06ac5cfa-62ef-4477-a957-e240ddeb85d8 tags:
+
+### 1D case
+
+%% Cell type:code id:7f3cb74d-72f0-460b-a978-a47f87de8557 tags:
+
+``` python
+import numpy as np
+import matplotlib.pyplot as plt
+```
+
+%% Cell type:markdown id:9fd89d0b-50df-4711-a32b-ec1026e6b122 tags:
+
+<!-- BEGIN QUESTION -->
+
+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$ ?
+
+%% Cell type:markdown id:d1f8bd46 tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:markdown id:aed76fe2-8cb5-4b28-bd3b-fa0bea113a46 tags:
+
+<!-- END QUESTION -->
+
+<Your answers here>
+
+%% Cell type:markdown id:cceac118-c43d-4782-9f7e-4392e5ecc451 tags:
+
+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 id:d41105d8-6933-4ef8-b2f0-c6ac053e1f55 tags:otter_answer_cell
+
+``` python
+def lapl_matrix(N):
+    ...
+```
+
+%% Cell type:code id:756d9370 tags:
+
+``` python
+grader.check("q2")
+```
+
+%% Cell type:markdown id:4cea41dc-6f12-42f2-8644-b51e4918504d tags:
+
+<!-- BEGIN QUESTION -->
+
+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 tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:markdown id:418bc83f-7d46-41c2-8a7e-df892857320c tags:
+
+<!-- END QUESTION -->
+
+4. Implement a function that creates matrix $B$
+
+%% Cell type:code id:1d294eda-7edc-42c2-9640-13bec779bd48 tags:otter_answer_cell
+
+``` python
+def b_matrix(N, idx_list): # idx_list is a list of valid indices, e.g. N=5, idx_list=[1,3]
+    ...
+```
+
+%% Cell type:code id:d592616d tags:
+
+``` python
+grader.check("q4")
+```
+
+%% Cell type:markdown id:b77f215e-e4d0-4976-9ed6-7ee210c689a6 tags:
+
+<!-- BEGIN QUESTION -->
+
+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 tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:code id:2eaf5a8c-b83a-4f53-9c4e-6192bf4e6fa4 tags:otter_answer_cell
+
+``` python
+def c_matrix(D, B):
+    ...
+```
+
+%% Cell type:code id:d77116d4 tags:
+
+``` python
+grader.check("q4")
+```
+
+%% Cell type:markdown id:10a3d136-f5f3-4c1b-9c0a-ece58f891dfe tags:
+
+<!-- END QUESTION -->
+
+<!-- BEGIN QUESTION -->
+
+5. What explicit formula can you use to compute the pseudo-inverse of $C$ (justify)?
+
+%% Cell type:markdown id:ce8aa19a tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:markdown id:3726eb03-0c82-4aed-b179-77ff947583c3 tags:
+
+<!-- END QUESTION -->
+
+<!-- BEGIN QUESTION -->
+
+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 id:79f8c6d1-7b92-48b6-bb90-e888aaa41243 tags:otter_answer_cell
+
+``` python
+# v is a vector of size N
+# u is a vector of size len(idx_list)
+def solve_poisson(N, v, idx_list, u):
+    ...
+```
+
+%% Cell type:markdown id:ff6184e1-e308-47d1-8b72-ab4559b48adb tags:
+
+<!-- END QUESTION -->
+
+<!-- BEGIN QUESTION -->
+
+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 tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:markdown id:4a33db39-077d-4775-a55d-ae5bd7400b23 tags:
+
+<!-- END QUESTION -->
+
+<!-- BEGIN QUESTION -->
+
+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 id:96381701-e737-41c1-8642-fa1a5c13924c tags:otter_answer_cell
+
+``` python
+N = 50
+v = np.zeros(N)
+idx_list = [10, 20, 30, 40]
+u = ...
+sol = solve_poisson(N, v, idx_list, u)
+plt.plot(sol) # plot the result
+```
+
+%% Cell type:markdown id:04e6fc55-6c8a-4dfa-8ca0-9b45ed2f157c tags:
+
+<!-- END QUESTION -->
+
+<!-- BEGIN QUESTION -->
+
+9. Are the results conform to what we expected ? What happens near the boundaries ?
+
+%% Cell type:markdown id:e1c9c95f tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:markdown id:ecf46b87-accd-4043-8ab6-06835830699b tags:
+
+<!-- END QUESTION -->
+
+<!-- BEGIN QUESTION -->
+
+10. Let us use a step function for $f$ and try to solve the system
+
+%% Cell type:code id:d81b4501-b2a5-4412-a266-d04ffa8dd2c7 tags:
+
+``` python
+idx_list = [10, 20, 30, 40]
+u2 = [1, 1, -1, -1]
+sol2 = solve_poisson(N, v, idx_list, u2)
+plt.scatter(idx_list, u2, marker='+', color='r')
+plt.plot(sol2)
+```
+
+%% Cell type:markdown id:4918d31c-40a2-4df3-88ec-d18a80e879e7 tags:
+
+What do you observe and what is the reason for that ?
+
+%% Cell type:markdown id:1482c041 tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:markdown id:b520b2c6-7c84-4c3e-be16-8eb866648a97 tags:
+
+<!-- END QUESTION -->
+
+<!-- BEGIN QUESTION -->
+
+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 id:96932aa1-74aa-4f4e-bbf5-df0f0d85c47c tags:
+
+``` python
+N = 50
+v3 = np.zeros(N)
+v3[10] = 1
+v3[40] = 1
+plt.plot(v3)
+```
+
+%% Cell type:markdown id:42861911-53e1-4cb3-8dc7-aa617c01a7e8 tags:
+
+- 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).
+- Plot the analytical solution
+- Compute and plot the numerical solution
+- What do you observe ?
+
+%% Cell type:markdown id:1b19e6d2 tags:otter_answer_cell
+
+_Type your answer here, replacing this text._
+
+%% Cell type:code id:f42fc23b-c412-4791-b5b6-bef54fa3bab3 tags:otter_answer_cell
+
+``` python
+...
+```
+
+%% Cell type:code id:f2f9bd66-067e-4590-8088-87b7969621bb tags:otter_answer_cell
+
+``` python
+# Plot analytical solution
+...
+```
+
+%% Cell type:code id:22eae7ab-4c91-4000-b123-7105cc4c023a tags:otter_answer_cell
+
+``` python
+idx_list = [10, 20, 30, 40]
+u3 = [1, 1, 1, 1]
+sol3 = solve_poisson(N, v3, idx_list, u3)
+plt.scatter(idx_list, u3, marker='+', color='r')
+plt.plot(sol3)
+```
+
+%% Cell type:code id:b0310f39-9bbf-4518-be7d-af0e3f6724a7 tags:otter_answer_cell
+
+``` python
+# check the 2nd derivative of the solution computed
+...
+```
+
+%% Cell type:markdown id:c0f0755f tags:
+
+<!-- END QUESTION -->
+