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/week02/Linear transforms.ipynb

-%% Cell type:code id:e039b976 tags:
+%% Cell type:code id:60552e92 tags:
 
 ``` python
 # Initialize Otter
 import otter
 grader = otter.Notebook("Linear transforms.ipynb")
 ```
 
 %% Cell type:markdown id:73bbe3f2-6e0a-463e-abc6-01213d4e5811 tags:
 
 # Matrix Analysis 2025 - EE312
 ## Week 2  - Linear transforms
 [N. Aspert](https://people.epfl.ch/nicolas.aspert) - [LTS2](https://lts2.epfl.ch)
 
 The first week notebook (introduction to Python, Numpy and Matplotlib) can be used as a help.
 
-## Important
-You need to submit *individually* your answers on moodle before the next exercise session. For the theoretical questions you can either fill the notebook or write it on a separate sheet (if you are not comfortable with Markdown/TeX) and upload a scanned version.
-
 ## Objective
 
 The end goal is to understand purely algebraic, matrix based, view of a few linear transforms. You will use those linear transform to perform some basic time-frequency analysis of signals.
 
 %% Cell type:code id:9e522683-69ce-458f-a78c-9780e07a865e tags:
 
 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 ```
 
 %% Cell type:markdown id:5901baaa-4eb1-446a-8c0b-88369a9e213e tags:
 
 ## Part I - Fourier
 
 %% Cell type:markdown id:7e81811f-e9c9-4d57-bbb9-bee71b7df722 tags:
 
 <!-- BEGIN QUESTION -->
 
 1. Prove that any set of orthogonal vectors $v_i \in \mathbb{C}^N, \, i=1, \ldots , M \leq N$ such that $v_i^H v_j = C \delta_{i,j}$ is linearly independent (where $C$ is some constant).
 
-%% Cell type:markdown id:c6cd60d1 tags:otter_answer_cell
+%% Cell type:markdown id:b6087a12 tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:markdown id:accea561-acef-429f-aa68-cc002cba19b8 tags:
 
 <!-- END QUESTION -->
 
 <!-- BEGIN QUESTION -->
 
 2. Compute $a_r = \sum_{n=0}^{N-1}e^{j2\pi r\frac{n}{N}}$, where $r$ is an integer (discuss the result depending on the value of $r$).
 
-%% Cell type:markdown id:aed5c02d tags:otter_answer_cell
+%% Cell type:markdown id:3f0ff7a0 tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:markdown id:b591574c-ba6d-4a15-8b42-0cb62508da73 tags:
 
 <!-- END QUESTION -->
 
 <!-- BEGIN QUESTION -->
 
 3. Let $v_k \in \mathbb{C}^N$ be such that $v_k[n] = e^{-j 2 \pi \frac{kn}{N}}$, for $k,n = 0, \ldots , N-1$.
 - Prove that these vectors are mutually orthogonal, hence linearly independent.
 - Compute the norm of $v_k$.
 
-%% Cell type:markdown id:db9f3b40 tags:otter_answer_cell
+%% Cell type:markdown id:d8031ee7 tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:markdown id:efc7a6a3-1602-4100-8ec8-d82f9e01ad86 tags:
 
 <!-- END QUESTION -->
 
 4. Implement the function `get_fourier_matrix` that returns a normalized Fourier matrix of size $N\times N$. Do not make use of the builtin DFT/FFT functions in `numpy` or `scipy`. Raise a `ValueError` exception if a ngetive or zero $N$ value is supplied.
 
 %% Cell type:code id:f58a90c7-d4e3-48d9-a21b-42a4d79e2344 tags:otter_answer_cell
 
 ``` python
 def get_fourier_matrix(N):
     ...
 ```
 
-%% Cell type:code id:a6772e5f tags:
+%% Cell type:code id:2283adf4 tags:
 
 ``` python
 grader.check("q4")
 ```
 
 %% Cell type:markdown id:a5f6b819-1e53-456a-bf67-b76258922d6e tags:
 
 Let us now generate two test signals.
 
 The first one $x_1(t)$ is made of four piecewise sinusoids, of different frequencies:
 
 $$
 x_1(t) = \cos(2\pi 5t), 0\le t < 2\\
 x_1(t) = \cos(2\pi 10t), 2\le t < 4\\
 x_1(t) = \cos(2\pi 25t), 4\le t < 6\\
 x_1(t) = \cos(2\pi 50t), 6\le t < 8\\
 $$
 
 %% Cell type:code id:081e39c8-1874-4274-a7dc-59c8bdab84e4 tags:
 
 ``` python
 Fs = 256 # sampling frequency
 t = np.arange(0, Fs*8)/Fs
 x1 = np.zeros(Fs*8)
 x1[0:Fs*2] = np.cos(2*np.pi*5*t[0:Fs*2])
 x1[Fs*2:Fs*4] = np.cos(2*np.pi*10*t[Fs*2:Fs*4])
 x1[Fs*4:Fs*6] = np.cos(2*np.pi*25*t[Fs*4:Fs*6])
 x1[Fs*6:Fs*8] = np.cos(2*np.pi*50*t[Fs*6:Fs*8])
 ```
 
 %% Cell type:markdown id:b7f7475c-bfb0-4e9f-9dd6-ac9eeb5abef1 tags:
 
 The second signal $x_2(t)$ is the sum of the same sinusoids over the complete time interval, with a scaling term s.t. the amplitude of both signals is identical.
 
 %% Cell type:code id:8e0199b4-baf4-4376-a61c-b778aabe59f7 tags:
 
 ``` python
 x2 = 0.25*(np.cos(2*np.pi*5*t) + np.cos(2*np.pi*10*t) + np.cos(2*np.pi*25*t) + np.cos(2*np.pi*50*t))
 ```
 
 %% Cell type:markdown id:c52f3daf-0017-4252-b297-b5d6bb6c8151 tags:
 
 <!-- BEGIN QUESTION -->
 
 5.
 - Display the generated signals using `plt.plot`.
 - Compute their Fourier transforms using the Fourier matrix generate at the previous question.
 - Display the amplitude of their Fourier spectrum. What do you observe ?
 
-%% Cell type:markdown id:8389b07e tags:otter_answer_cell
+%% Cell type:markdown id:cbde7de1 tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:code id:7d4aee96-f61e-4e51-84e8-56b4e357aa0f tags:otter_answer_cell
 
 ``` python
 # plot x1
 ...
 ```
 
 %% Cell type:code id:7a504349-03a3-4efd-bc3b-8249c09453fb tags:otter_answer_cell
 
 ``` python
 # plot x2
 ...
 ```
 
 %% Cell type:code id:7bf1bf4a-04b8-43f6-8128-29647e1f964a tags:otter_answer_cell
 
 ``` python
 # Compute the Fourier transform of x1 and x2
 ...
 ```
 
 %% Cell type:code id:159e693c-4427-4a8b-9f53-42dc63bafbc9 tags:otter_answer_cell
 
 ``` python
 # Plot the spectrum of x1
 ...
 ```
 
 %% Cell type:code id:f2c3061b-900e-44d9-bb73-02c98334a818 tags:otter_answer_cell
 
 ``` python
 # Plot the spectrum of x2
 ...
 ```
 
 %% Cell type:markdown id:19060be9-6494-4b4e-9803-1a4cb6fc36ac tags:
 
 <!-- END QUESTION -->
 
 In order to have a better compromise between time and frequency, the input signal will be split in smaller non-overlapping blocks of length $p$, and we will perform the DFT of each block.
 
 6. Using the `get_fourier_matrix` implemented previously, fill the `get_block_dft_matrix` function below s.t. it returns a $N\times N$ matrix that will perform the block Fourier transform when applied to an input vector. Raise a `ValueError` if $p$ does not divide $N$.
 
 Hint: [numpy.pad](https://numpy.org/doc/stable/reference/generated/numpy.pad.html#numpy.pad) and/or [numpy.kron](https://numpy.org/doc/stable/reference/generated/numpy.kron.html) might prove useful.
 
 %% Cell type:code id:903005b2-e249-48a0-b28b-24daef16cdb7 tags:otter_answer_cell
 
 ``` python
 def get_block_dft_matrix(N, p):
     ...
 ```
 
-%% Cell type:code id:eae601f6 tags:
+%% Cell type:code id:58978e43 tags:
 
 ``` python
 grader.check("q6")
 ```
 
 %% Cell type:markdown id:4e894db3-b5e6-4b44-a7d1-988c59542b74 tags:
 
 <!-- BEGIN QUESTION -->
 
 We will now use this block Fourier transform to how the frequencies of the signal evolve through time.
 
 7.
 - Using the `reshape` and `plt.imshow` functions, display the amplitude of the result of the block Fourier transform of $x_1$ and $x_2$. Is the result improved when compared to the one observed in question 5 ?
 - What is the influence of parameter $p$ ?
 
-%% Cell type:markdown id:655fc363 tags:otter_answer_cell
+%% Cell type:markdown id:6239185e tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:code id:6da33694-7bd6-4257-9a5f-9bb300cd5c69 tags:otter_answer_cell
 
 ``` python
 # Compute the block DFT matrix Wb
 ...
 ```
 
 %% Cell type:code id:e1eb36fe-3079-4feb-86df-10596293e550 tags:otter_answer_cell
 
 ``` python
 # Plot the block DFT of x1
 ...
 ```
 
 %% Cell type:code id:4de97281-75ba-49d0-bc02-2573cf9791ec tags:otter_answer_cell
 
 ``` python
 # Plot the block DFT of x2
 ...
 ```
 
 %% Cell type:markdown id:76cc3589-b7d9-44ec-98ea-07fc6550a53a tags:
 
 <!-- END QUESTION -->
 
 <!-- BEGIN QUESTION -->
 
 8. In a real-world application, would generating a $N\times N$ matrix to perform the block Fourier transform be a good way to implement it ? (Justify)
 
-%% Cell type:markdown id:1ad05860 tags:otter_answer_cell
+%% Cell type:markdown id:ff1a761a tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:markdown id:87372bd1-4055-4742-ab6a-f7b6ba8750af tags:
 
 <!-- END QUESTION -->
 
 ## Part II - Haar transform
 
 In this part we will study another approach to study the time/frequency properties of signals.
 
 Let us consider a vector $x\in\mathbb{R}^N$, with $N$ being even. The single-level Haar transform decomposes $x$ into two vectors $a^1$ and $d^1$ belonging to $\mathbb{R}^{\frac{N}{2}}$.
 
 The coefficients of $a^1$ are defined as follows: $a^1_n = \frac{1}{\sqrt{2}}(x_{2n-1} + x_{2n}), n=1, ..., \frac{N}{2}$. $a^1$ is referred to as the *average coefficients vector*.
 
 The coefficients of $d^1$ are defined as follows: $d^1_n = \frac{1}{\sqrt{2}}(x_{2n-1} - x_{2n}), n=1, ..., \frac{N}{2}$. $d^1$ is referred to as the *detail coefficients vector*.
 
 %% Cell type:markdown id:170a3e62-8b25-482c-a8dd-91988b1d08ed tags:
 
 <!-- BEGIN QUESTION -->
 
 9. Let us represent the single-level Haar transform by a matrix $H_1$ s.t.
 $$
 H_1 x = \begin{pmatrix}a^1 \\ d^1\end{pmatrix}
 $$
 Prove that $H_1$ is orthonormal.
 
-%% Cell type:markdown id:7a80f125 tags:otter_answer_cell
+%% Cell type:markdown id:56fd9eda tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:markdown id:92f1b78b-3017-4129-bc95-f132717bc7b9 tags:
 
 <!-- END QUESTION -->
 
 10. Write a function that returns the single-level Haar transform matrix $H_1$ for a given $N$. Raise a `ValueError` if $N$ is invalid.
 
 %% Cell type:code id:5d123a1c-7f98-4edf-b896-7d773bf48e76 tags:otter_answer_cell
 
 ``` python
 def get_sl_haar_matrix(N):
     ...
 ```
 
-%% Cell type:code id:7fdfd203 tags:
+%% Cell type:code id:86b9b019 tags:
 
 ``` python
 grader.check("q10")
 ```
 
 %% Cell type:markdown id:e2e7a4af-0f4e-4972-a051-468a12967e79 tags:
 
 The multi-level Haar transform is defined by recursively applying the single-level transform **to the average coefficient parts**.
 
 For instance constructing 2-level Haar transform over $N$ points start with the previously defined $H_{1,N}$ matrix of size $N\times N$ and the corresponding $\frac{N}{2}\times\frac{N}{2}$ version denoted by $H_{1,\frac{N}{2}}$.
 
 $H_{1,N}$ can be written as
 $$
 \begin{pmatrix} H_{1, N}^a \\ H_{1,N}^d \end{pmatrix},
 $$
 where $H_{1, N}^a$ and $H_{1,N}^d$ are respectively the average and detail coefficient matrices, both of size $\frac{N}{2}\times N$.
 
 Following these notations, the 2-level Haar transform matrix $H_{2,N}$ can be written as:
 $$
 \begin{pmatrix} H_{1,\frac{N}{2}}H_{1, N}^a \\ H_{1,N}^d \end{pmatrix},
 $$
 
 11. Implement a function that returns the $H_{p,N}$ matrix of size $N\times N$ that performs a $p$-level haar transform. Raise a `ValueError` if the size and the level are incompatible, or if the level is smaller than 1.
 
 %% Cell type:code id:d73eb83b-8d4b-4bf9-9c19-b3d15ee927c1 tags:otter_answer_cell
 
 ``` python
 def get_haar_matrix(N, level):
     ...
 ```
 
-%% Cell type:code id:e6d9e517 tags:
+%% Cell type:code id:c3138843 tags:
 
 ``` python
 grader.check("q11")
 ```
 
 %% Cell type:markdown id:918653e3-22c1-4b25-abbc-2461a6736557 tags:
 
 <!-- BEGIN QUESTION -->
 
 12. Prove that $H_{p,N}$ is orthonormal.
 
-%% Cell type:markdown id:ee4201b9 tags:otter_answer_cell
+%% Cell type:markdown id:02e21ab7 tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
 %% Cell type:markdown id:8224d668-7cdc-4175-a587-0647df59b2b4 tags:
 
 <!-- END QUESTION -->
 
 ### Haar transform visualization
 
 In order to make the visualization of the Haar decomposition easy, we provide you the `plot_haar_coeffs` function that will display the average and detail coefficients of the different levels.
 
 The function takes 2 arguments:
 - the input signal
 - the number of levels
 
 %% Cell type:markdown id:b6cac5b0-a173-499c-ad89-9e4b82166c7a tags:
 
 <!-- BEGIN QUESTION -->
 
 13. Display the Haar transform of $x1$ and $x2$.
 
 %% Cell type:code id:5ca2ec21-222e-430c-a7bd-d2af54e2de47 tags:
 
 ``` python
 from display_helper import plot_haar_coeffs
 ```
 
 %% Cell type:code id:6a946638-18dc-44fa-8ffc-0b44a3f653e1 tags:otter_answer_cell
 
 ``` python
 # display the decomposition of x1
 ...
 ```
 
 %% Cell type:code id:b01b9267-9609-474e-b5e0-8b7adb687164 tags:otter_answer_cell
 
 ``` python
 # display the decomposition of x2
 ...
 ```
 
 %% Cell type:markdown id:81e23ac6-de1a-4cb4-98f2-89545052f9f3 tags:
 
 <!-- END QUESTION -->
 
 ## Part III - Denoising
 
 We will now use the different transforms defined in part I and II to perform denoising.
 Let us create a noisy signal for this purpose.
 
 %% Cell type:code id:e5583ca2-13ea-49f8-baf6-4d6a223f3b2e tags:
 
 ``` python
 angle1 = np.linspace(0, 5*np.pi/2, 300)
 wave1  = np.sin(angle1)
 angle2 = np.linspace(0, 3*np.pi/2, 300)
 wave2  = np.sin(angle2)
 angle3 = np.linspace(np.pi/2, 2*np.pi, 424)
 wave3  = np.sin(angle3)
 wave   = np.append(wave1, wave2)
 wave   = np.append(wave, wave3)
 wave_noisy = wave + 0.2*np.random.normal(0, 1, 1024)
 ```
 
 %% Cell type:code id:d3278fc1-e814-47b9-83f0-61f4fad8d4d7 tags:
 
 ``` python
 plt.plot(wave_noisy, 'r')
 plt.plot(wave)
 ```
 
 %% Cell type:markdown id:33489d2d-b038-4896-87e7-0af6568ad702 tags:
 
 The noise is usually located in the higher frequencies. However, the signal we created is a bit special as it has two discontinuities which also generate high frequencies components (remember the Fourier transform of a rectangle function is a sinc).
 
 %% Cell type:markdown id:fbd9ce45-7214-4ca7-a5e3-cf0c340458a3 tags:
 
 <!-- BEGIN QUESTION -->
 
 14. Implement a function `denoise_signal` that perform denoising of the input signal by using a supplied orthonormal transform matrix, and by setting the transformed coefficients having an amplitude smaller than a given threshold to 0. You might want to use the [numpy.where](https://numpy.org/doc/stable/reference/generated/numpy.where.html) function for this. When denoising using the Haar transform, you can preform the thresholding only on the detail coefficients. Implement the function `denoise_signal_haar` that performs this operation.
 
 NB: The result returned should be real, in order to be displayed.
 
 %% Cell type:code id:cf7270a7-de1d-4c24-913c-916479ee8ab5 tags:otter_answer_cell
 
 ``` python
 def denoise_signal(W, x, threshold=0.1):
     """
     W: NxN input orthonormal matrix (Fourier, block-Fourier or Haar)
     x: input signal (of length N)
     """
     ...
 ```
 
 %% Cell type:code id:c029ab19-11a5-4bc9-b348-31bf8ca07bca tags:otter_answer_cell
 
 ``` python
 def denoise_signal_haar(W, x, threshold=0.1, detail_start_index=256):
     """
     W: NxN input orthonormal matrix (Fourier, block-Fourier or Haar)
     x: input signal (of length N)
     detail_start_index: thresholding is performed on x[detail_start_index:]
     """
     ...
 ```
 
 %% Cell type:code id:1b47fbc5-070b-49c4-a616-d8e753adede1 tags:otter_answer_cell
 
 ``` python
 # Perform denoising with the full Fourier transform and display the result.
 # Make sure you choose a good threshold
 ...
 ```
 
 %% Cell type:code id:d06fa927-ca94-4387-b987-a4e3ffeabeb2 tags:otter_answer_cell
 
 ``` python
 # Perform denoising with the block Fourier transform and display the result
 # Make sure you choose a good threshold and block size
 ...
 ```
 
 %% Cell type:code id:98e6e4a6-a289-4552-b5a9-27f40c40cb40 tags:otter_answer_cell
 
 ``` python
 # Perform denoising with the Haar transform and display the result
 # Make sure you choose a good threshold and an appropriate number of levels
 ...
 ```
 
-%% Cell type:code id:4ab6578e tags:
+%% Cell type:code id:d5cb5707 tags:
 
 ``` python
 grader.check("q14")
 ```
 
 %% Cell type:markdown id:e6b05f2f-34a1-4689-82ef-c3be92c842ba tags:
 
 <!-- END QUESTION -->
 
 <!-- BEGIN QUESTION -->
 
 15. Compare the three denoising methods (Fourier, block Fourier and Haar). Which one performs better (in terms of noise removal but also in terms of discontinuity preservation) ? Was that expected (justify) ?
 
-%% Cell type:markdown id:9a63dd83 tags:otter_answer_cell
+%% Cell type:markdown id:a999fadd tags:otter_answer_cell
 
 _Type your answer here, replacing this text._
 
-%% Cell type:markdown id:43d32956 tags:
+%% Cell type:markdown id:e705a9ec tags:
 
 <!-- END QUESTION -->
 

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 -->
+