<< Statistics`DescriptiveStatistics`

 (*    Hermite spline interpolation of the Perfect Spline of degree 2*
        k      *) 
 (*    on [0, 1] with 2  
          k - 4 interior knots                                     *) 
 (*    generated as order statistics of 2  k - 
        4 i . i . d .  uniform (0, 
            1)            *) 
 (*    random variables                                                      *) 
 (*    Restriction :   k  ≥  
          3                                                   *) 

Print["k: (2k-1) degree spline interpolation"]
k = 4
Print["d = number of digits of accuracy"]
d = 30
f[x_, j_] := x^j
g[x_, y_, 
        k_] := ((x - y)*UnitStep[x - y])^(2*k - 
            1)
fp[x_, j_] := j*x^(j - 1)
gp[x_, y_, 
        k_] := (2*k - 
            1)*((x - y)*UnitStep[x - y])^(2*k - 
              2)

PS[x_, vec_List,  
        k_] := (1/Factorial[2*k])* (x^(2*k) + 
            2*Sum[(-1)^(j)*(x - vec[[j]])^(2*k)*
                  UnitStep[x - vec[[j]]],  {j,  1, 
                  2*k - 4}])
PSp[x_, vec_List, 
        k_] := (2*k)*(1/Factorial[2*k])*(x^(2*k - 1)  + 
            2*Sum[(-1)^(j)*(x - vec[[j]])^(2*k - 1)*
                  UnitStep[x - vec[[j]]],  {j,  1, 
                  2*k - 4}])

Basis[x_, vec_List, k_] := 
      Join[{1}, Table[f[x, j], {j, 1, 2*k - 1}], 
        Table[g[x, vec[[j]], k], {j, 2, 
            2*k - 3}]]
HermiteInterp[x_, coeff_List, vec_List, k_] := 
        coeff . Basis[x, vec, k];
mxlist[coeff_List, vec_List,  k_] := 
      Table[HermiteInterp[i/1000, coeff, vec, k] - PS[i/1000, vec,  k], {i, 
          1, 1000}]

Print["Interior knots"]
VectauUnif = 
      Sort[Table[
          Random[Real, {0, 1}, 40], {j, 1, 
            2*k - 4}]]
Print["All knots"]
Vectau = 
        Append[Prepend[VectauUnif, 0], 1];
MM = 
        Join[
Table[
            Join[{1}, 
Table[
                f[Vectau[[j]], i], {i, 1, 2*k - 1}], 
              Table[g[Vectau[[j]], Vectau[[i]], k], {i, 2, 
                  2*k - 3}]], {j, 1, 2*k - 2}], 
Table[
            Join[{0, 1}, 
Table[
                fp[Vectau[[j]], i], {i, 2, 2*k - 1}], 
              Table[gp[Vectau[[j]], Vectau[[i]], k], {i, 2, 
                  2*k - 3}]], {j, 1, 2*k - 2}]];
VecPS = 
        Join[Table[PS[Vectau[[j]], VectauUnif, k], {j, 1, 2*k - 2}], 
          Table[PSp[Vectau[[j]], VectauUnif, k], {j, 1, 
              2*k - 2}]];
HermiteCoeff = LinearSolve[MM, VecPS];
mtl = mxlist[HermiteCoeff, Vectau, k];
Print["Maximal interpolation error"]
mx = Max[mtl]
Print["plot of the perfect spline"] 
P1  =  
      Plot[PS[x, VectauUnif, k], {x, 0, 1}, 
        PlotRange -> All]
Print["plot of interpolant"]
P2 = 
      Plot[HermiteInterp[x, HermiteCoeff, Vectau, k], {x, 0, 1}, 
        PlotRange -> All]
Print["plot of the perfect spline and its interpolant"]
Show[P1, P2]
Print["Plot of the interpolation error"]
P3=Plot[
        HermiteInterp[x, HermiteCoeff, Vectau, k] - PS[x, VectauUnif, k], {x, 
          0, 1},   PlotRange -> 
          All]