<< Statistics`DescriptiveStatistics`

 (*    Hermite spline interpolation of f_  
          0 (x)  =  
        x^(2*k)            *) 
 (*                          PART1:                                  *) 
 (*         on [0, 1] with equidistant 2  k - 
        4 interior knots            *) 
 (*    Restriction :   k  ≥  
          3                                           *) 

Print["k: (2k-1) degree spline interpolation"]
k = 4


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)
f0[x_, k_] :=  x^(2*k)
f0p[x_, 
        k_] := (2*
            k)*x^(2*k - 
              1)

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[Abs[HermiteInterp[i/1000, coeff, vec, k] - f0[i/1000,  k]], {i, 
          1, 1000}]


Print["Equidistant interior knots "]
VectauEqui = 
      Table[j/(2*k - 3), {j, 1, 2*k - 4}]
Print["All knots"]
Vectau = Append[Prepend[VectauEqui, 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}]];
VecMS = 
        Join[Table[f0[Vectau[[j]], k], {j, 1, 2*k - 2}], 
          Table[f0p[Vectau[[j]], k], {j, 1, 
              2*k - 2}]];
HermiteCoeff = LinearSolve[MM, VecMS];
mtl = mxlist[HermiteCoeff, Vectau,  k];
Print["Maximum interpolation error"]
mx = Max[mtl]
Print["Value given by the Bernoulli numbers formula"]
bernform =  (2/(2*k - 3)^(2*k))*
        Abs[(1 - 2^(-2*k))*
            BernoulliB[2*k]]
P1  =  
      Plot[f0[x, 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 function and interpolant"]
Show[P1, P2]
Print["Plot of the monospline of degree 2*k"]
P3=Plot[
        HermiteInterp[x, HermiteCoeff, Vectau, k] - f0[x, k], {x, 0, 
          1},   PlotRange -> 
          All]

<< Statistics`DescriptiveStatistics`

 (*    Hermite spline interpolation of f_  
          0 (x)  =  
        t^(2*k)            *) 
 (*                          PART2:                                  *) 
 (*    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)
f0[x_, k_] :=  x^(2*k)
f0p[x_, 
        k_] := (2*
            k)*x^(2*k - 
              1)

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[Abs[HermiteInterp[i/1000, coeff, vec, k] - f0[i/1000,  k]], {i, 
          1, 1000}]

Print["Interior knots "]
VectauUnif = 
      Sort[Table[
          Random[Real, {0, 1}, d], {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}]];
VecMS = 
        Join[Table[f0[Vectau[[j]], k], {j, 1, 2*k - 2}], 
          Table[f0p[Vectau[[j]], k], {j, 1, 
              2*k - 2}]];
HermiteCoeff = LinearSolve[MM, VecMS];
mtl = mxlist[HermiteCoeff, Vectau,  k];
Print["Maximum interpolation error"]
mx = Max[mtl]
Print["plot of function being interpolated, f_0 (x) = x^(2*k)"]
P1  =  
      Plot[f0[x, 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 function and interpolant"]
Show[P1, P2]
Print["Plot of the monospline (the interpolation error) of degree 2*k\
"]
P3=Plot[
        HermiteInterp[x, HermiteCoeff, Vectau, k] - f0[x, k], {x, 0, 
          1},   PlotRange -> 
          All]