Lagrange trigonometric interpolating polynomial

The basic idea is the same as with the common Lagrange interpolating polynomial:
find a “polynomial” of exponentials or trigonometric functions which is periodic and
it equals 1 at one point of the grid and zero everywhere else

qj[x_]:=π(x-xj)/dx; (* auxiliary function to shorten expressions *)

(* Lagrange trigonometric interpolating polynomial as a sum of exponentials *)
hj[x_,n_]:=Sum[Exp[I(i-n/2)2qj[x]/(n+1)],{i,0,n}]/(n+1);
Print["Simplifying the sum of exponentials gives ",FullSimplify[hj[x,n]]]
Print["  and trying to convert it to trogonometric functions results in ",FullSimplify[ExpToTrig[hj[x,n]]]," :-("]

(* as I was not able to persuade Mathematica to transform it into the simplest form, I defined this simple version directly *)
myhj[x_,n_]:=Sin[qj[x]]/Sin[qj[x]/(n+1)]/(n+1);
Print["The simplest expression is ",myhj[x,n]]

(* checking that it is really the same result *)
Print["If both expressions are the same, we should get zero here: ",FullSimplify[ExpToTrig[hj[x,n]]-myhj[x,n]]]

Plotting the Lagrange trigonometric interpolating polynomial

(* grid *)
a=-1;b=2;nn=20;
h=(b-a)/(nn+1);
xx=Range[a,b-h,h];

(* which (Subscript[h, j])(x) to plot *)
jj=5; (* from 0 to nn *)
xxj=a+jj*h;
delta=If[#==xxj,1,0]&/@xx; (* δ-function on the grid *)
Show[
Plot[Evaluate[myhj[x,nn]/.{xj->xxj,dx->h}],{x,a,b},PlotRange->All,AxesLabel->{"x","(Subscript[h, j])(x)"}],
ListPlot[Transpose[{xx,delta}],PlotStyle->{Black}]
]

Derivatives of Lagrange trigonometric interpolating polynomial

First Derivative

dhj[x_,n_]:=Evaluate[FullSimplify[D[myhj[x,n],x]]]
Print["The first derivative of (Subscript[h, j])(x) is ", dhj[x,n]]

(* Evaluation at grid points *)
(* Diagonal elements *)
Print["Trying to calculate the diagonal elements of the differentiation matrix directly:"]
FullSimplify[dhj[xj,n]] (* this should not work *)
dmjj=Limit[dhj[x,n],x->xj];(* but a proper limit should *)
Print["Diagonal elements as a proper limit are ", dmjj]

(* Off-diagonal elements *)
dm=FullSimplify[Evaluate[dhj[xi,n]/.{xj->j*dx,xi->i*dx}],i∈Integers&&j∈Integers];
Print["Off-diagonal elements are ", dm]

(* a column of differentiation matrix contains derivatives of (Subscript[h, j])(x) at all grid points - beware of indices *)
dmj=Table[N[dm/.{i->ii-1,j->jj,n->nn,dx->h}],{ii,1,Length[xx]}]; 
dmj[[jj+1]]=dmjj/.{i->ii-1,j->jj,n->nn,dx->h};
Show[
Plot[{Evaluate[myhj[x,nn]/.{xj->xxj,dx->h}],Evaluate[dhj[x,nn]/.{xj->xxj,dx->h}]},{x,a,b},PlotRange->All,AxesLabel->{"x","f(x)"},PlotLegends->{"(Subscript[h, j])(x)","h(Subscript[', j])(x)"}],
ListPlot[{Transpose[{xx,delta}],Transpose[{xx,dmj}]},PlotStyle->{Black,Brown}]
]

Second Derivative

d2hj[x_,n_]:=Evaluate[FullSimplify[D[myhj[x,n],x,x]]]
Print["The second derivative of (Subscript[h, j])(x) is ", d2hj[x,n]]

(* Evaluation at grid points *)
(* Diagonal elements *)
Print["Trying to calculate the diagonal elements of the differentiation matrix directly:"]
FullSimplify[d2hj[xj,n]] (* this should not work *)
dm2jj=Limit[d2hj[x,n],x->xj];(* but a proper limit should *)
Print["Diagonal elements as a proper limit are ", dm2jj]

(* Off-diagonal elements *)
dm2=FullSimplify[Evaluate[d2hj[xi,n]/.{xj->j*dx,xi->i*dx}],i∈Integers&&j∈Integers];
Print["Off-diagonal elements are ", dm2]

(* a column of differentiation matrix dm2 contains the second derivatives of (Subscript[h, j])(x) at all grid points - again, beware of indices *)
dm2j=Table[N[dm2/.{i->ii-1,j->jj,n->nn,dx->h}],{ii,1,Length[xx]}]; 
dm2j[[jj+1]]=dm2jj/.{i->ii-1,j->jj,n->nn,dx->h};

Show[
Plot[{Evaluate[50hj[x,nn]/.{xj->xxj,dx->h}],10Evaluate[dhj[x,nn]/.{xj->xxj,dx->h}],Evaluate[d2hj[x,nn]/.{xj->xxj,dx->h}]},{x,a,b},PlotRange->All,AxesLabel->{"x","f(x)"},PlotLegends->{"(Subscript[h, j])(x) (50 x)","h(Subscript[', j])(x) (10 x)","h(Subscript['', j])(x)"}],
ListPlot[{Transpose[{xx,50delta}],Transpose[{xx,10dmj}],Transpose[{xx,dm2j}]},PlotStyle->{Black,Brown,Green}]
]