# Minimal distance between point and plane

Consider the case that $V = \mathbb R^3$ and $S$ is the plane spanned by the vectors

$$
b_0 = (-1,-1,0), \quad b_1 = (1,0,1).
$$

Let $u = (0,0,3)$ and let us find the closest point $s \in S$ to $u$. 
There are $c_0, c_1 \in \mathbb R$ such that  

$$
s = c_0 b_0 + c_1 b_1 = B c,
$$

where $B$ is the matrix with columns $b_0$ and $b_1$ and $c$ is the vector with entries $c_0$ and $c_1$. Perhaps the most natural way to find the closest point is to solve the least squares problem to minimize $|Bc - u|^2$.


```python
import numpy as np
from scipy import linalg as la
```


```python
u = np.array([0, 0, 3], dtype=float)
b0 = np.array([-1, -1, 0], dtype=float)
b1 = np.array([1, 0, 1], dtype=float)

B = np.stack((b0, b1), axis=-1)
c, _, _, _ = la.lstsq(B, u)
s = B @ c
c
```


```python
from matplotlib import pyplot as plt

import plot_vec3d as p3d
def prettify():
    p3d.plot_plane(s, c[0] * b0)
    p3d.plot_vec(u - s, 'k--', o=s)
    p3d.plot_coeffs(c, b0, b1, 'k:')
    p3d.set_equal_aspect()
    p3d.hide_ticks()
```


```python
plt.figure(dpi=120)
plt.axes(projection='3d')
p3d.plot_vec(u, 'b')
p3d.plot_vec(s, 'g')
p3d.plot_vec(b0, 'r')
p3d.plot_vec(b1, 'r')
prettify()
```

The finite element method leads to problems that are similar to the above distance minimization, but 

* $V$ will then be an infinite dimensional vector space
* The corresponding least squares problem can not be solved directly

However,

* $S$ will still be finite dimensional
* We will be able to form a finite dimensional system of linear equations 

As a predule, let's demonstrate this idea in the context of the above distance minimization.

By the _orthogonality implies minimality_ lemma, the closest point $s$ is characterized by 

\begin{align}\tag{1}
(s,v) = (u, v) \quad \text{for all $v \in S$}.
\end{align}

In particular, equation (1) needs to hold for $v=b_0$ and $v=b_1$.
On the other hand, if (1) holds for $v=b_0$ and $v=b_1$ then it holds for all $v \in S$ by taking linear combinations. Writing again $s = c_0 b_0 + c_1 b_1$, we get the following system of equations for vector $c$ with entries $c_0$ and $c_1$

$$
\sum_{j=0}^1 c_j (b_j, b_i) = (u, b_i), \quad i=0,1.
$$

We write $K$ for the matrix with entries $(b_i, b_j)$ and $F$ for the vector with entries $(u, b_i)$. Then (1) is equivalent with $Kc = F$.


```python
# Assemble matrix K
K = np.zeros((2,2))
for i in range(2):
    for j in range(2):
        K[i, j] = np.dot(B[:, i], B[:, j])

# Assemble vector F
F = np.zeros(2)
for i in range(2):
    F[i] = np.dot(u, B[:, i])
    
c_alt = la.solve(K, F)

# Compare to the solution using least squares
eps = np.finfo(float).eps # Machine epsilon
la.norm(c_alt - c) / la.norm(c) < 10*eps

```

# P1 finite element

## Definition: P1 basis functions on reference domain 
>$$
\psi_0(x) 
= \begin{cases}
1 - x & x \in [0,1]
\\
0 & \text{otherwise},
\end{cases} 
\qquad
\psi_1(x) 
= \begin{cases}
x & x \in [0,1]
\\
0 & \text{otherwise}.
\end{cases}  
$$

Here the reference domain means the unit interval $[0,1]$.
Viewing $\psi_0$ and $\psi_1$ as functions on $[0,1]$, they form a basis of $\mathbb P_1$.
Indeed, the standard basis of $\mathbb P_1$ is obtained as 

$$
1 = \psi_0(x) + \psi_1(x), \quad x = \psi_1(x).
$$



```python
def cutoff(x):
    return (x > 0) * (x < 1)

def psi0(x):
    return cutoff(x) * (1 - x)
def psi1(x):
    return cutoff(x) * x

xs_plot = np.linspace(-0.2,1.2, 200)
plt.plot(xs_plot, psi0(xs_plot), 'r', label='$\psi_0$')
plt.plot(xs_plot, psi1(xs_plot), 'b--', label='$\psi_1$')
plt.legend();
```

