The objective of this article is to present the application of the spectral method to solve parabolic partial differential equations (PDE). The theoretical resolution will be done jointly with a practical study of a typical problem.

Let $a\frac{\partial^2 u}{\partial x^2} + b\frac{\partial^2 u}{\partial xy} + c\frac{\partial^2 u}{\partial y^2} + d\frac{\partial u}{\partial x} + e\frac{\partial u}{\partial y} + fu = g$ a second oder linear partial differential equations with $a,b,c,d,e,f$ constants. Such equation is said to be parabolic when $b^2-4ac=0$ and elliptic, respectively hyperbolic, when $b^2-4ac<0$, respectively $b^2-4ac>0$.

The boundary value problem that we will study is
$$ (P)\begin{cases} \frac{\partial u}{\partial t} - \Delta u = f(x,y,t)\in L^2(\Omega ), \Omega= ]0,1[\times]0,1[\times ]0,T[\cr u(x,0,t)=\frac{x^2}{2}\cr u(x,1,t)=\frac{-(x-1)^2}{2}e^{-t}\cr \frac{\partial u}{\partial n}(0,y,t) = -ye^{-t}\cr \frac{\partial u}{\partial n}(1,y,t) = 1-y\cr u(x,y,0) = \sin(2\pi y)-x^2y+xy-\frac{y}{2}+\frac{x^2}{2} \end{cases} $$
The problem models the heat diffusion over time in a squared homogeneous material, heated via a given heat source $f$ of finite energy. The boundary conditions are of Dirichlet type on vertical sides and Neumann on the horizontal sides.

Finally, the initial condition gives us information on the initial heat distribution on the domain.

The solution is broken down into four steps:

1. Homogenizing the problem, i.e. transforming the initial formulation such that the boundary conditions are equal to $0$,
2. Solving the Spectral problem associated to the Laplacian operator, i.e. finding its eigenvalues and eigenfunctions,
3. Rewriting the homogenized problem into the spectral basis
4. Solving the Ordinary Differential Equation that appeared with the previous step.

Homogenization via Ruled Surface

The boundary conditions are not homogeneous such that we first need to transform the original problem $(P)$ into an homogeneous one. For that, we define $w(x,y,t) = \phi(x,y) + u(x,y,t)$. By looking at the boundary conditions, we can deduce that we are looking for a ruled surface  $\phi(x,y)$ such that $(P)$ is homogeneous.

The Ruled Surface

We are looking $\phi(x,y) = g_1(x) + yg_2(y)$.

Using the initial condition on the boundaries, we obtain:
$$\begin{cases}w(x,0,t) = \phi(x,0) + u(x,0,t) = 0 \Rightarrow \boxed{g_1(x) = -\frac{x^2}{2}}\\ w(x,1,t) = \phi(x,1) + u(x,1,t) = 0 \Rightarrow \boxed{g_2(x) = (x^2 - x + \frac{1}{2})e^{-t}}\end{cases}$$
Therefore, the function $\phi$ that we are looking for is as follows:
$$\boxed{\phi(x,y) = -\frac{x^2}{2} + y(x^2 - x + \frac{1}{2})e^{-t}}$$

Writing the homogeneous problem

We have $w(x,y,t) = \phi(x,y) + u(x,y,t)$, such that, by linearity of the Laplacian operator, we obtain $- \Delta u(x,y,t) = \Delta \phi(x,y) - \Delta w(x,y,t)$.
Similarly, $\frac{\delta u}{\delta t}(x,y,t) = \frac{\delta w}{\delta t}(x,y,t) - \frac{\delta \phi}{\delta t}(x,y,t)$
Therefore, $\frac{\delta w}{\delta t}(x,y,t) - \Delta w(x,y,t) = f(x,y,t) + \Delta \phi(x,y,t) - \frac{\delta \phi}{\delta t}(x,y,t)$.
By definition of $\phi$, we have:
$$\frac{\delta u}{\delta t}(x,y,t) - \Delta \phi = 1-y(1+e^{-t}(x^2-x+\frac{3}{2}))$$
which allows to write the following equivalence:
$$(P) \Leftrightarrow (P_h)\begin{cases} \frac{\partial w}{\partial t} - \Delta w = f(x,y,t) + 1-y(1+e^{-t}(x^2-x+\frac{3}{2})), f\in L^2(\Omega ), \Omega=  ]0,1[\times]0,1[\times ]0,T[\cr w(x,0,t)=0\cr w(x,1,t)=0\cr \frac{\partial w}{\partial n}(0,y,t) = 0\cr \frac{\partial w}{\partial n}(1,y,t) = 0\cr w(x,y,0) = \sin(2\pi y)-x^2y+xy-\frac{y}{2}+\frac{x^2}{2} \end{cases}$$

Solving the spectral problem associated to $(P_h)$

At first, we will focus only on the spatial part of the problem. The Laplacian being a compact linear operator, we are looking for its eignvalues and eigenfunctions, $\lambda_k$ and $w_k$,  solutions to the equation:
$$-\Delta w = \lambda w~~(PS)$$
This solution exists everywhere in $L^2(\Omega)$ thanks to the following spectral decomposition theorem:

The proof of this theorem is quite long and uses the theorem of spectral decomposition of a compact operator in a Hilbert space. This is also the reason why we need the compactness of the inclusion operator and its continuity. Details can be found in [1].

To take into account the Dirichlet conditions on an edge, we choose $V=H_0^1(\Omega)$ and $H=L^2(\Omega)$. A variational approach applied to $-\Delta w = \lambda w$, naturally leads to finding $a(u,v) = \int_{\Omega} \nabla u . \nabla v dx = \lambda \int_{\Omega} uv dx= \lambda <u,v>_{L^2(\Omega)}$. According to Rellich's theorem, $H_0^1(\Omega)$ is indeed compactly included in $L^2(\Omega)$. As the space of test functions $\cal{D}(\Omega)$ is dense in both spaces, we deduce that $H_0^1(\Omega)$ is dense in $L^2(\Omega)$ and we can apply the spectral theorem above, demonstrating the existence of a solution to the eigenvalue problem.

To return to the original problem, a Green formula suffices, taking into account the boundary conditions and using the trace theorem to conclude. We can more simply go through a weak form via the distributions, having already verified the compact and continuous injection:

$$\begin{matrix}~  \int_{\Omega} \nabla u . \nabla v dx = \lambda \int_{\Omega} uv dx\cr \Leftrightarrow < \nabla u , \nabla v > = \lambda < u, v>\cr \Leftrightarrow -< \Delta u , v > = \lambda < u, v>\end{matrix}$$

Hence $-\Delta u = \lambda u$ in the sense of distributions. From the density of $\cal D$ in $H_0^1$ and $L^2$ and by compact canonical injection, we deduce the equivalence in the sense of $L^2$.

Now that we have been convinced of the existence of the solutions of $(PS)$, we need to explain these solutions. To do this, we will consider a separation of variables by setting: $w(x,t)=X(x)Y(y)$.

From the conditions at the edge we can deduce that:
$$\begin{cases}w(x,0) = 0 = X(x)Y(0) \Rightarrow \boxed{Y(0)=0}\\ w(x,1) = 0 = X(x)Y(1)\Rightarrow \boxed{Y(1)=0}\end{cases}$$
$$\frac{\delta w}{\delta x}(0,y) = X_(x)Y(y) \Rightarrow \boxed{\begin{cases} X'(0)=0\\X'(1)=0\end{cases}}$$
Therefore:
$$\begin{matrix}~ (PS) & \Leftrightarrow & -X"(x)Y(y) - X(x)Y"(y) = \lambda X(x)Y(y) \cr & \Leftrightarrow &-\frac{X"(x)}{X(x)} - \frac{Y"(y)}{Y(y)} = \lambda \cr &\Rightarrow & \begin{cases} - \frac{X"(x)}{X(x)}=\alpha~~~(*)\cr - \frac{Y"(y)}{Y(y)}=\beta~~~(**)\cr \alpha+\beta = \lambda \end{cases} \end{matrix}$$
The solution poses no issue since $(*)$ and $(**)$ are second order linear differential equations with constant coefficients. A simple analysis of the characteristic equation gives us the form of the solution and the conditions on $X'$ and $Y$ deduced from the conditions at the edge of the homogeneous problem make it possible to eliminate solutions which are identically zero and therefore not interesting nor from a mathematical point of view and even less from a physical point of view.

Solving $(*)$

$$- \frac{X"(x)}{X(x)}=\alpha \Rightarrow -X"(x) - \alpha X(x) = 0$$
From the characteristic equation we find:
$$ \Delta = -4\alpha$$
We then distinguish 3 cases:

$\alpha > 0$

Thus $\Delta < 0$ and $X$ is of the form $X(x) = A_x \cos(\sqrt \alpha x) + B_x \sin(\sqrt \alpha x)$.
We deduce that $X'(x) = -\sqrt \alpha A_x \sin(\sqrt \alpha x) + \sqrt \alpha B_x \cos(\sqrt \alpha x)$.
From the boundary conditions we deduce that:
$$\begin{cases}X'(0) = \sqrt \alpha B_x = 0 \Rightarrow B_x = 0\\ X'(1) = -\sqrt \alpha A_x \sin(\sqrt \alpha) = 0 \Rightarrow \sin(\sqrt \alpha) = 0\end{cases}$$
We deduce that $\alpha_k = k^2\pi^2, \forall k>0$ and therefore $ X_k(x) = A_x\cos(k\pi x), \forall k>0$.

$\alpha = 0$

So $\Delta = 0$ and $X$ is of the form $X(x) = A_xx + B_x$. We deduce that $X'(x) = -\sqrt \alpha A_x \sin(\sqrt \alpha x) + \sqrt \alpha B_x \cos(\sqrt \alpha x)$.
From the boundary conditions we deduce that:
$$X'(0) = A_x = 0$$
So $X(x) = B_x$, constant solution.

$\alpha < 0$

So $\Delta < 0$ and $X$ is of the form $X(x) = A_xe^{-\sqrt(-\alpha)x} + B_xe^{\sqrt(-\alpha)x}$.
From the boundary conditions we deduce that $A_x = B_x = 0$. So $X \equiv 0$, which is not interesting.

Final formulation

Finally, we can group the two writings into one, in the following form:
$$X_k(x) = C_x\cos(k\pi x), \forall k \geq 0$$
With :
$$C_x = \begin{cases} B_x & \text{ if } k=0 \cr A_x & \text{ if } k>0 \end{cases}$$

Solving $(**)$

In the same way we solve the equation $(**)$. From the characteristic equation we find:
$$\Delta = -4\beta$$
We then distinguish 3 cases:

Case $\beta > 0$

Thus $\Delta < 0$ and $Y$ is of the form . From the boundary conditions we deduce that:
$$\begin{cases}Y(0) = 0 \Rightarrow A_y = 0\\ Y(1) = 0 \Rightarrow \sin(\sqrt \beta) = 0\end{cases}$$
We deduce that $\beta_k = k^2\pi^2, \forall k>0$ and therefore $Y_k(y) = A_y\sin(k\pi y), \forall k>0$.

Case $\beta = 0$

So $\Delta = 0$ and $Y$ is of the form $Y(y) = A_yy + B_y$.
From the boundary conditions we deduce that:
$$\begin{cases}Y(0) = B = 0\ Y(1) = A = 0\end{cases}$$
So $Y \equiv 0$, which is not interesting.

Formulation Finale

The solution to $(**)$ is therefore written:
$$Y_k(y) = B_y\sin(k\pi y), \forall k > 0$$

Rewriting the problem in the Spectral basis

By solving $(*)$ and $(**)$ we deduce that the eigenvalues and eigenvectors of the Laplacian operator are:
$$\begin{cases} w_{k,l}(x,y) = C\sin(l\pi y)\cos(k\pi x)\\ \lambda _k = \pi^2(l^2 + k^2) \end{cases}, \forall k \geq 0$$
According to the Riesz-Fredholm theorem cited above, there exists a Hilbertian basis (orthonormal by definition) of eigenfunctions of $L^2(\Omega)$ belonging to $H^1_0(\Omega)$ and a sequence increasing positive eigenvalues which tend towards infinity.

Normalization of basis vectors

Before rewriting the problem in the basis of eigenfunctions thus found, let us determine the constant $C$ such that the basis is orthonormal in $L^2(\omega)$.
Let's solve:
$$\begin{matrix} ~ \vert \vert w_{l,k}\vert \vert_{L^2(\Omega)} = 1 & \Leftrightarrow & \int_0^1 \int_0^1 C^2\sin^2(l \pi y) \cos^2(k \pi x) dxdy \cr & \Leftrightarrow & C^2 \underset{\frac{1}{2}}{\underbrace{\int_0^1 \sin^2(l \pi y) dy}} \underset{\frac{1}{2}}{\underbrace{\int_0^1 \cos^2(k \pi x) dx}}\cr & \Leftrightarrow & \begin{cases} C = \sqrt 2 & \text{ if } k=0 \cr C = 2 & \text{ otherwise } \end{cases} \end{matrix}$$

Rewriting the problem in the spectral basis

We can write our solution in the orthonormal basis found:
$$ w(x,y) = \sum_{l,k}\xi_{l,k}(t)w_{l,k}(x,y)~~ ~ (S)$$
According to the problem statement, $f\in L^2(\omega)$, which implies that it can also be expressed in the found basis:
$$f(x,y) = \sum_{l,k}f_{l,k}(t)w_{l,k}(x,y)$$
With $f_{l,k}(t) = <f(x,y,t), w_{l,k}(x,y)>_{L^2(\omega)}$, i.e. the projection of $f$ on the vectors of the Hilbertian basis ${w}_{l,k}$.
We can therefore rewrite the initial problem in $w$ as follows:
$$\begin{matrix} ~ & \frac{\partial w}{\partial t} - \Delta w = f(x,y,t) + 1-y(1+e^{-t}(x^2-x+\frac{3}{2})) & \cr \Leftrightarrow & \frac{\partial }{\partial t}\sum_{l,k}\xi_{l,k}(t)w_{l,k}(x,y) - \Delta \sum_{l,k}\xi_{l,k}(t)w_{l,k}(x,y) & \cr & = \sum_{l,k}f_{l,k}(t)w_{l,k}(x,y) + 1-y(1+e^{-t}(x^2-x+\frac{3}{2})) \end{matrix}$$
As normal, for a completely explicit resolution, it would be necessary to project $1-y(1+e^{-t}(x^2-x+\frac{3}{2})$ on the basis of ${L^2(\omega)}$ found, but this is of little interest.
We will simply put $\overset{\sim }{f}(x,y) = f(x,y) + 1-y(1+e^{-t}(x^2-x+\frac{3}{2})$. This allows us to rewrite the problem as follows:
$$\frac{\partial }{\partial t}\sum_{l,k}\xi_{l,k}(t)w_{l,k}(x,y) - \Delta \sum_{l,k}\xi_{l,k}(t)w_{l,k}(x,y) = \sum_{l,k}\overset{\sim }{f}_{l,k}(t)w_{l,k}(x,y)$$
We use the linearity of the differentiation operators to bring the Laplacian and the time derivative into the sum and we use the relation $\Delta w = -\lambda w$ to factor like this:
$$\sum_{l,k}[\xi'_{l,k}(t) - \lambda_{l,k} \xi_{l,k}(t)]w_{l,k}(x,y) = \sum_{l,k}\overset{\sim }{f}_{l,k}(t)w_{l,k}(x,y)$$
Which amounts to solving an ordinary differential equation:
$$\xi'_{l,k}(t) - \lambda_{l,k} \xi_{l,k}(t) = \overset{\sim }{f}_{l,k}(t)$$

Solving the ordinary differential equation

The resolution of the ODE is done very simply, using the initial conditions to obtain the $\xi_0$ and the explicit expression of the $\overset{\sim }{f}_{l,k}(t) $ to $0$.

The technique that works every time is obviously variation of the constant. I leave it to the reader to make these calculations in order to obtain the $\xi_{l,k}(t)$.

Final solution

Once the EDO is resolved, the job is complete. If an explicit expression is really desired, simply replace the $w_{l,k}$ with the $\xi_{l,k}$, which explicitly depends on the $\overset{\sim }{f}_{l ,k}$ in the expression of $(S)$.

We could possibly work a little more to give, in a more general framework, the link between the $\xi_{l,k}$, the $\overset{\sim }{f}_{l,k}$ and the own values. Likewise, if the source term $f$ is given explicitly, we must calculate its projection on our basis in order to obtain a completely explicit solution.

References

[1]: Analyse Numérique et Optimisation, G. Allaire, Chapter 7.

‌

Introduction

In this article, we are interested in solving Initial Value Problem (IVP) using a special isomorphism between the algebraic fractions in $\delta'$ and the usual algebraic fractions. For that, we will use the space of distributions with compact and positive support and observe that the convolution of two distributions behaves like the multiplication over a field.

Usually, this approach to solve an IVP simply presents Laplace transforms as a tool, without much explanation about where it comes from and why it works. I took a different approach with the point of view of distributions which is the natural space to perform convolution. In my opinion, it leads to naturally observe the similar structure between this particular subspace of distributions and the polynomials. On top of that, the concrete example presented at the end of the article, which is solving a non-homogeneous ordinary differential system with constant coefficients, would not even require the mention of Laplace transform, nor to compute explicitly the transformation, but using only a symbolic table deduced directly from a simple convolution equation and observing the effect of convolving with the derivatives of the Dirac distribution.

This article starts by presenting the space of distributions and in particular the regular distributions, as well as the convolution operator.

Space of Distributions

We denote by $\mathcal D$ the set of test functions. A test function $\phi \in \mathcal D$ is a function from $\phi: \mathbb{R}^n \mapsto \mathbb{R}$. The set of test functions is associated with an extremely strong notion of convergence:

$$ \phi_n \underset{\mathcal D}{\to} \phi  \Leftrightarrow  \begin{cases} \exists[a,b], ~ \forall n, ~~  \phi_n = 0~~\text{outside of } [a,b]\ \phi = 0~~\text{outside of } [a,b] \ \forall k, ~~ \phi_n^{(k)} \underset{u}{\to} \phi^{k} \end{cases} $$

In other words, the tests functions have a compact support and are $C^{+\infty}$.

The space of distribution is the topological dual of $\mathcal D$, denoted $\mathcal{D}^*$, i.e., a distribution $T$ is a continuous linear map on $\mathcal D$. In other words,

$$T \colon \begin{aligned} \mathcal{D} &\to \mathbb{C}\\\\ \phi & \to <T,\phi> \end{aligned}$$

with the following properties:

• $<T, a\phi + b\varphi> = a<T, \phi> + b<T,\varphi>$
• Si $\phi_n \to \phi$ alors $<T, \phi_n> \to <T, \phi>$

We say that $T_n \overset{\mathcal{D}'}{\to} T$ if and only if $\forall \phi \in \mathcal{D}, ~~ <T_n, \phi> \to <T,\phi>$.

Regular Distributions

We denote by  $L^1_{\text{loc}}$  the set of measurable and locally integrable functions, i.e. $\forall a < b \in \mathbb{R},~~ \int_a^b \vert f(x) \vert dx < \infty$.

It is relatively easy to see that if $f \in L^1_{\text{loc}}$ et $\phi \in \mathcal D$ then the product $f\phi$ is integrable. Then, any function $f \in L^1_{\text{loc}}$ uniquely defines a distribution, denoted $f$ by abuse of notation:

$$f \colon \phi \in \mathcal{D} \to <f,\phi> \overset{def}{=} \int_{-\infty}^{+\infty} f(x)\phi(x)dx$$

Therefore, we can identify $L^1_{\text{loc}}$ as a subset of the space of distributions $\mathcal{D}'$. This subset defines regular distributions. A distribution that is not regular is singular.

Convolution, the multiplication of distributions

Given $f$ and $g$, two measurable functions, then the $f$ and $g$ can be convolved if and only if for almost all $x$

$$h(x) = \int \vert f(u)g(x-u)\vert du < +\infty$$

Then, we define convolution between two functions by:

$$f \star g(x) = \int f(u)g(x-u) du, \text{almost everywhere}$$

Convolution obviously has the commutativity property. It has regularization property on $f$. For instance, if we take $g(x) = \frac 1 {2h} 1_{[-h,h]}(x)$, $f \star g$ represents the moving average of $f$.

We can generalize the convolution operator for the space of distributions. We say that $S$ and $T$ are convolvable if and only if $\forall \phi \in \mathcal D$, the expression $<S_x, <T_y, \phi(x+y)>>$ is well-defined.

At first sight, convolution of distributions looks like composition of functions. Obviously, one might wonder when this expression is well-defined. Unfortunately, we would have to introduce the space of Schwartz distributions, which is outside the scope of this article. Therefore, we will only invoke $D_+'$, the space of distributions with positive and compact support, i.e a subset of $[0,+\infty[$.

The good news is that the convolution is well-defined for any function of $D_+'$. In particular, the Dirac distribution is the neutral element for $\star$. Indeed,

$$<\delta \star T, \phi> = <\delta_x, <T_y, \phi(x+y)>> = <T_y, \phi(0+y)> = <T,\phi>$$

Finally, let us consider the effect of the derivative of the Dirac by the convolution operator:

$$<T\star \delta', \phi> = <T_x,<\delta'_y, \phi(x+y)>> = <T_x, -<\delta_y, \phi'(x+y)>> = <T', \phi> $$

That is to say $T\star \delta' = T'$ and $T \star \delta^{(k)} = T^{(k)}$.

The convolution of a distribution $T$ by the $k$-th derivative of the Dirac distribution is the $k$-th derivative of the distribution!

Laplace Transform and the algebra over $D_+'$

$(D_+', +, \times, \star)$ is an associative, commutative and unitary algebra. Precisely, it is a Banach $\star$-algebra usually denoted $L^1(G)$. $\times$ denotes the scalar multiplication. We denote by $T^{\star n}$ the $n$ times composition of the convolution of $T$ by itself. For any locally integrable function $f$, $Hf$ belongs to $D_+'$ where $H$ is the Heaviside distribution^[1].

Notice that $\delta'' = \delta' \star \delta' = {(\delta')}^{\star^2}$ and that a polynomial of $\delta'$ is defined by $P(\delta') = \delta^{(n)} + a_{n-1}\delta^{(n-1)}+\ldots+a_0\delta$. Therefore, if $T\in D_+'$ then $P(\delta')\star T = T^{(n)}+ a_{n-1}T^{(n-1)}+\ldots+a_0T$.

We can create an isomorphism $\mathcal L$ between algebraic fractions in $\delta'$ and the usual algebraic fractions such that:

$$\begin{array}{c  c  c} \star & \leftrightarrow & \times\\ \delta & \leftrightarrow & 1\\ c\delta & \leftrightarrow & c \\ \delta' & \leftrightarrow & p \\ \delta^{(n)} & \leftrightarrow & p^n\\ P(\delta') & \leftrightarrow & P(p) \end{array}$$

In particular, $\mathcal L$ needs to satisfy the following conditions: $\mathcal L{ f \star g} = F(s)G(s)$ and $\mathcal L^{-1}{ f g} = (f \star g)(t)$ where $F = Hf$ and $G = Hg$.

This transformation is known as Laplace transform and can be explicited by:

$$\forall f\in L^1_{\text{loc}}, ~ \mathcal L[f](\lambda) = \int_0^{+\infty} f(t)e^{-\lambda t}dt$$

However, we do not need to compute any transformation analytically, and we can just work with the property of our newly defined algebra.

It is easy to show that $P(\delta')^{\star -1} \leftrightarrow \frac 1 {P(p)}$. Hence, $Q(\delta')\star P(\delta')^{\star -1} \leftrightarrow \frac {Q(p)} {P(p)}$.

Another remarkable relation is given by $He^{\lambda t} = (\delta'-\lambda \delta)^{\star -1} \leftrightarrow \frac 1 {p - \lambda}$

If the product of two distribution is not defined in general, the product of a distribution by a $C^{\infty}$ function is not a problem. ↩︎

Some examples of correspondance

From the previous relation, we can deduce the following one:

$H\frac {t^{n-1}}{(n-1)!} e^{\lambda t} \leftrightarrow \frac 1 {(p - \lambda)^n}$

In particular $\lambda = 0$ : $H\frac {t^{n-1}}{(n-1)!} \leftrightarrow \frac 1 {p^n}$

Obviously, with the exponential decomposition, we can obtain the relations with all trigonometric functions:

• $H \frac {sin (\omega t)}{\omega } \leftrightarrow \frac 1 {p^2 + \omega^2}$
• $H cos (\omega t) \leftrightarrow \frac p {p^2 + \omega^2}$
• $H \frac {sh (\omega t)}{\omega } \leftrightarrow \frac 1 {p^2 - \omega^2}$
• $H ch (\omega t) \leftrightarrow \frac p {p^2 - \omega^2}$

A list of Laplace transform for some usual functions can be found here.

Application to the Initial Value Problem (IVP)

Consider the following system of ordinary differential equations:

$$\left\{\begin{aligned}
x'+y' = f\\
x+y''=g
\end{aligned}\right.$$

with $f$ and $g$ in $L^1_{\text{loc}}$, given, and the initial values $x(0)$, $y(0)$ and $y'(0)$ also given.

This system is non-homogeneous because of $f$ and $g$, coupled because the dependent variables appear in more than one equation, and linear.

To solve this, the first step is to transform our system in $D_+'$ by defining $X = Hx$, $Y=Hy$, $F=Hf$ et $G=Hg$ , définie par $0$ sur $\mathbb{R}^{-}$ et $1$ sur $\mathbb{R}^{+}$.

To completely rewrite the system, we need to obtain $X'$ as well as $Y'$ and $Y''$. By applying the derivative of Heaviside, we obtain $X'=Hx'+x(0)\delta$ and $Y'=Hy'+y(0)\delta$.
Doing the same for $Y'$ leads to $Y''=Hy''+y'(0)\delta + y(0)\delta'$.

Thus, our system is equivalent to the following in $D_+'$:

$$\left\{\begin{aligned}
X'+Y'& = F + (x(0)+y(0))\delta\\
X+Y''& = G + y(0)\delta'+y'(0)\delta
\end{aligned}\right.$$

The second step is to apply Laplace transform to obtain a pure algebraic system:

$$\left\{\begin{aligned}
pX + pY & = F + x(0)+y(0)\\
X+p^2Y & = G + y(0)p+y'(0)
\end{aligned}\right.$$

Which is equivalent to:

$$\begin{aligned}
\begin{pmatrix}
p & p\\
1 & p^2
\end{pmatrix} & \begin{pmatrix}
X \\
Y
\end{pmatrix} & = & \begin{pmatrix}
F + x(0)+y(0) \\
G + y(0)p+y'(0)
\end{pmatrix}
\end{aligned}$$

In other words, our initial system is equivalent to a the linear system $AZ = B$.

The determinant of $A(p)$ is $det(A(p)) = p(p^2-1)$ such that we can calculate the inverse of $A$ and solve this sytem very easily.

$$A(p)^{-1} = \frac 1 {p(p^2-1)} \begin{pmatrix}
p^2 & -p\\
-1 & p
\end{pmatrix}$$

Therefore,

$$\begin{aligned}
\begin{pmatrix}
X \\
Y
\end{pmatrix} & = & \frac 1 {p(p^2-1)} \begin{pmatrix}
p^2 & -p\\
-1 & p
\end{pmatrix} &
\begin{pmatrix}
F + x(0)+y(0) \\
G + y(0)p+y'(0)
\end{pmatrix}
\end{aligned}$$

To obtain $X$ and $Y$, it is enough to apply the inverse of Laplace transform.
For $x$:
$X = (x(0)+y(0))\frac p {p(p^2-1)} - y(0)\frac p {p^2 -1} - y'(0)\frac 1 {p^2 -1} + F \frac p {p^2 -1} - G \frac 1 {p^2 -1}$

For $t>0$, we have $x(t) = x(0)ch(t) - y'(0)sh(t)+\int_0^tf(s)ch(t-s)-g(s)sh(t-s)ds$.

Using a partial fracton decomposition, we obtain $y(t) = -x(0)(ch(t) -1)+y'(0)sh(t)+\int_0^tg(s)sh(t-s)-f(s)(ch(t-s)-1)ds$.

A basic verification consists in taking $t=0$ to observe that the initial conditions hold.

To summarize, the method consists in rewriting the problem in $D_+'$, then expressing the system in a symbolic form, which is implicitely applying Laplace transform. The new problem is basically a linear system that can easily be solved for $p$. Once the solution is obtained, the solution to the original problem can be found thanks to the inverse transform.

Conclusion

In this article, I tried to show another point of view on Laplace transform, and in particular, how it is possible to leverage the convolution algebra to avoid to have to explicitely compute Laplace transforms (or to refer to tables).
The isomorphism between the algebraic fractions in $\delta'$ and the usual algebraic fractions is absolutely remarkable and beautiful because it allows to solve analytically complex non-homogeneous, coupled and linear systems of ordinary differential equation without manipulating a single integral, but with basic algebra only.

Notice that the method holds for systems without constant coefficient, but requires to solve more complex convolution equations to determine the correspondance. In some cases, it might be easier to refer to tables of Laplace transform. But it is obviously less beautiful!

Introduction

Numerical computations are always tainted by errors. A typical example is when scientists simulate a physical system by a numerical method for which a mathematical study gives a bound on the error depending on some parameters. For example, with the Euler method, the error compared to the real solution decreases with the time step. If we did the numerical computation by hand, a time step that goes towards $0$ would allow us to find the exact solution. Obviously, it is impossible, even for a machine, to have a time step of $0$, but one might think that reducing it as much as our time budget or computing power allows is a good thing. Absolutely not!

Indeed, a second category of errors, not connected to the method but much more general, is the representation  error . The problem is far more trivial, almost crude in the simplicity of its statement: it is impossible to represent an infinite quantity in memory space of finite size.

This leads us to consider the fact that, whatever the representation we chose, i.e., the way of translating a real number for the hardware, there will always exist numbers for which the representation will not be possible. An army of engineers and researchers worked to find and formalize a practical and intelligent representation, as a compromise between precision and ease of handling. This standardization process resulted in the IEEE-754 standard followed by most of the world's computer hardware.

A question that naturally arises is: can I numerically control the errors which are induced by the representation error and then propagated during the calculation? This is what we will try to answer positively thanks to the CESTAC method, which stands for STochastic Control and Estimation of Rounding of Calculations (fr.: Contrôle et Estimation STochastique des Arrondis de Calculs).

An implementation of the method described in this article can be found in the following repository:

Some floating point arithmetic reminders

Let us first see how to represent a real number in scientific notation and in a any base

We denote by $b$ the base of the arithmetic in which we will work, with generally $b=2$ or $b = 16$ for modern computation units. Then, any number $x \in \mathbb {R}$ can be written by:

$$x = \pm mb^E$$

With $\frac{1}{b} \leq m < 1$ and $m$ the significand, possibly having an infinite number of digits after the comma, and $e$ an exponant which is an integer expressed in base $b$.

We can rewrite the significand in base $b$ such that $\sum_{i}^n m_ib^{-i}, 0 \leq m_i < b$ with $n \in \mathbb{N} \cup +\infty$.

As we mentioned in the introduction, since a computer has only a finite amount of memory, it is impossible for it to store an infinite amount of information. Worse, whatever the base $b$ chosen, there exists an infinite number of numbers whose represensation include a significand having an infinite number of digits. In other words, it is impossible to perfectly represent the set of reals with a computer. Real numbers are generally approximated by numbers called floating point numbers.

Let us give a succinct proof.

This is how for a machine, a real number $x$ is represented by a floating number $X$ which can itself be written as follows:

$$X = \pm Mb^E$$

With $\frac{1}{b} \leq M < 1$ and $M$ the significand encoded on a finite number of digits $n$ and $E$ the exponent, also encoded on a finite number of digits. We can therefore write $M$ in base $b$ by: $\sum_{i=1}^n M_ib^{-i}, 0 \leq M_i < b$, where this time the number of elements to be summed is always finite.

As these are only reminders, I am not going into all the intricacies of the IEEE-754 standard, and these explanations are sufficient to continue the article.

We consider the assignment operation ($:=$): $\mathbb {R} \to \mathcal {F}$, where $\mathcal{F}$ is the set of floating point numbers. That is to say the operation which associates to a real number its machine representation.

To concretely illustrate the banality of the thing via C++:

double x = 0.1;

For a given real $x$, there exists a float $X^+$ and a float $X^-$ such that $X^- \leq x \leq X^+$ and such that there is no float $Y$ and $Z$ such that $X^- < Y < x < Z < X^+$. In other words, $X ^ +$ and $X ^-$ are the floats immediately greater than and less than $x$.

If $x$ is representable in an exact way, then we have equality between the three terms and the assignment operation $X: = x$ is unambiguously defined.

If this is not the case, as in the above example with $0.1$ in base 2, we must choose a representative between $X^+$ and $X^-$. At first sight, none of them are more legitimate to represent $x$.

This is where the IEEE-754 standard comes in to offer four rounding modes to remove ambiguity about the assignment operation. Here is a brief description:

• Rounding towards $+ \infty$ (or by excess): we return $X ^ +$ except if $x = X ^ -$.
• Rounding to $- \infty$ (or by default): we return $X ^ -$ unless $x = X ^ +$.
• Round to $0$: returns $X ^ -$ for positives and $X ^ +$ for negatives.
• Rounding to nearest: returns the machine number closest to $x$.

An essential property of the IEEE-754 standard is that it guarantees that the result of a floating point operation is equal to the result of the corresponding real operation to which the rounding mode is applied rightafter. In other words, if we choose a rounding mode $\text{Arr}$, $a$ and $b$ two real numbers whose floating point representations are $A$ and $B$, $+$ a real operation, and $\oplus$ its machine counterpart, then, the standard guarantees us that $A \oplus B = Arr(a + b)$.

This property, known as correct rounding, is essential because it makes it possible to prove on a numerical algorithm or to obtain proven bounds for numerical results.

Finally, we can define the relative assignment error by the following formula: $\alpha = \frac {X - x}{X}$. It is precisely this initial error of representation which propagates during the computation.

The different problems caused by numerical computation

s mentioned earlier, any algorithm that performs floating point computations gives an error-ridden result. When an algorithm is finite (where finite is understood as a finite number of steps to find the solution to a problem), then the numerical error is the result of the propagation of rounding or truncation errors during floating point operations.

In the case of iterative algorithms, for example Newton's method, it is also necessary to stop the algorithm after a certain number of iterations, optimal if possible, which is a problem considering that:

• if the algorithm is stopped too early, the solution obtained will not be a good approximation of the real solution. It is the rate of convergence which informs us about this, that is to say the mathematics behind a specific method;
• if the algorithm is stopped too late, additional steps will not bring more precision to the solution obtained, and worse, the propagation of errors can degrade the solution, until, for pathological cases, returning a result which has no meaning.

What is important here is that given an iterative numerical method, the objectives of the mathematician and the engineer are in a way opposed: the former would like to continue to iterate as much as possible (that is, as long as the time constraints allows it) because he knows that this leads to a better solution in theory, while the engineer tells us that we must stop at some point.

In reality, there are at least four interesting and central issues:

1. For the mathematician: given some hardware and a system of representation, how can I obtain a better approximation of the solution to my problem?
   Answer: find algorithms with higher convergence rate or methods to speed up convergence! In this regard, one might cite the method of Aitken's Delta-2 or the $\epsilon$-algorithme.
2. For the engineer: for a given algorithm AND some hardware with a system of representation, how can we get a better solution approximation?
   Answer: reorganize the operations within the algorithm to limit the error propagation while not changing the convergence rate! A generic technique of reorganizing the terms of a sum is called Kahan's summation algorithm.
3. For everyone: for a given algorithm and its implementation, how can we get the best out of it?
   Answers: Choose the most suited representation of real numbers (symbolic system, decimal32, etc. which generally requires better hardware performance) or increase the encoding size of reals (from simple to double precision or quadruple precision, etc. which only consists in increasing the number of bits to encode the significand and the exponent to represent real numbers, which again requires better hardware performance).
4. For the numericist: how to determine the optimal number of iterations to be performed by an algorithm, whatever the input data? How far is my numeric solution from its real equivalent?
   Answer: Find methods to estimate the numerical precision of a result, which involves estimating the propagation of rounding errors!

CESTAC method attempts to solve the last problem and will be presented in the following section. However, as we will see below, we cannot do it without some knowledge about the other problems.

No, yes, and no. Let us take an extremely simple pathological case: $x_n = ax_ {n-1} - b$, which is an extremely simple computation. Here is also a C++ implementation with a particular initialization:

#include 
#include 
#include 

int main() {
     using namespace std;

     double    b=4095.1;
     double    a=b+1;
     double    x=1;

     for (int i = 0; i < 100; ++i) {
          x = (a*x) - b;
          cout << "iter " << i << " - " << setprecision(numeric_limits::max_digits10) << x << '\n';
     }
     return 0;
}

Which gives for output (directly available here):

iter 0 - 1.0000000000004547
iter 1 - 1.0000000018630999
iter 2 - 1.0000076314440776
iter 3 - 1.0312591580864137
iter 4 - 129.04063743775941
iter 5 - 524468.25500880636
iter 6 - 2148270324.2415719
iter 7 - 8799530071030.8047
...
iter 88 - 3.519444240677161e+305
iter 89 - inf

While the expected mathematical result is $1$, constant for each iteration, we observe that after a few iterations on the machine there is a fast divergence towards infinity. After only 4 iterations, the number of significant digits between the exact actual result and its floating bias is 0!

Stochastic control and estimation of rounding

CESTAC core and error propagation

The algebra $\mathcal{F}$ over the field of floating point numbers is not associative nor distributive. In other words, the order in which we perform our arithmetic operations has an impact on the result.

The correct rounding property ensures the commutativity.

From now on, consider $f$, a procedure acting on $\mathbb {R}$, and its image $F$, a procedure acting on $\mathcal{F}$. Because of the non-associativity, the image of $f$ is actually not unique and there are several procedures which mathematically transcribe $f$ exactly. And obviously, these procedures, due to the roundings, will not return the same float.

As in practice, the algorithms are a succession of small computation steps, as illustrated above on the evaluation of a polynomial, the computation will propagate the errors operation after operation. In optimistic scenarios, the errors compensate each other or are too small and the result is remarkably close to what the precision of the representation allows (however, it is impossible to exceed 16 significant digits in double precision, by definition!), but in the worst case scenario, the result can be totally irrelevant.

In summary, from a procedure $f$ over the field of real numbers, there are several procedures $(F_i)_{0 \leq i <K}$ that we can obtain by permuting the elementary operations and, which theoretically offer in the worst case $K$ different results. On top of that, there is a perturbation of the result by the chosen rounding mode which further exacerbates the worst case ^[1].

The main idea behind CESTAC is to take advantage of the great variability of the results that can be obtained by a sequence of numerical computations. For this, we use some perturbations on the result of an operation and some permutations of the operands, in order to estimate the number of significant figures of a numerical result. By propagating the numerical errors in different random ways, we will be able to distinguish the variable part - the part stained with errors, or not representative -, from the fixed part - the exact part -.

Finding the number of significant digits

If we have a procedure $F$ that we run $N$ times with a random perturbation and permutation each time, we get a sample $R = (R_0, R_1, ..., R_ {N-1})$ of results. We can therefore consider $F$ as a random variable with values in $\mathcal{F}$, with an expected value $\mu$ and a standard deviation $\sigma$. The expectated value $\mu$ can be interpreted as the expected result of the algorithm, i.e. the floating point number that encodes our real solution $r$. The error compared to this expectation, that is to say $\alpha = |r - \mu|$ is the loss of precision that one is entitled to expect from performing numerical computations in floating point. The problem is that $\mu$ is not known, and therefore, we have to estimate it.

In this context, with the hypothesis that the elements of $R$ come from a Gaussian distribution (which is verified in practice), the best estimator of $\mu$ is the mean of the sample $R$:
$$\bar R = \frac 1 N \sum^N_i R_i $$

Similarly, the best estimator of $\sigma^2$ is given by:
$$S^2 = \frac 1 {N-1} \sum^N_i (R_i - R)^2 $$

A classic use of the central limit theorem allows us to build a confidence interval for the exact value $r$ for a threshold $p$:

$$\mathbb{P},(r \in [\bar R - t_{N-1}(p) \frac{S}{\sqrt{N}}; \bar R + t_{N-1}(p) \frac{S}{\sqrt{N}}]) = 1 - p$$

Where $t_{N-1} (p)$ is the value of the cumulative distribution function of Student for $(N-1)$ degrees of freedom and a threshold of $p$.

From this interval, it is possible to calculate the number of significant digits $C$ of our estimator $\bar R$:

$$C_{\bar R} = \log_{10}(\frac {|\bar R|} {S}) - K(N, p)$$

where $K$ depends only on $N$ and $p$, and such it that tends towards $0$ with $N$ increasing. The value of $p$ is fixed in practice at $0.05$, which makes it possible to obtain a confidence interval of $95%$. Here is now the value of $K$ obtained as a function of $N$, for $p = 0.05$:

This may seem surprising but using a sample of size $N = 3$ results in $K$ being less than $1$, i.e., on average, there is no significant digit loss for the sample $R$. In fact, increasing this number is useless because the length of the interval evolving in $\frac 1 N$, to obtain an additional significant figure, it would be necessary to multiply $N$ by 100 (because of the $\log_{10}$)!

Constructing the sample $R$

Now that we have the theory, we need to know how to construct a sample of results $R$ that is as representative as possible of the multitude of results obtainable from our procedure $F$.

For that, we have a perturbation function, pert, which for a particular float $X$, returns a disturbed float $X'$ such that $X'$ is $X$ for which we modified the last bit of its significand in a random and uniform way. In other words, we add to the last bit of significand $-1$, $0$ or $1$ with a probability of $\frac 1 3$.

This perturbation consists in choosing randomly among $X^+$ and $X^-$, which we mentioned in the first part, in order to simulate the propagation of rounding errors.

We use pert whenever an assignment is made, whether it is an initial assignment as a floating variable declaration, or the result of multiple computations.

We also have a permutation operator, perm, which for each assignment operator will randomly modify the term on the right by performing one of the permutations authorized by associativity and distributivity. In other words, it is a question of choosing between $F_1$, $F_2$ or $F_3$ in the example above.

Thus, we will not consider permutations in the following.

From there, there are two ways to use pert to create a sample $R$ and estimate the number of significant digits of a numeric result.

Asynchronous version

The asynchronous version consists in performing our perturbations at each assignment, and building our sample $R$ as the result of $N$ calls to the procedure $F$. In other words, the $N$ calls are independent, hence the name asynchronous. Once the sample is at our disposal, we compute the number of significant digits a posteriori.

Illustration of the asynchronous method with an iterative sequence defined by $X_n = F(X_ {n-1})$ and for initial term $X_0$ with $N = 3$:

$$
\begin{matrix}
& \nearrow X^0_1 = \text{pert}(F(X_0)) & \to X^0_2 = \text{pert}(F(X^0_1)) & \to \dots \to & X^0_n = \text{pert}(F(X^0_{n-1})) & \searrow & \\\
X_0 & \to X^1_1 = \text{pert}(F(X_0)) &\to X^1_2 = \text{pert}(F(X^1_1)) & \to \dots \to & X^1_n = \text{pert}(F(X^1_{n-1})) & \to & C((X_n^i)) \\\
& \searrow X^2_1 = \text{pert}(F(X_0)) & \to X^2_2 = \text{pert}(F(X^2_1)) & \to \dots \to & X^2_n = \text{pert}(F(X^2_{n-1})) & \nearrow & \\\
\end{matrix}
$$

Apparently logical, this method comes up against two major problems.

• As the propagation of the errors is not necessarily being done in the same way, it is possible that the executions of the procedure converge towards different real numbers, in which case the sample is inconsistent. This may be the case if the problem to be solved admits of several solutions for example.
• From one execution to another, as there are certainly conditional branches, two results can come from a series of different branches because of rounding errors. Therefore it is not relevant to compare these results.

For these two reasons, the asynchronous version is generally inapplicable.

Synchronous version

Conversely, the synchronous version consists in modifying the sample at each assignment operation and using the empirical average as value for the conditional branches. In practice, we can also systematically use $X_i$ for a given $i$, which avoids having to calculate the average each time. It is possible to give an estimate of the number of significant digits at any time because the sample is available at all times during the execution. In fact, this answers the two problems of the asynchronous version:

• At each step, the result is consistent with itself, it cannot be different solutions since there is never only one value which is used for the conditional structures.
• The series of branches will necessarily be unique by construction, which makes the final result consistent.

Illustration of the synchronous method on the same example as before:

$$
\begin{matrix}
& \nearrow X^0_1 = \text{pert}(F(X_0)) & \searrow & & \nearrow X^0_2 = \text{pert}(F({X^1_1})) \searrow & & \nearrow & X^0_n = \text{pert}(F(X^1_{n-1})) & \searrow & \\\
X_0 & \to X^1_1 = \text{pert}(F(X_0))& \to & \bar{X_1} &\to X^1_2 = \text{pert}(F(X^1_1)) \to & \dots &\to & X^1_n = \text{pert}(F(X^1_{n-1})) & \to & C((X^i_n))\\\
& \searrow X^2_1 = \text{pert}(F(X_0))& \nearrow & & \searrow X^2_2 = \text{pert}(F(X^2_1)) \nearrow & & \searrow & X^2_n = \text{pert}(F(X^2_{n-1})) & \nearrow &\\\
\end{matrix}
$$

Notice the synchronization points after each step, hence the method name.

CESTAC on an iterative algorithm

Now that we know how to determine the number of significant digits of a numerical result, we will focus on the optimal stopping problem. The exact solution to our problem is noted $x^\star$. We have an iterative algorithm, which gives at iteration $k$, the approximate solution $x_k$, and we know that this algorithm converges after a number of steps potentially infinite, that is to say that $x_k \to x^\star$.

Finally, we have an implementation of our algorithm which at each iteration provides an approximate solution $X_k$ tainted by numerical errors.

A classic stop criterion for iterative algorithms at a given precision is the test $|| x_k - x_ {k-1}|| <\epsilon$, where $\epsilon$ controls the precision. A variant is given by $|| x_k - x_ {k-1} || <|| x_k || \epsilon$. This test is extremely robust in usual arithmetic having infinite precision since it makes it possible to detect when an iteration no longer brings any significant gain in precision. Conversely, in floating point arithmetic, since all $X_k$ are tainted with errors, the subtraction of very close terms leads to small values ​​which may not be significant at all.

The worst possible scenario is the following: $\epsilon$ is chosen too small, and the computational errors that have accumulated are too large compared to $\epsilon$, so the algorithm will never converge, its solution will degrade to either converge towards an inconsistent value, or to diverge outright!

As CESTAC purpose is to determine the number of significant digits of a result, we can therefore use it to find the machine zeros and modify our stop test accordingly, which becomes the following one at the iteration $k$:

1. If $C(X_k) = 0$ and $X_k \neq 0$, then the result is tainted with an error greater than its own value and there is no point in continuing: we stop the algorithm.
2. If $|| X_n - X_{n-1} || = \bar 0$, we stop the algorithm since the difference between two iterations only represents computation errors.
3. If the number of iterations exceeds a certain limit $K$, the sequence is considered as non-convergent and the algorithm is stopped.

Going further: discussions on the validity of CESTAC

There are several hypotheses which have been formulated in order to arrive at the formula for the number of significant digits and which lead to the following question: can we reasonably use a Student's test in order to obtain the confidence interval? This is equivalent to wonder whether the estimator $\bar X$ used is biased or not.

The estimator bias

Without making a proper demonstration (we refer the reader to the studies by the creators of the CESTAC method), we give the broad outlines justifying that the mean estimator is unbiased.

he bias of the estimator $\bar X$ is the quantity $E[\bar X] - x$. Assuming that the $(\alpha_i)_i$ follow a uniform distribution over the “proper” interval (that depends on whether we are in truncation or rounding mode), it is enough to correctly choose the $(h_i)_i$ to re-center the $(\alpha_i)_i$ and thus obtain the following result, neglecting higher-order terms:

$$X' = x + \sum_{i = 1}^n d_i 2^{-t}(z_i)$$

Where the $z_i$ are identically distributed and centered variables, so that $E(X') = x$, i.e. the estimator is unbiased.

The hypothesis of the distribution of $\alpha_i$ is validated when there are enough operations in the procedure $F$, in other words that $n$ is large enough. In fact, very precisely, the bias is not zero but we can show that it is of the order of a few $\sigma$ which skews the final estimate by less than one significant figure.

Validity of the Student test

As we have seen, the hypothesis about the distribution of $\alpha_i$ is satisfied in theory and in practice. But on the other hand, to conclude that the estimator is unbiased, we made an additional assumption: the higher order terms are negligible.

However, while it is easy to see that addition or subtraction does not create an error of second order, this is not the case for multiplication or division, since by considering $X_1 = x_1 + \epsilon_1$ and $X_2 = x_2 + \epsilon_2$, these operators are written:

$$X_1X_2 = x_1x_2 + x_1 \epsilon_2 + x_2 \epsilon_1 + \epsilon_1 \epsilon_2 $$

$$\frac{X_1}{X_2} = \frac{x_1}{x_2} + \frac{\epsilon_2}{x_1} - \frac{x_1 \epsilon_2}{x_2^2} + O(\frac{\epsilon_2}{x_2})$$

The problem for multiplication is that if the respective errors $\epsilon$ for the two operands are preponderant over the exact values ​​$x$, then the second order term becomes preponderant. For the division, the higher order terms become preponderant if $\epsilon_2$ is preponderant w.r.t. $x_2$.

Consequently, there are two complementary ways to ensure the hypothesis behind CESTAC are valid:

• Increase the precision of the representation, i.e., encode each real on more bits, which will reduce the $\epsilon$. In other words, the answer to problem 3.
• Limit the propagation of calculation errors as much as possible, i.e., apply the recipes in response to problem 2.

One possible systematic approach is to include dynamic control over the outcome of multiplication or division operations, at a significant cost.

Conclusion: CESTAC, but for what and for who?

We have seen the different problems that appear when we perform floating point computations and we have given a robust technique to control the propagation of errors induced by these computations. After having presented the foundations of the method and the modalities of use, we applied CESTAC to the optimal stopping problem of an iterative algorithm. Finally, a last part was devoted to the discussion on the validity of the method and a sketch of proof.

The only question that has not been addressed, and which will serve as a conclusion, is: when to use CESTAC? It is obvious that the method presents a significant cost in terms of computation time. Therefore, it must be justified by a gain at least as important. The need for evaluation and error control is crucial, for example, for critical systems such as airplanes. For this reason, the method is widely used in the field of aeronautics, to control both simulations but also directly computations within on-board systems.

Bibliography

• Ingénierie du contrôle de la précision des calculs sur ordinateur, Michèle Pichat et Jean Vignes.
• Validité du logiciel numérique, Jean-Marie Chesnaux, support de cours dispensé à l'UPMC.
• Approche stochastique de l'analyse des erreurs d'arrondi : méthode CESTAC, Jean Vignes et René Alt.
• Handbook of Granular Computing, Witold Pedrycz, Andrzej Skowron et Vladik Kreinovich.

I stress that this is the worst case scenario, which is far from being the common practical scenario. However, standard deterministic studies reason mainly about the worst case scenario, which leads to an overestimation of the computation errors, sometimes rendering these methods obsolete. On the contrary, CESTAC makes it possible not to fall into this trap due to the very construction of the method. ↩︎

If you have been reading my article on the Tuxedo InfinityBook Pro 16 - Gen 8 then you know that I received a laptop with no SSD and therefore, no OS pre-installed. With 2x TB Crucial P3 Plus PCIe 4.0 my goal is to install Ubuntu 22.04 LTS on LVM on a RAID 1 array. However, as you can imagine, this is not as easy as simply using the graphical installer and that is the reason why I write this guide.

A first look at the Tuxedo InfinityBook Pro 16 - Gen8

I was looking at a replacement for my MSI Stealth GS65 that served me well the last 5 years, when I heard, via The Linux Experiment, of Tuxedo, a German company based in Augsburg. They focuse on Linux laptops and in particular, they offer laptops with the Intel Management Engine

EndomorphismAlexandre Quemy

I will assume that you already have a bootable Ubuntu Live with the latest LTS version. Once booted, you should be left with a choice Install Ubuntu or Try Ubuntu

• Select Try Ubuntu.
• Start a terminal by pressing Ctrl-Alt-T.

First, you should identify the name of your disk using lsblk:

AME                    MAJ:MIN RM   SIZE RO TYPE  MOUNTPOINTS
...
nvme1n1                 259:0    0   3.6T  0 disk                       
nvme0n1                 259:3    0   3.6T  0 disk  
...

In my case, the disks are nvme0n1 and nvme1n1 . You should adjust the commands according to the name of your disks.

Create partitions on the physical disks

Remove any previous partition schema using the following destructive operation:

sudo sgdisk --zap-all /dev/nvme0n1
sudo sgdisk --zap-all /dev/nvme1n1

Create two partitions on each drive; one for EFI (type ef00) and one for the RAID device (type fd00).

sudo sgdisk -n 1:0:+512M -t 1:ef00 -c 1:"EFI System" /dev/nvme0n1
sudo sgdisk -n 2:0:0 -t 2:fd00 -c 2:"Linux RAID" /dev/nvme0n1
sudo sgdisk -n 1:0:+512M -t 1:ef00 -c 1:"EFI System" /dev/nvme1n1
sudo sgdisk -n 2:0:0 -t 2:fd00 -c 2:"Linux RAID" /dev/nvme1n1

Create a FAT32 system for the EFI partition on the first drive.

Create the RAID device

Ensure that you have internet connection, either wired, or by configuring a wireless connection which can be done using the graphical interface.

Then install mdadm :

  sudo apt update -y
  sudo apt install mdadm -y

Create the RAID array using the second partition of each disk:

 sudo mdadm --create /dev/md0 --bitmap=internal --level=1 --raid-disks=2 /dev/nvme0n1p2 /dev/nvme1n1p2

You can safely ignore the warning about the metadata.

Partition the RAID device

We now have a RAID device created ( /dev/md0 ) that we can partition using sgdisk . Create a single partition with type LVM (code E6D6D379-F507-44C2-A23C-238F2A3DF928):

  sudo sgdisk --zap-all /dev/md0
  sudo sgdisk -n 1:0:0 -t 1:E6D6D379-F507-44C2-A23C-238F2A3DF928 -c 1:"Linux LVM" /dev/md0

Create the LVM devices

Create a physical volume on the first and unique partition of the RAID device:

sudo pvcreate /dev/md0p1

Create a volume group on the physical volume:

sudo vgcreate vg0 /dev/md0p1

Create logical volumes (partitions) on the new volume group. For my purpose, I will create an extremely simple layout with a single volume that takes the whole RAID device partition that itself takes the whole disk using 100%FREE.

sudo lvcreate -Z y -L 100%FREE --name root vg0

If you want multiple logical volumes, for instance, to separate the root from the home, you will have to use something like this:

sudo lvcreate -Z y -L 100GB --name root vg0
sudo lvcreate -Z y -L 100%FREE --name home vg0

Install Ubuntu 22.04

Double-click on the Install Ubuntu 22.04 LTS icon on the desktop.

Follow the instruction according to your needs until the storage configuration step.

On the Installation type page, select Something else. This will present you with a list of partitions called /dev/mapper/vg0-root, etc.

Double-click on each partition starting with /dev/mapper/vg0-. Select Use as: Ext4, check the Format the partition box, and choose the appropriate mount point (/ for vg0-root, /home for vg0-home, etc.).

Select the first disk /dev/nvme0n1 for the boot loader.

Press Install Now and continue the installation.

When the installation is finished, select Continue Testing.

If every is correct, by typing lsblk in a terminal, you should see something like:

NAME               MAJ:MIN RM   SIZE RO TYPE  MOUNTPOINT
...
nvme1n1                 259:0    0   3.6T  0 disk  
├─nvme1n1p1             259:1    0   512M  0 part  
└─nvme1n1p2             259:2    0   3.6T  0 part  
  └─md0                   9:0    0   3.6T  0 raid1 
    └─md0p1             259:6    0   3.6T  0 part  
      └─vg0-workstation 253:0    0   3.6T  0 lvm   /target
                                                   
nvme0n1                 259:3    0   3.6T  0 disk  
├─nvme0n1p1             259:4    0   512M  0 part
└─nvme0n1p2             259:5    0   3.6T  0 part  
  └─md0                   9:0    0   3.6T  0 raid1 
    └─md0p1             259:6    0   3.6T  0 part  
      └─vg0-workstation 253:0    0   3.6T  0 lvm   /target

Installing mdadm on the system

To install mdadm , we need to chroot into the system. To do so, we need to mount the partitions and bind the devices.

The root partition is mount on /target. If you have more than  one logical volume, you should mount them. For instance, assuming a separate partition for home:

sudo mount /dev/mapper/vg0-home /target/home

Then you can bind some devices

cd /target
sudo mount --bind /dev dev 
sudo mount --bind /proc proc
sudo mount --bind /sys sys

Finally, chroot into the system:

sudo chroot .

Help the system to resolve the repositories' DNS:

echo "nameserver 1.1.1.1" >> /etc/resolv.conf 

We can now install mdadm:

apt install mdadm -y

Configuring the RAID on the target system

We need to check if the array has properly been discovered. Check the file /etc/mdadm/mdadm.conf. It should contain a line near the end similar to

ARRAY /dev/md0 metadata=1.2 name=ubuntu:0 UUID=e1a1650d:138d43bd:d53e1e18:29aa74ea

If the line is present, then simply delete name=ubuntu:0.
If the line is not present, then you can generate it using the following:

sudo mdadm -Db /dev/md0

Check that the output looks like the line above and if so, output it to the end of /etc/mdadm/mdadm.conf:

sudo mdadm -Db /dev/md0 >> /etc/mdadm/mdadm.conf

Then, we need to edit the file and remove the name=ubuntu:0.
Update the module list the kernel should load at boot.

echo raid1 >> /etc/modules

Update the boot ramdisk:

update-initramfs -u

Finally, exit from chroot:

exit

And that's it! We now have a fresh new installation of Ubuntu with RAID 1 and LVM.
While the RAID array is protecting the main partition from a disk failure, the EFI boot partition is not protected via RAID as it resides only on the first partition of the first disk.

Instead of relying on RAID, we will clone the partition and insert clone into the boot chain. Be certain to copy the first partition of the first disk to the first partition of the secon disk:

sudo dd if=/dev/nvme0n1p1 of=/dev/nvme1n1p1 bs=4096

You can check the result with

sudo blkid /dev/nvme[01]n1p1

which should return something like

/dev/nvme0n1p1: UUID="F41E-A163" BLOCK_SIZE="512" TYPE="vfat" PARTLABEL="EFI System" PARTUUID="a10fa1d6-0079-4370-9152-1761b2282c07"
/dev/nvme1n1p1: UUID="F41E-A163" BLOCK_SIZE="512" TYPE="vfat" PARTLABEL="EFI System" PARTUUID="ac3839c7-cf53-47c4-a376-905646e0a193"

Then, to insert the second EFI partition into the boot chain, identify the location of the EFI file with

sudo efibootmgr -v

The path is located after the File in the output:

oot0000* ubuntu        HD(1,GPT,a10fa1d6-0079-4370-9152-1761b2282c07,0x800,0x100000)/File(\EFI\ubuntu\shimx64.efi)

Insert the EFI file into the first partition of the second disk. Note that the -d (for disk) expect a disk and not directly a partition. The partition is specified using the flag -p:

sudo efibootmgr -c -d /dev/nvme0n1 -p 1 -L "ubuntu2" -l '\EFI\ubuntu\shimx64.efi'

In case the first disk fails, the data it protected by the RAID array, and the system will still boot from the EFI partition of the healthy disk. The RAID array would be in degraded mode.

In my case, after a reboot, here is the output of lsblk:

aquemy@ws:~$ lsblk
NAME                    MAJ:MIN RM   SIZE RO TYPE  MOUNTPOINTS
nvme0n1                 259:0    0   3.6T  0 disk  
├─nvme0n1p1             259:2    0   512M  0 part  
└─nvme0n1p2             259:3    0   3.6T  0 part  
  └─md0                   9:0    0   3.6T  0 raid1 
    └─md0p1             259:6    0   3.6T  0 part  
      └─vg0-workstation 253:0    0   3.6T  0 lvm   /
nvme1n1                 259:1    0   3.6T  0 disk  
├─nvme1n1p1             259:4    0   512M  0 part  /boot/efi
└─nvme1n1p2             259:5    0   3.6T  0 part  
  └─md0                   9:0    0   3.6T  0 raid1 
    └─md0p1             259:6    0   3.6T  0 part  
      └─vg0-workstation 253:0    0   3.6T  0 lvm  /

The status of the RAID device:

aquemy@ws:~$ sudo mdadm --detail /dev/md0
/dev/md0:
           Version : 1.2
     Creation Time : Wed Sep 13 17:10:02 2023
        Raid Level : raid1
        Array Size : 3906361152 (3.64 TiB 4.00 TB)
     Used Dev Size : 3906361152 (3.64 TiB 4.00 TB)
      Raid Devices : 2
     Total Devices : 2
       Persistence : Superblock is persistent

     Intent Bitmap : Internal

       Update Time : Sat Sep 16 16:45:40 2023
             State : active 
    Active Devices : 2
   Working Devices : 2
    Failed Devices : 0
     Spare Devices : 0

Consistency Policy : bitmap

              Name : ubuntu:0
              UUID : e1a1650d:138d43bd:d53e1e18:29aa74ea
            Events : 31227

    Number   Major   Minor   RaidDevice State
       0     259        3        0      active sync   /dev/nvme0n1p2
       1     259        5        1      active sync   /dev/nvme1n1p2

The status of the LVM physical volume:

aquemy@ws:~$ sudo pvdisplay 
  --- Physical volume ---
  PV Name               /dev/md0p1
  VG Name               vg0
  PV Size               <3.64 TiB / not usable <4.80 MiB
  Allocatable           yes (but full)
  PE Size               4.00 MiB
  Total PE              953700
  Free PE               0
  Allocated PE          953700
  PV UUID               8ocB3w-R3Dp-76Yr-PE2H-mfC2-D4Wf-cmAG8S

And finally, the status of the logical volume:

aquemy@ws:~$ sudo lvdisplay 
  --- Logical volume ---
  LV Path                /dev/vg0/workstation
  LV Name                workstation
  VG Name                vg0
  LV UUID                kq9vhC-OqZA-QaTv-JkBe-OU1i-PJsA-Sukflj
  LV Write Access        read/write
  LV Creation host, time ubuntu, 2023-09-13 17:22:46 +0200
  LV Status              available
  # open                 1
  LV Size                <3.64 TiB
  Current LE             953700
  Segments               1
  Allocation             inherit
  Read ahead sectors     auto
  - currently set to     256
  Block device           253:0

I was looking at a replacement for my MSI Stealth GS65 that served me well the last 5 years, when I heard, via The Linux Experiment, of Tuxedo, a German company based in Augsburg. They focuse on Linux laptops and in particular, they offer laptops with the Intel Management Engine deactivated! If you are as paranoid as I am, this might sound very tempting. And tempting it was since I finally decided to order the new version of their flagship, namely the InfinityBook Pro 16 - Gen8.

InfinityBook Pro 16 - Gen8 - TUXEDO InfinityBook Pro 16 - Gen8 - TUXEDO Computers

TUXEDO InfinityBook Pro 16 - Gen8 - 16-inch 16:10 Omnia display + Intel Core i7-13700H + optional: NVIDIA GeForce RTX 4070 + 80 Wh battery + magnesium chassis + max 64GB RAM + Thunderbolt 4

TUXEDO ComputersTUXEDO Computers, Germany

They offer up to a GeForce RTX 4070 in an magnesium casing for a total weight of 1.6kg. The elegant design reminds me of a crossover between a MacBook and a Dell XPS. It is possible to personalize the keyboard layout or the logo on the machine. Beyond the hardware, they offer to pre-configure the OS and eventually a virtual machine. They offer full encryption and propose several distributions, including their own, based on Ubuntu LTS called TuxedoOS.

As I still need to have a discrete NVidia GPU for work, I went for a GeForce RTX 4060. However, the goal with this laptop is to have enough autonomy to travel with it in a non-stationary setup, yet enough power to work properly. None of the other model would achieve this compromise as far as I understand.

Here is the full configuration that I ordered:

• NVIDIA GeForce RTX 4060 8GB | Silver casing
• Intel Core i7-13700H
• 64 GB (2x 32GB) DDR5 4800MHz Samsung
• without M.2 SSD
• without Wifi & Bluetooth

The reason why I chose without any SSD is because I was able to grab 2x TB Crucial P3 Plus PCIe 4.0 with a significant discount. However, it has to be noted that the price indicated on their website correspond roughly to retailer price which is a good point for the consumer.

I also decided to order without any Wifi & Bluetooth module because they do not offer Wifi 6E. Therefore, I preferred to ordered such Wifi card separately.

The downside is that I will not be able to provide a feedback on the services related to the OS installation but I will definitely be providing a feedback on the build and upgrade since I will have to open right way!

First, the service: I placed the order on a Sunday, on Monday the order was approved and on Friday the laptop was sent. Meanwhile, the order status was updated throughout the different stages of the production. No need to say that for a custom laptop build on order this is extremely fast! Note however that as I did not requested the installation of any OS, this might speed-up the order.

There is nothing special to say about the packaging: it is robust and does the job. Delivered without a scratch.

The chassis and connectivity

There is usually the Tuxedo logo on the top. For mine, I decided to go without any logo. The chassis feels sturdy. The weight is balanced and it is possible to open the laptop with one finger which is a good point of accessibility.

The only noticeable flex on the chassis is on the edge, between the screen and the bottom part. There is roughly half a millimeter such that a pressure on the screen part, right in the middle, will see it bent toward the bottom part. Nothing to worry about, but it is slightly more pronounced than on my MSI for instance.

The left side is equipped with

• a Kensington Slot,
• a USB-C 3.2 Gen2 that is connected to the NVidia and supports DisplayPort 1.4a and is compatible with G-SYNC. However, it does not support PowerDelivery DC-in such that it is not possible charge the laptop via this port. It supports DC-out up to 15W which is not really fast for modern smartphones, so it is mostly interesting for USB devices that can be powered via USB.
• a USB-A 3.2 Gen1,
• a SD card reader,
• a 2-in-1 audio jack (headphone + mic)

The right side is equipped with

• a Thunderbolt 4 that is connected to the integrated GPU and supports DisplayPort 1.4a and both way PowerDelivery: DC-in up to 100W as per TB4 norm and DC-out up to 15W.
• a USB-A 3.2 Gen1
• a HDMI 2.0b that is connected to the integrated GPU and supports HDCP 2.2.
• a Power plug (DC-In)

The fact that the TB4 and HDMI ports are connected to the iGPU instead of the dGPU is for me absolutely optimal and how it should be on any laptop. Indeed, my laptop is connected to a docking station with two additional screens over Thunderbolt. If the TB4 port was connected to the discrete GPU, it means that I would not be able to disable the discrete GPU while I am docked. However, I don't use the discrete GPU most of the time, only if I want to play or have to run an experiment to run on GPU. Therefore, having the discrete GPU running all the time would be a complete waste of energy, on top of inducing unnecessary heat. On the contrary, with the current wiring, even if I am docked and connected only to the integrated GPU, I am able to switch on the NVidia whenever I need it. Perfect!

The bottom part has a large mesh surface for the airflow. Although it is not visible during the usage, I find the screws particularly visible and ugly. Or maybe ugly because very visible?

Keyboard and trackpad

The keyboard is a full sized keyboard with a numpad! I think it is the first time I see a numpad on a laptop. In my case, having a numpad has no utility and makes it looks crowded. It also forces to adjust my right hand position. On the other hand, having the numpad adds some space between del/ins key and the power button which might be a blessing if, like me, you are prone to fat fingers.

Regarding the keyboard itself, it feels amazing and I immediately adjusted to the keycaps and activation force. I am no expert on keyboards so my words have little interest here, but it is rare that I can adapt that fast to a new keyboard. If I had to compain, I would say that the third function on some of the keys (e.g ~ or ! or @) is printed so small that at first I thought it was a defect on the keycap itself.

The trackpad is really large and reactive. It feels good although, as for the numpad, I rarely use it and don't have a serious point of comparison.

The biggest asset of this part of the chassis in my opinion is the extremely comfortable texture of the material around the keyboard. It is not slippery and soft as glass but rather a sort of velvet feeling that gives the proper grip while typing and prevent finger prints and sweat. Truly amazing.

Opening and upgrading

Opening the laptop is extremely easy: removing the 9 screws is enough for the lid to come easily without force. One small annoyance is the fact that there are two different screw lengths and as a result, one needs to remember carefully where are the longer ones to avoid any bad surprise when screwing them back. When screws cannot be the same length, the screws and corresponding threaded holes should have a color code somehow.

After removing the lid, we have a direct access to the motherboard and all components that are replaceable. This is a good thing, in particular in comparison to my MSI Stealth GS65 on which the motherboard is flipped! The build is extremely tight and clean. Note that the two connectors for the antenna have their endings protected in a small plastic tube and are then taped together.

The two M.2 slots are positioned on the opposite sides. The left one is positioned vertically and crosses the horizontal plan of the battery. If it was possible to find another position on the motherboard, it would probably free enough space to increase the battery capacity from  80Wh to something closer to the 99Wh limit allowed in an airplane.

Weight and mobility

Once everything has been installed, the Tuxedo InfinityBook Pro 16 - Gen8 weights, as promised, 1.6kg and 2.1kg with the entire stock power supply. Note that the power supply is rated 180W and is compact given the hardware in this laptop.

To improve the portability even more, I bought a Anker 737 GaNPrime 120W with 140W cable. It is quite pricey, but on the other hand, the weight of the charger with the cable is half of the stock one, i.e. 150g for a total of 1.85kg. Because a single port can output only 100W, keep in mind that it is not possible to fully use the NVidia GPU or at least, Tuxedo advises to disable it.

A word on the UEFI and software

The laptop shipped with a custom BIOS/UEFI with few interesting features:

1. The possibility to deactivate Intel ME,
2. Wake-on-Lan support
3. The possibility to select a performance versus silence profile
4. SecureBoot support

I am not able to give any feedback on TuxedoOS, a distro based on Ubuntu LTS as I installed myself a Ubuntu LTS. However, I've been able to install Tuxedo Control Center using their repository and the .deb package. TCC is a neat tool to monitor the hardware and manage the power profile, among other. Even on a vanilla Ubuntu, all the hardware components were properly recognized and manageable with the TCC. I would definitely recommend anyone using a Tuxedo laptop to use the Tuxedo Control Center.

More on the Tuxedo Control Center here:

TUXEDO Control Center - TUXEDO Computers

TUXEDO Control Center: [txt=product-info-scrolleffects] Your notebook. Your control. Ihr Browser kann dieses Video nicht wiedergeben. TUXEDO Control Center Dashboard All data at a glance Performance profiles for all use cases Create profiles freely…

TUXEDO ComputersTUXEDO Computers, Germany

If you are interested in installing Ubuntu with a RAID 1 and LVM, you might want to have a look at this article:

Installing Ubuntu 22.04 LTS with RAID 1 and LVM

If you have been reading my article on the Tuxedo InfinityBook Pro 16 - Gen 8 then you know that I received a laptop with no SSD and therefore, no OS pre-installed. With 2x TB Crucial P3 Plus PCIe 4.0 my goal is to install Ubuntu 22.04 LTS

EndomorphismAlexandre Quemy

Conclusion

The Tuxedo InfinityBook Pro 16 - Gen8 is probably the best linux laptop available on the market at the moment I write these lines if you are looking for a good compromise between portability and performances. It has an impressive hardware and the build quality is really high. The customization options are also welcomed if it is something that you are looking for.

From my point of view, the price point is correct, especially given the size of the company and the quality of the laptop.

My only complaint would be the battery. I would have prefered a larger battery as the norm allows for 99Wh in airplanes, but I understand that it is rather hard to fit everything is such a small chassis.

For my personal usage, I've been looking for a main homelab server with a low power consumption (almost) without any compromise. By main, one should understand where most of the services run, as opposed to a storage dedicated server.

The build presented in this article is in production for around 6 months now and has been drawing an average 33W in a 24/7 setup, which, given the specs and the usage, is pretty satisfying. Note that this figure is obtained without any particular tweak (yet) to lower the power consumption and therefore, I am sure it is possible to reach a sub 30W easily, and probably near 20W average.

Requirements

First, here are my requirements:

1. At least one SFP+ for a 10GbE connection,
2. A rackmount case, preferably 2U,
3. As silent as possible because it will be located in my office,
4. As low power as possible.

The difficulty with building a low power server is that the information about the idle power consumption or the PSU's power efficiency are rarely available. In addition, the effective consumption of each component may vary depending on the build, i.e. the other components.

To select the components, I used the following database with low power builds. The power draw of these builds are measured at the wall and therefore gives a good idea of the entire system draw rather than individual components.

-Die sparsamsten Systeme (<30W Idle)-

All Idle Watt,Last Watt,C-State,CPU,Mainboard,Netzteil,GPU,Messgerät,HDD / SSD,RAM,Betriebssystem,Maus, Tastatur, Monitor,Lüfter,User1,5,63,C8,Intel Core i5-9500T,HP 859C​,HP 90W PPP012C-S,iGPU,ina219​,Kingston A2000 500GB NVME​,1 x 8GB F4-2666C19-8GRS DDR4​,Fedora 33​,Dell Maus, HP Keyboard, D…

Google Docs

Components

The CPU I selected is the Intel Core i3-10105 because the benchmark for its non-graphic version provides an average power consumption of 8.60W which represents a good comprise performance/power. The database above even reports a build at 5.4W in C8 state!

CPU avg. power consumption [W] in idle – lower is better

The motherboard I selected is the ASRock H510M-HVS R2.0. The main reason is because it has 2x PCIe 3.0 slots (one full size x16, and one x1) and 4 SATA ports. The PCIe 3.0 x16 will be used for the network card while the additional x1 could be use for future upgrade (e.g. cold storage with a SATA controller card).

In terms of storage, the purpose of this build is to be a service node rather than a storage node. For this reason, I decided to use 3 SSD with a relatively small capacity: one for for the installation of the hypervisor, and two additional in RAID 1 for the containers and VMs. I had a spare (very old) 128GB SSD that I used as boot device and two Western Digital Red SA500 500Gb. Without this spare SSD, I would probably go for a third WD SA500 500Gb for practicity. Of course, if you don't plan to use another node for storage or don't have a NAS, you might want to take larger SSDs.

As network card, I selected the TEG-10GECSFP (Version v2.0R) because it has the lowest maximal power consumption that I could find at 1.5W. I have no way to check this assertion accurately, however, given the fact that it is based on Aquantia NIC, my assumption is that 1.5W can be true only with a Directly Attached Cable (DAC) and not a Fiber Transceiver.

In case you are not interested in a 10GbE connection, you might still want to a NIC in the PCIe x16 slot, either to have a 2.5GbE which is now pretty common or some additional ethernet ports.

Depending on how many ports has the NIC, you could turn the server into a switch or at least have a dedicated management network.

The case I selected is the SilverStone SST-RM23-502-MINI which offers several interesting features:

1. several PSU compatibilities, so depending on your budget and needs, you can go from a standard ATX to a 2U redundant.
2. two 5.25" bays to expand the server in the future, for instance with a SATA case.
3. a 40cm short depth compared to standard rack.

The 2U rackmount case is the biggest constraint because it implies a low-profile CPU fan. For the CPU fan my choice goes to the Alpenfohn Black Ridge Low Profile because the reviews seem to agree on the fact that it provides the lowest temperature among the most popular low-profile CPU coolers in early 2023.
However, there a two important things to notice with this CPU fan:

1. the stock fan should be reversed as it blows in the wrong direction for most efficient setting. Contrarily to the reviews I read all over the internet, the stock fan of the Black Ridge is extremely silent such that I did not have to exchange it.
2. because of the dimensions, only Very Low Profile (VLP) RAM can be used and a fan above will not fit due to the case height.

In fact, the constraint of the VLP drastically restrict the landascape of compatible RAM sticks. I went for 2x32GB of Corsair Vengeance LPX (CMK64GX4M2A2666C16) which have the following timings: 16-18-18-35 (CL-tRCD-tRP-tRAS).

In terms of additional fans, as the case supports 2 standards 80mm x 25mm fans, I went for the classic and silent Noctua NF-A8 PWM.

Now that we have all our components, we can select the PSU accordingly the the peak power consumption. Because we want a server as silent as possible, a natural choice is a PicoPSU. With only 3 SSD, this build can use a 90W PicoPSU. However, if you plan to use a SATA cage with the remaining PCIe slot, I recommend to use a 160W PicoPSU.

The PicoPSU needs to be coupled with a power supply. For this build, I recommend a Leicke but the choice is rather wide and the important thing is to check the power to match the PicoPSU.

To be practical, note that the case supports only one 2.5" SSD and therefore, using a SSD mounting frame that can be mounted in a 5.25" bay such as the Sharkoon BayExtension.

Summary and cost

Here is a summary of the components with the approximate price as of September 2023. The final build is around 850€, accounting for delivery costs and additional costs such as thermal paste, etc.

The build in practice

The build is pretty simple to setup. The fact that we use a PicoPSU and microATX leaves a lot of space in the case for a good airflow. In particular, I appreciate this case because the additional fans are not obstructed by anything and the storage is really separated from the motherboard.

As I did not think about ordering a PCI power bracket with a barrel jack, I used some cable management clips. For an easier maintenance, I use a different color of SATA cable for the storage drives (black) and the OS drive (red). Ideally, a sticker with a proper number would be better, but with only 3 drives, the risk of confusion is low.

I installed the server into the 9U rack that is directly inserted in the desk of my office. The server has its independent socket with a smart plug for power consumption measurement. The 1GbE RJ45 is used for for the management network while the 10GbE SFP+ is used for everything else. Unfortunately, I could not use a DAC since the server is roughly 15 linear meters from the router which is the limit before signal degradation.

I have no tool to measure the noise but I am extremely impressed by how silent it is. Basically, the laptop on which I work on a daily basis is louder. The server is right next to me the whole day, every single day, and yet, I have never heard it. This is even surprising because many people recommend to avoid 2U for silent builds. My assumption is that the key is the PicoPSU for two reasons:

1. it takes little space which leaves a huge opening behind the case such that the airflow is really maximized,
2. the PicoPSU itself does not make any noise and does not heat up which in return lower the need for the fans to spin as fast as they might need to in a build with an ATX PSU.

Finally, here is a summary of the sensors from the smart plug that I use to measure the power consumption. It is not as rigorous as a Kill A Watt but it already give a good idea. The active power oscillates between 26W and 35W depending on the load.

The server is running on Proxmox and as shown on the screenshot below, it runs many services for such a low energy consumption and noise.

In conclusion, I am extremely satisfied with this build. It is not as low power as some of the ones in the spreadsheet attached at the beginning of the article, but it offers (relatively) modern hardware with a ton of capabilities, both for playing around and productionizing. In comparison to this database, I am using an active SFP+ connection which add some watts which is the price to pay for a 10GbE connection.
Despite being actively cooled, it is as silent as most passive builds which is a big plus for many homelabers, including me!

In summary, I could not recommend enough this build, although the pricepoint might be a bit high for some people. However, the case and the 10GbE NIC account for about 37% of the total cost, so if you don't need a 10GbE and can find an alternative case - maybe not even rack-mounted -, it might be possible to squeeze the budget to a mere 600€. In my opinion, both options are a bang for the buck!

What is a Homelab and why you should have one?

Technically, a homelab is a lab at home, or, in other words, a domestic infrastructure used for experimentation and learning. The difficulty to truly define what a homelab is comes from the fact that there are as many homelabs as homelabers, to serve as many purposes as one can imagine:

1. experimenting with networking or system administration
2. leveraging storage solutions to replace Google Drive, Box or the infamous copy on this external hard drive once in a while
3. replacing popular cloud tools and app with self-hosted alternatives
4. serving music, movies and other multimedia content to an entire home
5. automating home with a focus on privacy
6. learning professional skills using a pseudo-corporate environment
7. ...

In general, all these use cases share a couple of common goals: privacy and control. The privacy and control come from the fact that everything^[1] in your homelab runs on hardware that is physically located on your premises and for which you have a full access. In addition, the softwares that homelabs run are usually Free and Open-Source Softwares (FOSS) which remains the best guarantee when it comes to data collection and usage - and thus, privacy.

So in fact, a homelab does not have to be some expensive and highly specialized hardware, mounted in dozen of rack units, consumming hundreds of kilowatts per hour. A homelab can be as simple as a Raspberry Pi, or maybe an old re-purposed laptop, or a newly bought mini-PC that will be stand next to the TV and serves as a Home Theater PC (HTPC), or maybe all of these working together to offer different type of services!

What does it take to start?

Define your goals and expected maturity.

The most important at first is to set up goals to achieve in terms of functionalities to offer to the users and level of maturity.
By maturity, I mean the stability and reliability of the service: can the service breaks every now and then for an obscure reason or does the service need to be available 100% of the time to actually be worth it?

In my case, my homelab hosts several critical services:

• the entire home automation, which includes turning off and on the lights, fans, most home appliances but also the heating system,
• my company data and multiple services to help me with daily job.

These services need to be as reliable as any commercial offering.

On the other end of the spectrum, I run several services related to MLOps or experimenting with Large Language Models (LLMs) that could break in 15 minutes without any particular impact on anyone's life. This is just my playground!

Finally, some services are in between, such as the media players: no one will die from not being able to play a movie tonight, but it is still nice to have!

Build incrementally or per sub-function.

Of course, no one has to know perfectly all the different functionalities at first. A homelab can and, maybe even, should be built incrementally with the needs. Actually, it is very likely that one will start with a development environment, to learn, experiment and break, before migrating progressively towards a more stable and production-ready environment. There is not a single path to a homelab. One might even argue that a homelab is the path, literally.

You have an old laptop? Use it as a server for anything. Your services are stable enough and your users want more stability? Upgrade to a better hardware, maybe with uninterruptible power supply (UPC).

You have an old Rapberry Pi? Use it for Home Assistant and a couple of other services like a DNS server. You start to feel the hardware limitation with all your automations and devices in place? Scale to a mini PC.

Another strategy that might be more costly as an initial investment is to plan per sub-function rather than adding hardware from here and there. I personally identify three large sub-functions for any homelab:

1. network: what sort of network infrastrure do you need?
2. storage: how much and what sort of storage do you need? what is your backup and redundancy strategy?
3. compute: how much compute power do you need to run your services and applications?

Because the network is transversal, I would suggest to pay a particular attention to this sub-function at first. Changing a backbone infrastructure can become costly, especially if it implies to re-wire an entire house! On the other hand, for most people, a gigabyte or 2.5Gb internet might be more than enough and the focus should be on the actually router and switches.

The cost of a homelab: hardware, energy and time

The last question to answer is how much does it cost to have a homelab?. Again, this may vary from homelaber to homelaber because some would throw an unreasonable amount of money in what is their main hobby, while some are on a tight budget and do not want to spend more than they would on a commercial alternative.

There are several aspects to take into account:

• the initial investment
• the running cost, both in terms of time and electricity

It is easy to estimate the initial investment cost prior spending anything. On the other hand, estimating the running cost might be slightly more difficult as the overall price will depends on the consumption under idle and load, the usage schedule and the price of electricity where the homelab resides. In addition, the load curve of power supplies at low load - which is the main regime for a homelab - are almost never communicated by the manufacturers. In fact, the efficiency certifications focuses on high load and might be misleading for a homelaber. Chosing the right hardware to keep a balance between the compute power and low consumption is a topic by itself which unfortunately is one of the first to face when building a homelab.

I would also recommend to compute the cost marginally per service. For instance, Google Drive and Box Personal Pro are respectively 25 EUR/year and 108 EUR/year where I am located. Because I need beyond the capability of the free offer, I would anyway pay for one of these services. Since I am using my own service hosted in my homelab, I can say that I marginally save 25 to 108 EUR/year. This single service allows me to justify the energy consumption of my homelab. With a cost of 13.97 cents per kWh at the moment I write this, and about 40W for my storage node, the total annual cost of this particular service is roughly 45 EUR - privacy incuded -.

In the end, once combined and compared to the commercial alternatives, the cost of running a homelab can be vastly compensated by the marginal savings.

Finally, a factor that is vastly overlooked is time. How much time do I need to build and maintain a homelab? Regardless of how much you plan to be involved, the answer is: more than that. Running experimental services to learn by breaking things will obviously lead to spending a lot time fix things. Conversely, reaching a maturity level close to a commercial offering for critical services will require a tremendous amount of time ahead. Even with a homelab designed and implemented near professionally, issues will arise and when they come, it means solving them. This can sometimes be daunting or frustrating, especially when the something is critical like turning on the heater, and was working perfectly 10min ago.

Get a homelab!

Homelabs are great in the end. They can be simple or complex depending on your needs. There are always things to learn on the way, whether we build a service with high availability for professional purposes or a small experiment with old hardware.

Even with a high cost of electricity, they can help to save money on the long run when they are used to replace commercial applications. Otherwise, like any hobby, they require a variable investment in cash and time!