# General Lagrange basis formula (usual problem in finite element context)

It is easy to prove that,

$$\{p_1(x,y)=1-x-y\;,\;p_2(x,y)=x\;,\;p_3(x,y)=y\}$$

is a Lagrangian basis of $\mathbb{P}_1(\hat{T})$ (polynomials of total degree less that 1 living on $\hat T$), where $\hat T$ is the triangle with vertices $a_1=(0,0)$, $a_2=(1,0)$ and $a_3(0,1)$. It is clear that $p_i(a_j)=1$ if $i=j$ and $p_i(a_j)=0$ if $i\neq j$.

I also can calculate a basis function generating $\mathbb{P}_2(\hat T)$. But calculate basis functions of $\mathbb{P}_3(\hat T)$, $\mathbb{P}_4(\hat T)$ it is too dificult, so I think that there exists a general formula to write in my finite element code.

There exists a formula to compute the basis function sets of $\mathbb{P}_k(\hat T)$, with $k$ any integer?

I hope I'm correct in interpreting your question as asking for a constructive approach for generating the Lagrange bases on the reference triangle, apologies if this is more simplistic than you want.

Let's introduce a third (dependent) spatial variable, $\zeta = 1 - x -y$. Then your basis is just the set ${\zeta,x,y}$. In the $P_2$ case we're free to place the node points at the vertices of the triangle, and the midpoints of the edges. We can identify the corner nodes with

• a quadratic function which is 0 at $\zeta=0,\frac{1}{2}$ (i.e proportional to $\zeta (\zeta-\frac{1}{2}$))
• a quadratic function which is 0 at $x=0,\frac{1}{2}$ (proportional to $x (x-\frac{1}{2}$))
• a quadratic function which is 0 at $y=0,\frac{1}{2}$ (proportional to $y (y-\frac{1}{2}$))

Meanwhile the edge points can be identified with

• a quadratic function which is zero on $x = 0$ and $y=0$ (ie. proportional to $xy$)
• a quadratic function which is zero on $\zeta = 0$ and $y=0$ (ie. proportional to $\zeta y$)
• a quadratic function which is zero on $x = 0$ and $\zeta=0$ (ie. proportional to $\zeta y$)

With $P_3$ we just extend the concept again; a corner point takes a cubic function zero on $\zeta =0,\frac{1}{3},\frac{2}{3}$ hence proportional to $\zeta(\zeta-\frac{1}{3})(\zeta-\frac{2}{3})$, with the others cycling through. Now for the others, we carry on with a process of elimination; for example there are 3 points not on the lines $\zeta =0$ or $\zeta=1/3$. Two of these lie on the line $x=0$ (similarly $y=0$) hence the remaining point takes a function proportional to $\zeta (\zeta-1/3)x$ [or $\zeta (\zeta-1/3)y$ for its friend. Symmetric rotation takes us through all points except the internal point. This requires a function which vanishes on $\zeta = 0$, $x=0$ and $y=0$, thus proportional to $xy\zeta$.

In this way finding the functions reduces to an essentially graphical exercise of finding $k$ lines through all but one of $T_{k+1}$ points arranged in a triangle. Calculating the normalisation constants is left as an exercise for the reader.

There is no closed formula for the general case, but there are standard procedures. I'm assuming you are looking for triangular Lagrange elements, but the adaption to other cases is straightforward. Say you are looking for the nodal basis $\{\psi_1,\dots,\psi_m\}$ of $P_k(T)$, where $T$ is the reference triangle with vertices $(0,0)$, $(0,1)$, and $(1,0)$, and you have $m:=\frac12(k+1)(k+2)$ (the dimension of $P_k(T)$) different nodes $\xi_1,\dots,\xi_m\in T$ (where $\xi_i=(x_i,y_i)$). A nodal basis satisfies by definition $$\label{1}\psi_i(\xi_j) = \delta_{ij} = \begin{cases} 1 & \text{if } i=j, \\ 0 & \text{else.}\end{cases} \tag{1}$$ If the nodes are chosen correctly, this uniquely determines the $\psi_i$. (In fact, this condition is exactly what determines how many nodes have to lie on the different vertices, edges, and the interior of $T$). This means that you can use this condition to solve for the unknown coefficients of the $\psi_i$.

One way (an alternative way is given by origimbo) to do this is to write $$\psi_i(\xi_j) = \sum_{n=1}^m a_{in} p_n(\xi_j),$$ where $p_n$ are the monomials spanning $P_k(T)$, i.e., $1,x,y,x^2,xy,y^2,\dots$ Inserting this into \eqref{1} leads to the following system of linear equations: $$\begin{pmatrix} p_1(\xi_1) & p_2(\xi_1) &\dots \\ p_1(\xi_2) & p_2(\xi_2) & \dots\\ \vdots & \vdots &\ddots \end{pmatrix} \begin{pmatrix} a_{11}& a_{21} &\dots\\ a_{12} & a_{22} & \dots \\\vdots & \vdots &\ddots \end{pmatrix} = \begin{pmatrix} 1& 0 &\dots\\ 0 & 1 & \dots \\\vdots & \vdots &\ddots \end{pmatrix}$$ (The first matrix is called Vandermonde matrix.) This system you can solve for the unknown $a_{ij}$ either symbolically (e.g., using SymPy or Matlab's Symbolic Computing Toolbox -- no need to do this by hand!) or numerically (e.g., using Matlab's backslash, although Vandermonde matrices can get quite ill-conditioned for large $k$ and arbitrary $\xi_j$, although there are choices that minimize the condition number).