## Definition: local basis functions
>Let 
>
>\begin{align*}
0 = x_0 < x_1 < \dots < x_n = 1,
\end{align*}
>
> write $I_e = [x_{e-1}, x_e]$, $e=1,\dots,n$, and consider the affine function
>
>\begin{align*}
\Phi_e : \mathbb R \to \mathbb R, \quad \Phi_e(x) = \frac{x - x_{e-1}}{x_e - x_{e-1}},
\end{align*}
>
> mapping $I_e$ to $[0,1]$. The local basis functions on $I_e$ are
>
>\begin{align*}
\psi_{j,e}(x) = \psi_j (\Phi_e (x)), \quad j = 0,1.
\end{align*}

The linear interpolant $\mathcal I_h u$ of a function $u \in C(I)$, with $I = [0,1]$, can be written on $I_e \subset I$ as 

$$
\mathcal I_h u(x) = u(x_{e-1}) \psi_{0,e}(x) + u(x_{e}) \psi_{1,e}(x).
$$


```python
n = 4
xs = np.linspace(0,1, n+1)

def Phi(psi, e, x):
    return psi((x - xs[e-1]) / (xs[e] - xs[e-1]))

xs_plot = np.linspace(0,1, 200)
for e in range(1, n+1):
    plt.plot(xs_plot, Phi(psi0, e, xs_plot), label=f'$\psi_0^{e}$')
plt.gca().set_xticks(xs)
plt.legend();
```


```python
for e in range(1, n+1):
    plt.plot(xs_plot, Phi(psi1, e, xs_plot), label=f'$\psi_1^{e}$')
plt.gca().set_xticks(xs)
plt.legend();
```


```python
def u(x):
    return x * (1 - x)

def I_loc(u, x):
    out = np.zeros(np.shape(x))
    for e in range(1, n+1):
        out += u(xs[e-1])*Phi(psi0, e, x) + u(xs[e])*Phi(psi1, e, x)
    return out

plt.plot(xs_plot, u(xs_plot), label='u')
plt.plot(xs_plot, I_loc(u, xs_plot), label='$I_h u$')
plt.gca().set_xticks(xs)
plt.grid(axis='x')
plt.legend();
```

Recall that the basis functions  $\phi_i$ for the space $S = P_{h,0}^1$ were defined in the _P1 finite element space_ section of the lecture notes.
They can be written as

$$
\phi_i(x) = \psi_{1,i}(x) + \psi_{0,i+1}(x), \quad i = 1,\dots,n-1.
$$

For $u$ satisfying $u(0) = 0 = u(1)$ there holds

$$
\mathcal I_h u(x) = \sum_{i=1}^{n-1} u(x_i) \phi_i(x).
$$


```python
def phi(i, x):
    return Phi(psi1, i, x) + Phi(psi0, i+1, x)

for i in range(1, n):
    plt.plot(xs_plot, phi(i, xs_plot), label=f'$\phi_{i}$')
plt.gca().set_xticks(xs)
plt.grid(axis='x')
plt.legend();
```


```python
def I(u, x):
    out = np.zeros(np.shape(x))
    for i in range(1, n):
        out += u(xs[i]) * phi(i, x)
    return out

plt.plot(xs_plot, u(xs_plot), label='u')
plt.plot(xs_plot, I(u, xs_plot), label='$I_h u$')
plt.gca().set_xticks(xs)
plt.grid(axis='x')
plt.legend();
```

# Implementation of the P1 finite element method

Let us consider the bilinear and linear forms $a$ and $L$ defined in the _a boundary value problem_ section of the lecture notes, that is,

$$
a(u, v) = \int_0^1 u'(x) v'(x) dx, \quad L(v) = \int_0^1 f(x) v(x) dx.
$$

Following the proof the _Galerkin solution_ theorem in the lecture notes, we define

$$
u_S = \sum_{j=1}^{n-1} U_j \phi_j, 
\quad
K_{ij} = a(\phi_i, \phi_j),
\quad
F_i = L(\phi_i), 
\qquad i,j=1,\dots,n-1,
$$

and write $U$ and $F$ for the vectors with elements $U_i$ and $F_i$, and $K$ for the matrix with elements $K_{ij}$.

