How to implement a piecewise function and then plo

Question

I am actually attempting to write code for the cubic spline interpolation. Cubic spline boils down to a series of n-1 segments where n is the number of original coordinates given initially and the segments are each represented by some cubic function.

I have figured out how to get all the coefficients and values for each segment, stored in vectors a,b,c,d, but I don't know how to plot the function as a piecewise function on different intervals. Here is my code so far. The very last for loop is where I have attempted to plot each segment.

%initializations
x = [1 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6 7 8 9.2 10.5 11.3 11.6 12 12.6 13 13.3].';
y = [1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25].';

%n is the amount of coordinates
n = length(x);
%solving for a-d for all n-1 segments
a = zeros(n,1);
b = zeros(n,1);
d = zeros(n,1);

%%%%%%%%%%%%%% SOLVE FOR a's %%%%%%%%%%%%%
%Condition (b) in Definition 3.10 on pg 146
%Sj(xj) = f(xj) aka yj
for j = 1: n
    a(j) = y(j);
end

%initialize hj
h = zeros(n-1,1);

for j = 1: n-1
    h(j) = x(j+1) - x(j);
end

A = zeros(n,n);
bv = zeros(n,1); %bv = b vector

%initialize corners to 1
A(1,1) = 1;
A(n,n) = 1;

%set main diagonal
for k = 2: n-1
    A(k,k) = 2*(h(k-1) + h(k));
end

%set upper and then lower diagonals
for k = 2 : n-1
    A(k,k+1) = h(k); %h2, h3, h4...hn-1
    A(k,k-1) = h(k-1); %h1, h2, h3...hn
end

%fill up the b vector using equation in notes
%first and last spots are 0
for j = 2 : n-1
    bv(j) = 3*(((a(j+1)-a(j)) / h(j)) - ((a(j) - a(j-1)) / h(j-1)));
end

%augmented matrix
A = [A bv];



%%%%%%%%%%%% BEGIN GAUSSIAN ELIMINATION %%%%%%%%%%%%%%%
offset = 1;
%will only need n-1 iterations since "first" pivot row is unchanged
for k = 1: n-1
  %Searching from row p to row n for non-zero pivot
  for p = k : n
      if A(p,k) ~= 0;
          break;
      end
  end

  %row swapping using temp variable
  if p ~= k
      temp = A(p,:);
      A(p,:) = A(k,:);
      A(k,:) = temp;
  end


  %Eliminations to create Upper Triangular Form
  for j = k+1:n
     A(j,offset:n+1) = A(j,offset:n+1) - ((A(k, offset:n+1) * A(j,k)) / A(k,k));
  end
  offset = offset + 1;
end

c = zeros(n,1); %initializes vector of data of n rows, 1 column

%Backward Subsitution
%First, solve the nth equation
c(n) = A(n,n+1) / A(n,n);

%%%%%%%%%%%%%%%%% SOLVE FOR C's %%%%%%%%%%%%%%%%%%
%now solve the n-1 : 1 equations (the rest of them going backwards
for j = n-1:-1:1 %-1 means decrement
    c(j) = A(j,n+1);
    for k = j+1:n
        c(j) = c(j) - A(j,k)*c(k); 
    end
    c(j) = c(j)/A(j,j);
end

%%%%%%%%%%%%% SOLVE FOR B's and D's %%%%%%%%%%%%%%%%%%%%
for j = n-1 : -1 : 1
    b(j) = ((a(j+1)-a(j)) / h(j)) - (h(j)*(2*c(j) + c(j+1)) / 3);
    d(j) = (c(j+1) - c(j)) / 3*h(j);
end

%series of equation segments

for j = 1 : n-1
    f = @(x) a(j) + b(j)*(x-x(j)) + c(j)*(x-x(j))^2 + d(j)*(x-x(j))^3;
end
plot(x,y,'o');

Let's assume that I have calculated vectors a,b,c,d correctly for each segment. How do I plot each cubic segment such that they all appear graphed on a single plot?

I appreciate the help.

Answer 1

That's pretty easy. You've already done half of the work by defining an anonymous function that is for the cubic spline in between each interval. However, you need to make sure that the operations in the function are element-wise. You currently have it operating on scalars, or assuming that we are using matrix operations. Don't do that. Use .* instead of * and .^ instead of ^. The reason why you need to do this is to make generating the points on the spline a lot easier, where my next point follows.

All you have to do next is define a bunch of x points within the interval defined by the neighbouring x key points and substitute them into your function, then plot the result.... so something like this:

figure;
hold on;
for j = 1 : n-1
    f = @(x) a(j) + b(j).*(x-x(j)) + c(j).*(x-x(j)).^2 + d(j)*(x-x(j)).^3; %// Change function to element-wise operations - be careful
    x0 = linspace(x(j), x(j+1)); %// Define set of points
    y0 = f(x0); %// Find output points
    plot(x0, y0, 'r'); %// Plot the line in between the key points
end
plot(x, y, 'bo');

We spawn a new figure, then use hold on so that when we call plot multiple times, we append the results to the same figure. Next, for each set of cubic spline coefficients we have, define a spline function, then generate a bunch of x values with linspace that are between the current x key point and the one beside it. By default, linspace generates 100 points between a start point (i.e. x(j)) and end point (i.e. x(j+1)). You can control how many points you want to generate by specifying a third parameter (so something like linspace(x(j), x(j+1), 25); to generate 25 points). We use these x values and substitute them into our spline equation to get our y values. We then plot this result on the figure using a red line. Once we're done, we plot the key points as blue open circles on top of the curve.

As a bonus, I ran your code with the above plotting mechanism, and this is what I get: