Assume, that we have points $x_i$ with $i=1,...,N+1$. We construct the Lagrange basis polynomials as \begin{align} L_j(x) = \prod_{k\not = j} \frac{x-x_k}{x_j-x_k} \end{align}

Now according to my computation and the results by Yves Daoust here, the derivative of $L_i$ can be computed as \begin{align} L'_j(x) = L_j(x)\cdot \sum_{k\not = j}\frac{1}{x_k-x_j} \end{align}

I try to reproduce the numerical results of a paper, and for this results the authors use the derivative matrix $D$ with $D_{ij} =L_i'(x_j)$.

The authors use the Legendre Gauss Lobatto quadrature points, plotted below. Now I can easily construct the basis polynomials, also plotted below, together with the quadrature points.

Given the concrete choice of basis points, the plots suggest, that $L'_j(x_j)=0$ except for the first and last basis function. Additionally it seems, that $L'_j(x_i)\not =0$ for $i\not = j$. But using the derived formula from above, we obtain \begin{align} L'_j(x_i) = L_j(x_i)\cdot \sum_{k\not = j}\frac{1}{x_k-x_j}=0 \end{align} for $i\not = j$ since $L_j(x_i) =\delta_{ij}$, which contradicts my observation.

I used the following Matlab functions, to construct the matrix $D$.

```
function y=dl(i,x,z)
n = length(x);
y = 0;
for m=1:n
if not(m==i)
y = y + 1/(x(m)-x(i));
end
end
size(y)
y = y*l(i,x,z);
end
function y=l(i,x,z)
n = length(x);
% computes h_i(z)
y = 1;
for m=1:n
if not(m==i)
y = y.*(z-x(m))./(x(i)-x(m));
end
end
end
```

Where `dl`

and `l`

are $L'$ and $L$. $D$ is then constructed by

```
D = zeros(M+1,M+1);
for i=1:M+1
for j=1:M+1
D(i,j) =dl(i,X,X(j));
end
end
```

which gives the matrix

```
D =
10.5000 0 0 0 0 0 0
0 -0.0000 0 0 0 0 0
0 0 0.0000 0 0 0 0
0 0 0 -0.0000 0 0 0
0 0 0 0 0.0000 0 0
0 0 0 0 0 0.0000 0
0 0 0 0 0 0 -10.5000
```

Here I agree with the zeros on the diagonal except for the first and last element. The zeros everywhere else agree with the formula, but they don't agree with my observation.

If I use finite differences (which is no solution for the real implementation) in the following way: ` FD(i,j) =(l(i,X,X(j)+eps)-l(i,X,X(j)-eps))/(2*eps);`

I obtain the output

```
FD =
-10.5000 -2.4429 0.6253 -0.3125 0.2261 -0.2266 0.5000
14.2016 0.0000 -2.2158 0.9075 -0.6164 0.6022 -1.3174
-5.6690 3.4558 0.0000 -2.0070 1.0664 -0.9613 2.0500
3.2000 -1.5986 2.2667 0 -2.2667 1.5986 -3.2000
-2.0500 0.9613 -1.0664 2.0070 -0.0000 -3.4558 5.6690
1.3174 -0.6022 0.6164 -0.9075 2.2158 -0.0000 -14.2016
-0.5000 0.2266 -0.2261 0.3125 -0.6253 2.4429 10.5000
```

which is zero on the diagonal (except first and last) and nonzero everywhere else.

**So here are my questions:**

Did I make an error in the implementation?

Is my formula derived under the assumption that $x$ is not a basis point?

Can I modify my code to obtain the results that my finite difference code produces?

Is my assumption correct, that rather the results of the finite difference computation are "correct"?

**EDIT**

It seems, like the derivation of the formula goes wrong for $x=x_i$. If I compute the derivative without simplification, I get \begin{align} L_j'(x) = \sum_{l\not = j} \frac{1}{x_j-x_l}\prod_{m\not = (j,l)} \frac{x-x_m}{x_j-x_m} \end{align} Or in code

```
function y = alternative_dl(j,x,z)
y = 0;
n = length(x);
for l=1:n
if not(l==j)
k = 1/(x(j)-x(l));
for m=1:n
if not(m==j) && not(m==l)
k = k*(z-x(m))/(x(j)-x(m));
end
end
y = y + k;
end
end
end
```

Which agrees with the finite difference computation.

So it seems to me, that simplifying the above formula includes some "hidden division by zero" if $x=x_i$.