Fitting Data by Cubic Splines

SplinesJanHurt.nb
Reference: Späth, Helmuth: Algorithmen für elementare Ausgleichsmodelle. Oldenbourg Verlag, München Wien 1973 (pp. 59-61). Beware of a lot of errors in the formulas!!!

Notice that with size n of the input data the system of 3n linear equations must be solved. Generally, this is simple task for Mathematica since the matrix of the system is five-diagonal, but rendering graphics of the spline functions usually takes a long time. For PC Pentium 133MHz with 32MB memory for n = 200 it takes some 10 minutes. 

Compare also with the method of runnig medians (Tukey).

splinesJH::usage = "splinesJH[x, u, p, z] gives a list of cubic splines in the respective intervals given by the list x of independent observations, u being the list of dependent observations at points x, weights list p, and z as the argument of the spline functions. If p->∞, the method leads to usual numerical interpolation while for p->0 leads to the least squares method."

plotsplines::="plotsplines[x, u, p]  produces the plot of the spline function given the list x of independent variables, the list u of dependent variables with weights list p."

\(TraditionalForm
    \`splinesJH[x_, \ u_, \ p_, \ z_]\  := \ 
\ \ 
      Module[{pp\  = \ p, \ dz, \ dyy, \ yy2, \ eq1, \ eq2, \ eq3, \ \ eq4, \ 
          eq5, \ variables, \ sol, \ yyy, \ yyy2, \ dyyy, \ 

          \ \ \ \ \ \ \ \ \ \ dyyy2, aa, \ bb, \ cc, \ dd, \ xk, \ n}, \ 

        \ \ \ n\  = \ Length[x]; \ 
        dx\  = \ Take[RotateLeft[x]\  - \ x, \ n\  - \ 1]; \ 
\ \ \ \ 
        yy\  = \ Array[y, \ n]; \ 
\ \ \ \ 
        dyy\  = \ Take[RotateLeft[yy]\  - \ yy, \ n\  - \ 1]; \ 
\ \ \ \ 
        yy2\  = \ Array[y2, \ n]; \ tt\  = \ Array[t, \ n]; \ 
\ \ \ \ 
        eq1\  = \ 
          Table[dx[\([k\  - \ 1]\)]*yy2[\([k\  - \ 1]\)]\  + \ 
                dx[\([k]\)]*yy2[\([k\  + \ 1]\)]\  + \ 
\ \ \ \ \ \ \ \ 
                2*\((dx[\([k\  - \ 1]\)]\  + \ dx[\([k]\)])\)*
                  yy2[\([k]\)]\  == \ 
\ \ \ \ \ \ \ 
              6*\((dyy[\([k]\)]/dx[\([k]\)]\  - 
                    dyy[\([k\  - \ 1]\)]/dx[\([k\  - \ 1]\)])\), \ 

            \ \ \ \ \ \ {k, \ 2, \ n\  - \ 1}]; 
        eq2\  = \ {yy2[\([1]\)]\  == \ 0, \ yy2[\([n]\)]\  == \ 0}; \ 

        \ \ \ \ eq3\  = 
          Table[pp[\([k]\)]*\((u[\([k]\)]\  - \ yy[\([k]\)])\)\  == \ 
              tt[\([k]\)], \ {k, \ n}]; \ 
\ \ \ \ 
        eq4\  = {
            tt[\([1]\)]\  == \ 
              \((yy2[\([2]\)]\  - \ yy2[\([1]\)])\)/
\ \ \ \ \ \ \ \ 
                dx[\([1]\)], \ 
            tt[\([n]\)]\  == 
					
              \(-\((\((yy2[\([n]\)]\  - \ yy2[\([n\  - \ 1]\)])\)/
                    dx[\([n\  - \ 1]\)])\)\)}; 
        eq5\  = \ 
          Table[tt[\([k]\)]\  == 
              \((yy2[\([k\  + \ 1]\)]\  - \ yy2[\([k]\)])\)/dx[\([k]\)]\  - \ 
                
\ \ \ \ \ \ \ \ 
                \((yy2[\([k]\)]\  - \ yy2[\([k\  - \ 1]\)])\)/
                  dx[\([k\  - \ 1]\)]\  == \ 
              yy2[\([k\  - \ 1]\)]/dx[\([k\  - \ 1]\)]\  + \ 
                yy2[\([k\  + \ 1]\)]/dx[\([k]\)]\  - \ 
\ \ \ \ \ \ \ \ 
                \((1/dx[\([k]\)]\  + \ 1/dx[\([k\  - \ 1]\)])\)*
                  yy2[\([k]\)], {k, \ 2, \ n\  - \ 1}]; \ 
\ \ \ \ 
        variables\  = \ Join[yy, \ yy2, \ tt]; \ 
\ \ \ \ 
        sol\  = \ 
          Solve[Join[eq1, \ eq2, \ eq3, \ eq4, \ eq5], \ 
\ \ \ \ \ \ 
            variables]; \ yyy\  = Flatten[Partition[yy\  /. \ sol, \ 1]]; \ 

        \ \ \ \ yyy2\  = \ Flatten[Partition[yy2\  /. \ sol, \ 1]]; \ 

        \ \ \ \ dyyy\  = \ Take[RotateLeft[yyy]\  - \ yyy, \ n\  - \ 1]; \ 

        \ \ \ \ dyyy2\  = Take[RotateLeft[yyy2]\  - \ yyy2, \ n\  - \ 1]; \ 

        \ \ \ \ aa\  = \ dyyy2/\((6*dx)\); 
        bb\  = \ 1/2*Take[yyy2, \ n\  - \ 1]; \ 
\ \ \ \ 
        cc\  = \ dyyy/dx\  - 
            1/6*dx*\((
                4*bb\  + \ 
\ \ \ \ \ \ \ \ \ 
                  Take[RotateLeft[yyy2, \ 1], \ n\  - \ 1])\); \ 
\ \ \ \ 
        dd\  = \ Take[yyy, \ n\  - \ 1]; \ xk\  = \ Take[x, \ n\  - \ 1]; \ 

        \ \ \ \ aa*\((z\  - \ xk)\)^3\  + \ bb*\((z\  - \ xk)\)^2\  + 
          cc*\((z\  - \ xk)\)\  + \ dd]\)

\(TraditionalForm
    \`plotsplines[x_, \ u_, \ p_]\  := \ 
\ \ 
      Module[{s1, \ t, \ x1, \ tt, \ intervals, \ tab1, \ k, \ 
\ \ \ \ n, 
          plot1, plot2}, \ s1\  = \ splinesJH[x, \ u, \ p, \ t]; \ 
\ \ \ \ 
        n\  = \ Length[x]; \ x1\  = \ RotateLeft[x]; \ 
\ \ \ \ 
        tt\  = \ Table[t, \ {n}]; \ 
\ \ \ \ 
        intervals\  = Take[Transpose[{tt, \ x, \ x1}], \ n\  - \ 1]; \ 

        \ \ \ \ tab1\  = 
          Table[Plot[s1[\([k]\)], Evaluate[intervals[\([k]\)]], \ 

              \ \ \ \ \ \ \ DisplayFunction\  -> \ Identity], {k, \ 1, \ 
              n\  - \ 1}]; \ 
\ \ \ \ 
        plot1 = Show[tab1, \ DisplayFunction\  -> $DisplayFunction]; \ 

        \ \ \ plot2 = \ ListPlot[u, \ PlotJoined\  -> \ True]; 
	
        Show[plot1, plot2]]\)