"If we deliberately choose $f$ to be outside of the solution space, we obtain a poor approximation and our constraints are not met (which is expected)."
"If we deliberately choose $f$ to be outside of the solution space, we obtain a poor approximation and our constraints are not met (which is expected). \n",
"We can check that the Laplacian of computed solution and the Laplacian of the step function (from which the constraints are sampled) do not match"
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$.
Let us use the matrix $C=\begin{pmatrix}D\\B\end{pmatrix}$.
The following holds $Cu = z$, with $z=\begin{pmatrix}v_0\\v_1\\...\\v_{n-1}\\y_0\\y_1\\...\\y_{m-1}\end{pmatrix}$.
Let us suppose there exist non zero $\alpha_i$ s.t.
$\sum_{i=0}^{n-1} \alpha_i C_i = 0$, where $C_i$ is the $i$-th column of $C$.
Looking at the bottom part of $C$, we must have $\alpha_i=0, i\in\{i_1,...,i_m\}$. The $\alpha_i$ are can only be non-zero when the lower part of the vector only has 0s. Given a subset of columns of $D$ will have at least one row with a single non zero coefficient making the associated $\alpha_i = 0$, which then propagates to the rest (since columns have disjoint support), the rank of $C$ is $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).
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 ?
We obtain a solution that is linear except near it boundaries. The aspect of the solution near the boundaries is due to the Dirichlet conditions we impose (i.e. solution should be constant outside the interval displayed). The values we had set as constraints are matched closely.
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).
The analytical solution is a piecewise linear solution. Integrating once the 2-peaked will lead to picewise constant function, and piecewise linear if we integrate twice.
If you use the constraints we impose for the numerical example you can compute the theoretical solution of the problem satisfying the constraints (that was not asked in the question so the above answer was sufficient). If you impose the solution to be constant between 10 and 40, it means its derivative is 0 in this interval. Below 10 the derivative is then -1 and above 40 the derivative is 1.
We can see the solution computed is acceptable as its second derivative matches the original distribution (it is not fully constant between the peaks but close enough). The Dirichlet boundary conditions are once again interfering at $x=0$ and $x=50$, and the derivative not being continuous also has an impact on the constant part.
