[Math] Cubic spline interpolation – how to calculate second derivative

algorithmsinterpolationreal-analysis

I ask this qeustion on stackexchange sites: stackoverflow, codereview, and signal processing and no one can help and they send me here 🙂

So I implement cubic spilne interpolation in Java base on this document:
http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf

At first author shows how to calculate linear spline interpolation and I did this on my data, and receive this result:

enter image description here

It should be similar to this:

enter image description here

The overall shape is good but to receive better results I should use cubic spilne intepolation with is extend of linear interpolation but here problems starts. I don't understand how they calculate second derivative.

Could you explain that to me?

I also implement cubic spilne interpolation based on algorithm from wiki:
http://en.wikipedia.org/w/index.php?title=Spline_%28mathematics%29&oldid=288288033#Algorithm_for_computing_natural_cubic_splines

x(0) is points.get(0), y(0) is values[points.get(0)], α is alfa and μ is mi Rest is the same as in pseudocode.

    public void createSpline(double[] values, ArrayList<Integer> points){
    a = new double[points.size()+1];

    for (int i=0; i <points.size();i++)
    {
        a[i] = values[points.get(i)];

    }

    b = new double[points.size()];
    d = new double[points.size()];
    h = new double[points.size()];

    for (int i=0; i<points.size()-1; i++){
        h[i] = points.get(i+1) - points.get(i);
    }

    alfa = new double[points.size()];

    for (int i=1; i <points.size()-1; i++){
        alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]);
    }

    c = new double[points.size()+1];
    l = new double[points.size()+1];
    mi = new double[points.size()+1];
    z = new double[points.size()+1];

    l[0] = 1; mi[0] = z[0] = 0;

    for (int i =1; i<points.size()-1;i++){
        l[i] = 2 * (points.get(i+1) - points.get(i-1)) - h[i-1]*mi[i-1];
        mi[i] = h[i]/l[i];
        z[i] = (alfa[i] - h[i-1]*z[i-1])/l[i];
    }

    l[points.size()] = 1;
    z[points.size()] = c[points.size()] = 0;

    for (int j=points.size()-1; j >0; j--)
    {
        c[j] = z[j] - mi[j]*c[j-1];
        b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ;
        d[j] = (c[j+1]-c[j])/((double)3*h[j]);
    }


    for (int i=0; i<points.size()-1;i++){
        for (int j = points.get(i); j<points.get(i+1);j++){
            //                fk[j] = values[points.get(i)];
            functionResult[j] = a[i] + b[i] * (j - points.get(i)) 
                                + c[i] * Math.pow((j - points.get(i)),2)
                                + d[i] * Math.pow((j - points.get(i)),3);
        }
    }

}

but result look like (That's the other part of the signal) :

enter image description here

when is should be similar to this:

enter image description here

Best Answer

According to the paper, which is a much much better reference than the outdated version of the wikipedia article, you got to solve

$$\frac{x_{k}-x_{k-1}}6 a_{k-1}+\frac{x_{k+1}-x_{k-1}}3 a_k + \frac{x_{k+1}-x_k}6 a_{k+1}=\frac{y_{k}-y_{k-1}}{x_{k}-x_{k-1}}-\frac{y_{k+1}-y_k}{x_{k+1}-x_k}$$

for $k=1,...,N-1$ and set $a_0=0=a_N$. To this end, you can implement the Jacobi or Gauss-Seidel method for linear systems.

$$a_k^{new} = \frac3{x_{k+1}-x_{k-1}}\left(\frac{y_{k}-y_{k-1}}{x_{k}-x_{k-1}}-\frac{y_{k+1}-y_k}{x_{k+1}-x_k}-\frac{x_{k}-x_{k-1}}6 a_{k-1}-\frac{x_{k+1}-x_k}6 a_{k+1}\right)$$

where you either compute all new values before replacing the current vector of $a$, or you can set $a_k=a_k^{new}$ directly after its computation.

Correction: There was a sign error on the right side, it must be the difference of the slopes; the paper has it correctly.

Remark: The system is visibly diagonal dominant, convergence should be fast.

Added: How to get a spline approximation: First variant: Use B-splines as basis functions. Second variant: Use the parametrization

$$a+bx+cx^2+\sum_{k=0}^M d_k\,(x-x_k)_+^3,\quad \text{where }(x-x_k)_+=\max(0,x-x_k)$$

with a sparse distribution of $x_0,...,x_M$ relative to the amount of sample points, and then do a standard least squares fit. If $x_M$ is outside the interval of sample points, then leave out the last term and set $d_M=-\sum_{k=0}^{M-1}d_k$ to get a proper spline with quadratic ends.

Related Question