The equation for the Galerkin solution 

$$
a(u_S,v) = L(v) \quad \text{for all $v \in S$}.
$$

is equivalent with $K U = F$.

Let us assemble the matrix $K$. As the global basis functions $\phi_j$ are defined piecewisely on intervals $I_e$, we will consider each interval separately. On $I_e$ there holds

$$
\phi_e(x) = \psi_1^e(x), \quad \phi_{e-1}(x) = \psi_0^e(x), \quad \phi_i(x) = 0, \quad i \ne e, e-1.
$$

Differentiating the local basis functions on $I_e$, we get 

$$
\phi_e'(x) = \frac{1}{x_e - x_{e-1}}, \quad \phi_{e-1}'(x) =  \frac{-1}{x_e - x_{e-1}}, \quad \phi_i'(x) = 0, \quad i \ne e, e-1.
$$

In order to handle the first and last interval in the same way as the rest, we introduce the auxiliary basis functions

$$
\phi_0(x) = \psi_{0,1}(x), \quad \phi_n(x) = \psi_{1,n}(x).
$$


```python
plt.plot(xs_plot, Phi(psi0, 1, xs_plot), label=f'$\phi_0$')
plt.plot(xs_plot, Phi(psi1, n, xs_plot), label=f'$\phi_n$')
plt.gca().set_xticks(xs)
plt.grid(axis='x')
plt.legend();
```


```python
def assemble_K_demo(x):
    n = np.size(x) - 1
    K = np.zeros((n+1, n+1))
    for e in range(1, n+1):
        dx = x[e] - x[e-1]
        dphis = [-1 / dx, 1 / dx]
        ix = [e - 1, e]
        for i in range(2):
            for j in range(2):
                K[ix[i], ix[j]] += dphis[i] * dphis[j] * dx
    
    return K 

K = assemble_K_demo(xs)
# Condense K, that is, 
# throw away the values corresponding to the auxiliary functions
K = K[1:-1,1:-1]
print(K)
```

We turn to assembly of the vector $F$. We need to evaluate the integrals $(f, \phi_j)$ numerically. Let us use the trapezium rule separately on each interval $I_e$:

$$
\int_{I_e} f(x) \phi_e(x) dx \approx \frac{x_e - x_{e-1}}{2} f(x_e),
\quad
\int_{I_e} f(x) \phi_{e-1}(x) dx \approx \frac{x_e - x_{e-1}}{2} f(x_{e-1}),
$$

and

$$
\int_{I_e} f(x) \phi_{j}(x) dx = 0, \quad j\ne e, e-1.
$$


```python
def assemble_F_demo(x, f):
    n = np.size(x) - 1
    F = np.zeros(n+1)
    for e in range(1, n+1):
        dx = x[e] - x[e-1]
        F[e] += dx / 2 * f(x[e])
        F[e-1] += dx / 2 * f(x[e-1])
    return F

def f(x):
    return 1

F = assemble_F_demo(xs, f)
# Condense F
F = F[1:-1] 
print(F)
```

## Example: Poisson's equation in 1d with unit load

Let us solve the problem

$$
\begin{cases}
-u'' = 1 & \text{on $(0,1)$},
\\
u(0) = 0 = u(1).
\end{cases}
$$


```python
def assemble_poisson_demo(xs):
    K = assemble_K_demo(xs)
    K = K[1:-1,1:-1]
    F = assemble_F_demo(xs, f)
    F = F[1:-1] 
    return K, F

def solve_poisson_demo(xs):
    K, F = assemble_poisson_demo(xs)
    u = np.zeros(np.size(xs))
    u[1:-1] = la.solve(K, F)
    return u
```

The solution vector `u` contains the coefficients of the Galerkin solution $u_S$ in the global basis $\phi_0, \dots, \phi_n$.
Recall that the global basis functions satisfy 

$$
\phi_j(x_k) = \delta_{jk} = \begin{cases}
1 & j=k,
\\
0 & \text{otherwise}.
\end{cases}
$$

Hence we can interpret the coefficients `u` as the point values of $u_S$ at the nodes $x_j$, $j=0,\dots,n$.


```python
plt.plot(xs, solve_poisson_demo(xs), label='$u_S$ coarse')
xs_fine = np.linspace(0,1)
plt.plot(xs_fine, solve_poisson_demo(xs_fine), label='$u_S$ fine')
plt.gca().set_xticks(xs)
plt.grid(axis='x')
plt.legend();
```

