bernstein.tex

\begin{flushright} {\tiny {\color{gray} bernstein.tex}} \end{flushright}


The $n+1$ Bernstein basis polynomials of degree $n$ on the interval $[0,1]$
are defined as \footnote{\url{https://en.wikipedia.org/wiki/Bernstein_polynomial}}
\index{general}{Bernstein Polynomials}
\[
b_{m,n}(x) = \left( \begin{array}{c} n \\ m \end{array}\right) x^m(1-x)^{n-m}
\qquad m=0,1,...n
\]
The first few Bernstein polynomials are 
\begin{eqnarray}
b_{0,0}(x) &=& 1 \\
b_{0,1}(x) &=& 1-x \nn\\
b_{1,1}(x) &=& x \\
b_{0,2}(x) &=& (1-x)^2 \nn\\
b_{1,2}(x) &=& 2x(1-x) \nn\\
b_{2,2}(x) &=& x^2 
\end{eqnarray}

\includegraphics[width=5cm]{images/bernstein/b0.pdf}
\includegraphics[width=5cm]{images/bernstein/b1.pdf}
\includegraphics[width=5cm]{images/bernstein/b2.pdf}

We see that the zero-th and first order polynomials are the same as the linear basis functions defined in 
Section~\ref{sec:elts1D}. However the second order polynomials (and higher) differ from the second-order
basis functions. 

Also, the Bernstein polynomials have a lot of properties, but one that is of importance to us
is the following: $b_{m,n}(x) \geq 0 \quad  \forall x\in [0,1]$, i.e. the polynomials 
are positive. This is however not true for basis functions for $n\geq 2$.
Another important property shared with basis functions is that their sum over the interval is 
exactly 1, i.e. $\sum_m b_{m,n}(x)=1$.

In order to facilitate the comparison between the 2nd-order basis functions and Bernstein 
polynomials, I will express the latter as a function of the reduced coordinate
$r\in[-1,1]=2(x-1/2)$ (or $x=(r+1)/2$). We have then:

\begin{eqnarray}
b_{0,2}(r) &=& \frac{1}{4}(1-r)^2 \nn\\
b_{1,2}(r) &=& \frac{1}{2}(1-r^2) \nn\\
b_{2,2}(r) &=& \frac{1}{4}(1+r)^2
\end{eqnarray}

Both 2nd-order Bernstein polynomials and basis functions are plotted here under:
\begin{center}
\includegraphics[width=7cm]{images/bernstein/b2_.pdf}
\includegraphics[width=7cm]{images/bernstein/N2_.pdf}\\
{\captionfont Left: Second-order Bernstein polynomials; right: 2nd-order basis functions.}
\end{center}

Having reached this point, the burning question is why should we care?

\paragraph{Example 1}
In order to answer this question, let us carry out the following 
experiment: each node $i$ in the element carries a field value $f_i$ and for simplicity, 
we choose $f_0=f(r=-1)=1$, $f_1=f(r=0)=0$, $f_2=f(r=+1)=0$.
Then, we can compute the value of the field inside of the element 
as we usually do in the FE methodology:
\[
f^h(r) = \sum_{i=0}^2 \bN_i(r) f_i = f_0 \bN_0(r) = \bN_0(r) = \frac{1}{2}r(r-1)
\]
This means that although the field $f$ is always positive (or null) inside the element
its representation with the basis functions is negative over half (!) of the element
(see purple curve on the right panel above).
If we now turn to the Bernstein polynomials:
\[
f^h(r) = \sum_{i=0}^2 b_{i,2} f_i = b_{0,2} = \frac{1}{4}(1-r)^2
\]
which is {\it always} positive over the interval $[-1,+1]$, 
and looking at the purple curve on the left panel above, 
we see that the value decreases monotonously when we go away from node 1, and reaches 
zero at the other end of the element. 

Also:
\[
\text{Shape function: } \int_{-1}^{+1} f^h(r)dr = \int_{-1}^{+1} \frac{1}{2}r(r-1) dr = \frac{1}{3}  
\]
\[
\text{Bernstein polynomial: } \int_{-1}^{+1} f^h(r)dr =  \int_{-1}^{+1} \frac{1}{4}(1-r)^2 dr=  \frac{2}{3}  
\]
Analytical value for the integral can be obtained by splitting the integral as $\int_{-1}^0 + \int_0^{+1}$.
The left part can be represented by a line of equation $-r$ and the right part simply by 0, so that the 
integral is equal to 1/2. We then see that the Shape function-based interpolation underestimates
the integral while the Bernstein polynomial-based interpolation overestimates it.


\paragraph{Example 2}

We now choose $f_0=f(r=-1)=1$, $f_1=f(r=0)=1$, $f_2=f(r=1)=0$. Then 
\[
f^h(r) 
= \sum_{i=0}^2 f_i \bN_i(r) 
= f_0 \bN_0(r) + f_1 \bN_1(r) 
= \bN_0(r) + \bN_1(r) 
= \frac{1}{2}r(r-1) + 1-r^2 = -\frac{1}{2}r^2 -\frac{1}{2}r +1  
\]
Looking now at the Bernstein polynomials:
\[
f^h(r) 
= \sum_{i=0}^2 b_{i,2} f_i 
= b_{0,2} + b_{1,2} = \frac{1}{4}(1-r)^2 + \frac{1}{2}(1-r^2)
= -\frac{1}{4}r^2-\frac{1}{2}r+\frac{3}{4}
\]
If we now plot both approximations:
\begin{center}
\includegraphics[width=7cm]{images/bernstein/N2__.pdf}
\end{center}
We see that in this case the $Q_2$ basis function-based approximation yields values 
above 1 over half of the element while the Bernstein polynomial-based 
approximation remains between 0 and 1 as expected.

\paragraph{Approximation of polynomials}
Let us now explore another aspect of such an interpolation based on the Bernstein polynomials
and assume that $f(r)=C$. Then 
\[
f^h(r) = \sum_{i=0}^2 f_i b_{i,2}(r) = C \sum_{i=0}^2 b_{i,2}(r) = C \cdot 1 = C
\]
Such interpolation can exactly represent a constant field. 
Let us assume that $f(r)=ar+b$. Then 
\begin{eqnarray}
f^h(r) 
&=& \sum_{i=0}^2 f(r_i) b_{i,2}(r)  \nn\\
&=& \sum_{i=0}^2 (ar_i+b) b_{i,2}(r) \nn\\
&=& a\sum_{i=0}^2 r_i b_{i,2}(r) + b \sum_{i=1}^3 b_{i,2}(r) \nn\\
&=& a (-b_{0,2}(r)+b_{2,2}(r)) + b \cdot 1 \nn\\
&=& ar+b 
\end{eqnarray}
Such interpolation can exactly represent a linear field. 

Let us assume that $f(r)=ar^2+br+c$. Then 
\begin{eqnarray}
f^h(r) 
&=& \sum_{i=0}^2 f(r_i) b_{i,2}(r) \nn \\
&=& \sum_{i=0}^2 (ar^2+br+c) b_{i,2}(r) \nn\\
&=& \sum_{i=0}^2 r_i^2 b_{i,2}(r) + b \sum_{i=0}^2 r_i b_{i,2}(r) + c \sum_{i=0}^2 b_{i,2}(r) \nn\\
&=& a\sum_{i=0}^2 r_i^2 b_{i,2}(r) + br + c \nn\\
&=& a (b_{0,2}(r) + b_{2,2}(r)) + br + c \nn\\
&=& a\frac{1}{2}(1+r^2)  + br + c 
\end{eqnarray}
which is not equal to $f(r)$.

On the other hand it is trivial to show that 
\begin{eqnarray}
f^h(r) 
&=&  \sum_{i=0}^2 f(r_i) \bN_i(r) \nn \\
&=&  \sum_{i=0}^2   (ar_i^2+br_i+c)  \bN_i(r) \nn \\
&=&  a \sum_{i=0}^2  r_i^2 \bN_i(r) + b \sum_{i=0}^2 r_i  \bN_i(r) +  c \sum_{i=0}^2  \bN_i(r) \nn \\
&=&  a \sum_{i=0}^2  r_i^2 \bN_i(r) + b \sum_{i=0}^2 r_i  \bN_i(r) +  c \nn\\
&=&  a (N_0(r)+N_2(r)) + b (-\bN_0(r) +\bN_2(r)) +  c \nn\\
&=&  a r^2  + b r+ c 
\end{eqnarray}

To hammer the point once more: let $f(r)=r^2+r+1$.
Then 
\begin{eqnarray}
f^h_{Q_1} 
&=& f(-1) \frac{1}{2}(1-r) + f(+1) \frac{1}{2}(1+r) \nn\\
&=&  \frac{1}{2}(1-r) + 3 \frac{1}{2}(1+r) \nn\\
&=& 2-r \\
f^h_{Q_2} &=& r^2+r+1 \\
f^h_{B_2} &=& \frac{1}{2} (1+r^2)+r+1 
\end{eqnarray}
We see on the following figure that although Bernstein polynomials cannot 
represent $f(r)$ exactly they still do a better job than first order basis functions
($Q_1$ projection).
\begin{center}
\includegraphics[width=7cm]{images/bernstein/hammer.pdf}
\end{center}

As a conclusion there is a trade-off: 2nd-order Bernstein polynomials {\it always} yield positive 
values when the field is positive (as opposed to 2nd-order basis functions) but they cannot 
represent exactly a 2nd-order polynomial field (while basis functions can).

The positivity can be really critical in geodynamical simulations: a negative density makes no sense, 
and a negative viscosity even less!

%...............................................................................
\paragraph{Example 3}

The 2nd-order Bernstein polynomials are used in \stone~64. The actual context of this stone is not 
important. Fields such as density and viscosity are known on the 
9 nodes of the $Q_2$ element and need to be projected onto the 9 quadrature points. 
For instance, these nodal fields are given by:
\begin{center}
\includegraphics[width=7cm]{images/bernstein/rhonodal.png}
\includegraphics[width=7cm]{images/bernstein/etaeffnodal.png}
\end{center}
The resulting fields on the quadrature points are shown:
\begin{center}
\includegraphics[width=7cm]{images/bernstein/qrho.png}
\includegraphics[width=7cm]{images/bernstein/qetaeff.png}
\end{center}
The bottom row is obtained with the basis functions while the top row is obtained with the Bernstein 
polynomials as interpolants. The thin blue line actually indicated points with negative viscosity and 
on the left the colour bar shows densities below the value of 1890 (lowest density in the domain).

%...............................................................................
\paragraph{Example 4}

In this example I have created a small python code (in {\tt /images/bernstein})
which creates a single $Q_2$ element in 2D, i.e. it consists of 9 nodes.
9 different scenarios are created by assigning some of the nodes either a 
density value of 1 or 0. 
The internal numbering is as follows:
\begin{verbatim}
3--6--2
|     |
7  8  5
|     |
0--4--1
\end{verbatim}

I then proceed to generate 100,000 random points inside the element
and I use either the $Q_2$ basis functions to interpolate the density onto them
or the 2nd order Bernstein polynomials.
The basis functions used are as follows:
\begin{lstlisting}
def NNN(r,s):
    val = np.zeros(9,dtype=np.float64)
    val[0]= 0.5*r*(r-1.) * 0.5*s*(s-1.)
    val[1]= 0.5*r*(r+1.) * 0.5*s*(s-1.)
    val[2]= 0.5*r*(r+1.) * 0.5*s*(s+1.)
    val[3]= 0.5*r*(r-1.) * 0.5*s*(s+1.)
    val[4]=    (1.-r**2) * 0.5*s*(s-1.)
    val[5]= 0.5*r*(r+1.) *    (1.-s**2)
    val[6]=    (1.-r**2) * 0.5*s*(s+1.)
    val[7]= 0.5*r*(r-1.) *    (1.-s**2)
    val[8]=    (1.-r**2) *    (1.-s**2)
    return val 

def BBB(r,s):
    val = np.zeros(9,dtype=np.float64)
    val[0]= 0.25*(1-r)**2  * 0.25*(1-s)**2  
    val[1]= 0.25*(1+r)**2  * 0.25*(1-s)**2  
    val[2]= 0.25*(1+r)**2  * 0.25*(1+s)**2  
    val[3]= 0.25*(1-r)**2  * 0.25*(1+s)**2  
    val[4]= 0.5*(1-r**2)   * 0.25*(1-s)**2  
    val[5]= 0.25*(1+r)**2  * 0.5*(1-s**2) 
    val[6]= 0.5*(1-r**2)   * 0.25*(1+s)**2  
    val[7]= 0.25*(1-r)**2  * 0.5*(1-s**2) 
    val[8]= 0.5*(1-r**2)   * 0.5*(1-s**2) 
    return val 
\end{lstlisting}


 
In the following table I report the min/max values of the density on the 
swarm of points for both methods and all 9 cases:

\begin{center}
\begin{tabular}{lllll}
\hline
case & min($\rho^h_{Q_2}$) & max($\rho^h_{Q_2}$)  & min($\rho^h_{B_2}$) & max($\rho^h_{B_2}$)  \\ 
\hline
\hline
1 &-0.124978 &0.994484 &0.000000 &0.996320  \\
2 &-0.124999 &0.998679 &0.000000 &0.499537  \\
3 &-0.140625 &1.124868 &0.000000 &0.999317  \\
4 &-0.125000 &0.999997 &0.000000 &0.999998  \\
5 &-0.187699 &1.124808 &0.000000 &0.999978  \\
6 & 0.000000 &1.265620 &0.000000 &0.999980  \\
7 &-0.225000 &1.124911 &0.000000 &0.999996  \\
8 &-0.124418 &1.224999 &0.000000 &1.000000  \\
9 & 1.000000 &1.000000 &1.000000 &1.000000  \\
\hline
\end{tabular}
\end{center}

We find that using $Q_2$ basis functions yields values 
either above 1 or less than 0, which makes little sense, 
while the $B_2$ functions produce values exactly between 0 and 1. 
Only case 2 with $B_2$ produces values up to 0.5, which is understandable 
since the function $B_{2,7}$ (i.e. at node 7) reaches a maximum value of 1/2.

Also, the last case 9 corresponds to assigning all nodes $\rho=1$,
so that we expect the min/max values to be exactly 1 (sanity check).

I show below from left to right: the nodal values, the $Q_2$-based field,
where it is negative, and the $B_2$-based field for all nine cases:

\begin{itemize}
\item Case 1:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0000.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_0.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_0.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_0.png}\\
\end{center}
\item Case 2:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0001.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_1.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_1.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_1.png}\\
\end{center}
\item Case 3:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0002.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_2.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_2.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_2.png}\\
\end{center}
\item Case 4:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0003.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_3.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_3.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_3.png}\\
\end{center}
\item Case 5:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0004.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_4.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_4.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_4.png}\\
\end{center}
\item Case 6:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0005.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_5.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_5.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_5.png}\\
\end{center}
\item Case 7:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0006.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_6.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_6.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_6.png}\\
\end{center}
\item Case 8:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0007.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_7.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2neg_7.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_7.png}\\
\end{center}
\item Case 9:
\begin{center}
\includegraphics[width=4cm]{images/bernstein/nodes0008.png}
\includegraphics[width=4cm]{images/bernstein/rhoQ2_8.png}
\includegraphics[width=4cm]{images/bernstein/rhoB2_8.png}\\
\end{center}
\end{itemize}