Cubic Spline Program

Question 1

Cubic Spline Program

matlab matrix interpolation spline

Cloud Strife · Oct 4, 2011 · Viewed 34.1k times · Source

Answer

Answer

I wrote a cubic spline package in Mathematica a long time ago. Here is my translation of that package into Matlab. Note I haven't looked at cubic splines in about 7 years, so I'm basing this off my own documentation. You should check everything I say.

The basic problem is we are given n data points (x(1), y(1)) , ... , (x(n), y(n)) and we wish to calculate a piecewise cubic interpolant. The interpolant is defined as

   S(x) = {  Sk(x)   when x(k) <= x <= x(k+1)
          {  0       otherwise

Here Sk(x) is a cubic polynomial of the form

  Sk(x) = sk0 + sk1*(x-x(k)) + sk2*(x-x(k))^2 + sk3*(x-x(k))^3

The properties of the spline are:

The spline pass through the data point Sk(x(k)) = y(k)
The spline is continuous at the end-points and thus continuous everywhere in the interpolation interval Sk(x(k+1)) = Sk+1(x(k+1))
The spline has continuous first derivative Sk'(x(k+1)) = Sk+1'(x(k+1))
The spline has continuous second derivative Sk''(x(k+1)) = Sk+1''(x(k+1))

To construct a cubic spline from a set of data point we need to solve for the coefficients sk0, sk1, sk2 and sk3 for each of the n-1 cubic polynomials. That is a total of 4*(n-1) = 4*n - 4 unknowns. Property 1 supplies n constraints, and properties 2,3,4 each supply an additional n-2 constraints. Thus we have n + 3*(n-2) = 4*n - 6 constraints and 4*n - 4 unknowns. This leaves two degrees of freedom. We fix these degrees of freedom by setting the second derivative equal to zero at the start and end nodes.

Let m(k) = Sk''(x(k)) , h(k) = x(k+1) - x(k) and d(k) = (y(k+1) - y(k))/h(k). The following three-term recurrence relation holds

  h(k-1)*m(k-1) + 2*(h(k-1) + h(k))*m(k) + h(k)*m(k+1) = 6*(d(k) - d(k-1))

The m(k) are unknowns we wish to solve for. The h(k) and d(k) are defined by the input data. This three-term recurrence relation defines a tridiagonal linear system. Once the m(k) are determined the coefficients for Sk are given by

   sk0 = y(k)
   sk1 = d(k) - h(k)*(2*m(k) + m(k-1))/6
   sk2 = m(k)/2
   sk3 = m(k+1) - m(k)/(6*h(k))

Okay that is all the math you need to know to completely define the algorithm to compute a cubic spline. Here it is in Matlab:

function [s0,s1,s2,s3]=cubic_spline(x,y)

if any(size(x) ~= size(y)) || size(x,2) ~= 1
   error('inputs x and y must be column vectors of equal length');
end

n = length(x)

h = x(2:n) - x(1:n-1);
d = (y(2:n) - y(1:n-1))./h;

lower = h(1:end-1);
main  = 2*(h(1:end-1) + h(2:end));
upper = h(2:end);

T = spdiags([lower main upper], [-1 0 1], n-2, n-2);
rhs = 6*(d(2:end)-d(1:end-1));

m = T\rhs;

% Use natural boundary conditions where second derivative
% is zero at the endpoints

m = [ 0; m; 0];

s0 = y;
s1 = d - h.*(2*m(1:end-1) + m(2:end))/6;
s2 = m/2;
s3 =(m(2:end)-m(1:end-1))./(6*h);

Here is some code to plot a cubic spline:

function plot_cubic_spline(x,s0,s1,s2,s3)

n = length(x);

inner_points = 20;

for i=1:n-1
   xx = linspace(x(i),x(i+1),inner_points);
   xi = repmat(x(i),1,inner_points);
   yy = s0(i) + s1(i)*(xx-xi) + ... 
     s2(i)*(xx-xi).^2 + s3(i)*(xx - xi).^3;
   plot(xx,yy,'b')
   plot(x(i),0,'r');
end

Here is a function that constructs a cubic spline and plots in on the famous Runge function:

function cubic_driver(num_points)

runge = @(x) 1./(1+ 25*x.^2);

x = linspace(-1,1,num_points);
y = runge(x);

[s0,s1,s2,s3] = cubic_spline(x',y');

plot_points = 1000;
xx = linspace(-1,1,plot_points);
yy = runge(xx);

plot(xx,yy,'g');
hold on;
plot_cubic_spline(x,s0,s1,s2,s3);

You can see it in action by running the following at the Matlab prompt

 >> cubic_driver(5)
 >> clf 
 >> cubic_driver(10)
 >> clf
 >> cubic_driver(20)

By the time you have twenty nodes your interpolant is visually indistinguishable from the Runge function.

Some comments on the Matlab code: I don't use any for or while loops. I am able to vectorize all operations. I quickly form the sparse tridiagonal matrix with spdiags. I solve it using the backslash operator. I counting on Tim Davis's UMFPACK to handle the decomposition and forward and backward solves.

Hope that helps. The code is available as a gist on github https://gist.github.com/1269709

Question 2

I'm trying to write a cubic spline interpolation program. I have written the program but, the graph is not coming out correctly. The spline uses natural boundary conditions(second dervative at start/end node are 0). The code is in Matlab and is shown below,

clear all
%Function to Interpolate
k = 10;                    %Number of Support Nodes-1
xs(1) = -1;
for j = 1:k
    xs(j+1) = -1 +2*j/k;   %Support Nodes(Equidistant)
end;
fs = 1./(25.*xs.^2+1);     %Support Ordinates
x = [-0.99:2/(2*k):0.99];  %Places to Evaluate Function
fx = 1./(25.*x.^2+1);      %Function Evaluated at x

%Cubic Spline Code(Coefficients to Calculate 2nd Derivatives)

f(1) = 2*(xs(3)-xs(1));
g(1) = xs(3)-xs(2);
r(1) = (6/(xs(3)-xs(2)))*(fs(3)-fs(2)) + (6/(xs(2)-xs(1)))*(fs(1)-fs(2));
e(1) = 0;

for i = 2:k-2
    e(i) = xs(i+1)-xs(i);
    f(i) = 2*(xs(i+2)-xs(i));
    g(i) = xs(i+2)-xs(i+1);
    r(i) = (6/(xs(i+2)-xs(i+1)))*(fs(i+2)-fs(i+1)) + ...
           (6/(xs(i+1)-xs(i)))*(fs(i)-fs(i+1));
end
e(k-1) = xs(k)-xs(k-1);
f(k-1) = 2*(xs(k+1)-xs(k-1));
r(k-1) = (6/(xs(k+1)-xs(k)))*(fs(k+1)-fs(k)) + ...
         (6/(xs(k)-xs(k-1)))*(fs(k-1)-fs(k));

%Tridiagonal System

i = 1;
A = zeros(k-1,k-1);
while i < size(A)+1;
    A(i,i) = f(i);
    if i < size(A);
        A(i,i+1) = g(i);
        A(i+1,i) = e(i);
    end
    i = i+1;
end

for i = 2:k-1                             %Decomposition
    e(i) = e(i)/f(i-1);
    f(i) = f(i)-e(i)*g(i-1);
end

for i = 2:k-1                             %Forward Substitution 
    r(i) = r(i)-e(i)*r(i-1);
end

xn(k-1)= r(k-1)/f(k-1);
for i = k-2:-1:1                          %Back Substitution
    xn(i) = (r(i)-g(i)*xn(i+1))/f(i);
end

%Interpolation

if (max(xs) <= max(x))
    error('Outside Range'); 
end

if (min(xs) >= min(x))
    error('Outside Range'); 
end


P = zeros(size(length(x),length(x)));
i = 1;
for Counter = 1:length(x)
    for j = 1:k-1
        a(j) = x(Counter)- xs(j);
    end
    i = find(a == min(a(a>=0)));
    if i == 1
        c1 = 0;
        c2 = xn(1)/6/(xs(2)-xs(1));
        c3 = fs(1)/(xs(2)-xs(1));
        c4 = fs(2)/(xs(2)-xs(1))-xn(1)*(xs(2)-xs(1))/6;
        t1 = c1*(xs(2)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(1))^3;
        t3 = c3*(xs(2)-x(Counter));
        t4 = c4*(x(Counter)-xs(1));
        P(Counter) = t1 +t2 +t3 +t4;
    else
        if i < k-1
        c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
        c2 = xn(i+1)/6/(xs(i+1)-xs(i-1+1));
        c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
        c4 = fs(i+1)/(xs(i+1)-xs(i-1+1))-xn(i+1)*(xs(i+1)-xs(i-1+1))/6;
        t1 = c1*(xs(i+1)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(i-1+1))^3;
        t3 = c3*(xs(i+1)-x(Counter));
        t4 = c4*(x(Counter)-xs(i-1+1));
        P(Counter) = t1 +t2 +t3 +t4;
        else
        c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1));
        c2 = 0;
        c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6;
        c4 = fs(i+1)/(xs(i+1)-xs(i-1+1));
        t1 = c1*(xs(i+1)-x(Counter))^3;
        t2 = c2*(x(Counter)-xs(i-1+1))^3;
        t3 = c3*(xs(i+1)-x(Counter));
        t4 = c4*(x(Counter)-xs(i-1+1));
        P(Counter) = t1 +t2 +t3 +t4;
        end    
    end
end

P = P';
P(length(x)) = NaN;

plot(x,P,x,fx)

When I run the code, the interpolation function is not symmetric and, it doesn't converge correctly. Can anyone offer any suggestions about problems in my code? Thanks.

Cubic Spline Program

Answer

Related questions