# Sparsity

Most of the elements of $K$ vanish, that is, $K$ is [sparse](https://en.wikipedia.org/wiki/Sparse_matrix). (In fact, in this particular case, $K$ is [tridiagonal](https://en.wikipedia.org/wiki/Tridiagonal_matrix).)
Solving the system $KU = F$ is typically the most costly part of the finite element method, and efficiency of the method boils down to this system being sparse.


```python
K, F = assemble_poisson_demo(xs_fine)
K
```


```python
plt.spy(K);
```

We could (and should) reimplement `assemble_K_demo` using the [sparse matrix package](https://docs.scipy.org/doc/scipy/reference/sparse.html) of SciPy.

# Implementation using Scikit-fem

Let us now assemble the linear system for the same problem by using the Scikit-fem package.


```python
import skfem as fem
from skfem.helpers import dot, grad # helpers make forms look nice
```


```python
@fem.BilinearForm
def a(u, v, w):
    # Return the integrand in the bilinear form a(u, v)
    # using a notation that works in all dimensions
    return dot(grad(u), grad(v)) 

@fem.LinearForm
def L(v, w):
    x = w.x[0] # quadrature points for numerical integration
    # Return the integrand in the linear form L(v) 
    return f(x) * v # use the same f as before

def assemble_poisson(xs):    
    mesh = fem.MeshLine(xs) 
    basis = fem.Basis(mesh, fem.ElementLineP1())
    K = a.assemble(basis)
    F = L.assemble(basis)
    # Empty call to get_dofs() matches all boundary degrees-of-freedom
    boundary_dofs = basis.get_dofs()
    K, F, _, _ = fem.condense(K, F, D=boundary_dofs)
    return K, F
```


```python
K, F = assemble_poisson(xs)
# K is a sparse matrix, convert it to a standard matrix
print(K.toarray())  
print(F)
```

We get the same linear system as before. Let's check that this is also the case when we use a finer mesh.


```python
K1, F1 = assemble_poisson_demo(xs_fine)
K2, F2 = assemble_poisson(xs_fine)
print(la.norm(K1 - K2))
print(la.norm(F1 - F2))
```

We can solve the Poisson problem in even briefer way using Scikit-fem.


```python
from skfem.models.poisson import laplace, unit_load
import skfem.visuals.matplotlib as fem_plt

mesh = fem.MeshLine(xs) 
basis = fem.Basis(mesh, fem.ElementLineP1())
K = laplace.assemble(basis)
F = unit_load.assemble(basis)
fem_plt.plot(mesh, fem.solve(*fem.condense(K, F, D=basis.get_dofs())));
```

# Computational study of convergence

Let $I = [0,1]$ and let $u$ solve the problem 

$$
\begin{cases}
-u'' = f & \text{on $I$},
\\
u(0) = 0 = u(1),
\end{cases}
$$

and consider the corresponding Galerkin solution $u_S$ computed in a mesh

$$
0 = x_0 < x_1 < \dots < x_n = 1.
$$ 

It is typical to think that $u_S$ is parametrized by the mesh constant

$$
h = \max_{e=1,\dots,n} |x_e - x_{e-1}|,
$$

and write $u_h = u_S$. 

We will study computationally how the error $u_h - u$ converges in the space $L^2(I)$ as $h \to 0$.
In order to do this, we need to set up the problem so that we know the exact solution $u$.
A way to achieve this is to choose first $u$ satisfying the boundary conditions and then compute $f$.

We set 

$$
u(x) = \sin(\pi x).
$$

Then $u(0) = u(1) = 0$ and 

$$
-u''(x) = \pi^2 \sin(\pi x) =: f(x). 
$$


```python
def u_exact(x):
    return np.sin(np.pi * x)

def f(x):
    return np.pi**2 * np.sin(np.pi * x)

# We need to redefine this, as we didn't parametrize f
def solve_poisson_demo(xs):
    K = assemble_K_demo(xs)
    K = K[1:-1,1:-1]
    F = assemble_F_demo(xs, f)
    F = F[1:-1] 
    u = np.zeros(np.size(xs))
    u[1:-1] = la.solve(K, F)
    return u
```


```python
ns = [4, 8]
for n in ns:
    xs, h = np.linspace(0,1, n+1, retstep=True)
    u = solve_poisson_demo(xs)
    plt.plot(xs, u, label=f'{h = }')
xs_plot = np.linspace(0, 1)
plt.plot(xs_plot, u_exact(xs_plot), label='exact')
plt.legend();
```

Recalling that on $I_e$

$$
\phi_e(x) = \psi_1^e(x), \quad \phi_{e-1}(x) = \psi_0^e(x), \quad \phi_i(x) = 0, \quad i \ne e, e-1,
$$

we have

$$
\|u_h - u\|^2 
= \sum_{e=1}^n \int_{I_e} |u_h(x) - u(x)|^2 dx
= \sum_{e=1}^n \int_{I_e} |U_{e} \psi_1^e(x) + U_{e-1} \psi_0^e(x) - u(x)|^2 dx.
$$

Evaluating the integral using the trapezium rule yields

$$
\int_{I_e} |U_{e} \psi_1^e(x) + U_{e-1} \psi_0^e(x) - u(x)|^2 dx
\approx \frac{x_e - x_{e-1}}{2} \left( |U_{e-1} - u(x_{e-1})|^2 + |U_{e} - u(x_{e})|^2 \right).
$$

In principle, we could expand the square and compute the terms

$$
U_i U_j \int_{I_e} \psi_1^i(x) \psi_0^j(x) dx, \quad i,j=e,e-1
$$

analytically. In practice, we can use a numerical integration method, that yields an exact result for polynomials of degree 2 or less, to achieve the same accuracy as the analytic computation. Here we will use the trapezium rule for simplicity. 


```python
def error_sq_demo(x, u):
    '''Squared L2 error'''
    n = np.size(x) - 1
    error = 0
    for e in range(1, n+1):
        dx = x[e] - x[e-1]
        error += dx / 2 * ( (u[e-1] - u_exact(x[e-1]))**2 + (u[e] - u_exact(x[e]))**2 )     
    return error

ns = [2**n for n in range(2,6)]
N = len(ns)
errs_demo = np.zeros(N)
hs = np.zeros(N)
for i in range(N):
    n = ns[i]
    xs, h = np.linspace(0,1, n+1, retstep=True)
    u = solve_poisson_demo(xs)
    hs[i] = h
    errs_demo[i] = np.sqrt(error_sq_demo(xs, u))
```


```python
plt.loglog(hs, errs_demo, '-o', label='error')
plt.loglog(hs, hs**2, label='$h^2$')
plt.legend();
```

# Evaluation of error using Scikit-fem


```python
# We need to redefine this, as we didn't parametrize f
@fem.LinearForm
def L(v, w):
    x = w.x[0] # quadrature points for numerical integration
    # Return the integrand in the linear form L(v) 
    return f(x) * v 

@fem.Functional
def error_sq(w):
    x = w.x[0] 
    uh = w['uh'] # uh is given in the assemble call below
    # Return the integrand in the squared L2 error
    return (uh - u_exact(x))**2

def error_poisson(xs):    
    mesh = fem.MeshLine(xs) 
    basis = fem.Basis(mesh, fem.ElementLineP1())
    K = a.assemble(basis)
    F = L.assemble(basis)
    boundary_dofs = basis.get_dofs()
    u = fem.solve(*fem.condense(K, F, D=boundary_dofs))
    return np.sqrt(error_sq.assemble(basis, uh=basis.interpolate(u)))

errs = np.array([
    error_poisson(np.linspace(0,1, n+1)) for n in ns])
```


```python
plt.loglog(hs, errs, 'ob-', label='skfem')
plt.loglog(hs, errs_demo, '.r--', label='demo')
plt.legend();
```

The two errors are slightly different since Scikit-fem is using more accurate numerical integration method than the trapezium rule. 


```python

```

626d3737a962907de8b7a1e6503d71679d717488
ipynb
Jupyter Notebook
implementation.ipynb
uh-comp-methods2/notebooks
ce55d0e3e6b3d23341ca9d54630b37832a11d1d9
CC-BY-4.0
implementation.ipynb
uh-comp-methods2/notebooks
ce55d0e3e6b3d23341ca9d54630b37832a11d1d9
CC-BY-4.0
implementation.ipynb
uh-comp-methods2/notebooks
ce55d0e3e6b3d23341ca9d54630b37832a11d1d9
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES