Natural Spline Interpolation Matlab Coefficients

I made matlab code to find the natural cubic spline. The question wants me to evaluate a natural cubic spline at different S(x) values. i am able to do that and get correct responses but the question also asks for the aj,bj,cj,dj,xj (that are in the code) at the current S(x) value and i can not figure out how to find those values at the current S(x) value. Could anyone help me figure this out? The question is asking for

S(x), x=1.22 S(x),x=5.5 S(x)=2.2

This is the initial array

xi=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 ... 
    9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3]; 
a=[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 .25];

This is the code i have :

% Cubic Spline Interpolation - Natural Spline 
% INPUT: xi is the grid (points on x-axis) and a are points on y-axis. inter
% is the point on the x-axis you want to know the value of on the y-axis.

function [a_inter] = cubic_spline(xi,a,inter)

if length(xi) ~= length(a) 
  error('vectors xi and a must be of same length'); 
end

% Plotting points we want to interpolate between: 
grid on; 
hold on;     
title('Cubic Spline Interpolation'); 
plot(xi,a,'or');

n = length(xi);

% Vector h with subintervals: 
h = zeros(n-1,1); 
for j = 1:n-1
  h(j) = xi(j+1) - xi(j); 
end

% Coefficient matrix A: 
A = zeros(n);

% Natural Spline boundary conditions: 
A(1,1)= 1; 
A(n,n) = 1;

for i = 2:n-1 
  A(i,i-1) = h(i-1); 
  A(i,i) = 2*(h(i-1)+h(i)); 
  A(i,i+1) = h(i); 
end

% Vector b: 
b = zeros(n,1);
[it wasn't clear to me if this was supposed to be commented out or not..]:
for i = 2:n-1 
  b(i) = (3/h(i))(a(i+1)-a(i)) - (3/h(i-1))(a(i)-a(i-1));     
end

% Coefficient vector cj: 
cj = A\b;

% Coefficient vector bj: 
bj = zeros(n-1,1); 
for i = 1:n-1 
  bj(i) = (1/h(i))(a(i+1)-a(i)) - (1/3h(i))(2cj(i)+cj(i+1)); 
end

% Coefficient vector dj: 
dj = zeros(n-1,1); 
for i = 1:n-1 
  dj(i) = (1/(3*h(i))) * (cj(i+1)-cj(i)); 
end

% Making a matrix P with all polynomials 
P = zeros(n-1,4); 
for i = 1:n-1 
  P(i,1) = dj(i); 
  P(i,2) = cj(i); 
  P(i,3) = bj(i); 
  P(i,4) = a(i); 
end

% Plotting results: 
resolution = 100;

for i = 1:n-1 
  f = @(x) a(i) + bj(i).(x-xi(i)) + cj(i).(x-xi(i)).2 + dj(i).*(x-xi(i)).3; 
  xf = linspace(xi(i),xi(i+1),resolution); 
  plot(xf,f(xf),'b'); 
end

% Given a value on the x-axis, inter, that we wish to know the y-value of, 
% we must first find in which interval inter is. We will use bisection 
% search to accomplish this. Interval must be ascending or descending.

jl = 1; 
ju = n; 
while (ju - jl > 1) 
  jm = ceil((jl + ju)/2); 
  if (inter < xi(jm)) 
    ju = jm; 
  else 
    jl = jm; 
  end 
end

a_inter = polyval(P(jl,:), inter-xi(jl)); 
fprintf('\n The interpolated value is: %f \n', a_inter); 
plot(inter, a_inter, 'og'); 
fprintf('The value of bj is %f \n',bj); 
fprintf('The value of cj is %f \n',cj); 
fprintf('The value of dj is %f \n',dj);
end % [end of function]

[Image Analyst]

Read this: http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup and attach your m-files with the paper clip icon.

[John D'Errico]

Ye gawds! I can read it now after doing some repairs!

[John D'Errico]

So this clearly is not your code.

1. It is more well written than a student could have done.
2. The comments are not what a student would have written.
3. You don't seem to understand the code you have in your hands.

To that I would add that the code is not what I personally would call good, but more what an instructor might write who had some familiarity with splines and MATLAB, but not a huge amount of either. For example, use of bisection to determine an interval is silly if you have histc available in MATLAB. Histc will be faster & vectorized, so any number of points could theoretically be evaluated. Teach students to use the tools in MATLAB that will make their code efficient. As well, matrices are constructed using loops, where a simple one line assignment would suffice.

The polynomial coefficients are constructed in a somewhat kludgy way, although we are never really told what those coefficients represent. What do bj, cj, and dj mean in context? Code is useless unless it explains itself, and we lack your class notes, where presumably the instructor formulated a natural spline in terms of some coefficients. Vectors are just a list of arbitrary numbers unless some hint if given as to what they mean. And since there are MANY ways to parameterize a cubic spline, I'm not going to bother to read through this code that deeply to decide which one was used here.

At the same time, backslash is used to solve a linear system, so someone did know at least some of the basics.

Answer 1

[Mvrht Minute]

code from reddit