Spline Coefficients

Spline Coefficients

Getting the Spline Coefficients

  1. Set up the function file (inputs, outputs)
    2. n = length(x)
    3. Preallocate space for h (length n – 1), A (n by n), RHS (length n), b(length n – 1), and d (length n – 1).
    4. Fill in an array for h
    5. Set A(1; 1) and A(n; n)
    6. Fill in the matrix A and the right hand side vector RHS. You can do itwith just one for loop from i=2:n-1.
    7. Backslash solve. Put the solution in an array c.
    8. Use the c values to fill in arrays for b and d. If you used b as the variablefor the right hand side in the linear system, you will be overwriting thosevalues here.

Evaluating the Interpolant at z

At this point, we have all the coefficients for all the splines. To evaluate
the interpolant, we need to determine with spline to evaluate (which set of
coefficients to use). We need to know which interval z lies in. For example, if
x3≤ z ≤ x4, we evaluate S3(x) at z:
S3(z) = a3+ b3(z – x3) + c3(z – x3)2+ d3(z – x3)3
For one scalar z, we can simply use a while loop to determine which spline to
evaluate. If z is an array of values, we can put that while loop inside a for loop,
and loop over all values in z.
Or, the following code takes advantage of some tricks in Matlabwithout
while and for loops:
%preallocate k to be the same length as z
k = ones(size(z));
%fill in k so that, for example, k(5) equals the index
% of the spline that should be evaluated to get S(z(5))
for i = 2:n-1
k(x(i) <= z) = i;
end
%now, just evaluate the splines at the z values
%with the dot operator, we can evaluate at every z value at once
s = y(k) + b(k).*(z-x(k)) + c(k).*(z-x(k)).ˆ2 + d(k).*(z-x(k)).ˆ3;
You will need to make sure that your coefficients arrays a, b, c, and d are
the same size as your z array (that they are all columns or all rows).

Solution 

cubicinerp.m 

function s=cubicinerp(x,y,z)

n=length(x);

% Space h with subintervals:

h = zeros(n-1,1);

for i = 1:n-1

h(i) = x(i+1) – x(i);

end

% Coefficient matrix A (n by n)

A = zeros(n);

% A matrix 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);

for i = 2:n-1

b(i) = (3/h(i))*(y(i+1)-y(i))-(3/h(i-1))*(y(i)-y(i-1));

end

% Coefficient vector c:

c=A\b;

%c=c(1:end-1);

% Coefficient vector d:

d = zeros(n-1,1);

for i = 1:n-1

d(i) = (1/(3*h(i)))*(c(i+1)-c(i));

end

% Coefficient vector b:

bc = zeros(n-1,1);

for i = 1:n-1

bc(i) = (1/h(i))*(y(i+1)-y(i))-(1/3*h(i))*(2*c(i)+c(i+1));

end

%  matrix P with all polynomials

P = zeros(n-1,4);

for i = 1:n-1

P(i,1) = d(i);

P(i,2) = c(i);

P(i,3) = bc(i);

P(i,4) = y(i);

end

k = ones(size(z));

n=length(x);

%fill in k so that, for example, k(5) equals the index

% of the spline that should be evaluated to get S(z(5))

for i = 2:n-1

k(x(i) <= z) = i;

end

for i=1:length(z)

s(i) = y(k(i))+bc(k(i))*(z(i)-x(k(i)))+ c(k(i))*(z(i)-x(k(i)))^2 + d(k(i))*(z(i)-x(k(i)))^3;

end

end 

main1.m 

clear all;

clc

x=[0.06 0.47 1.01 1.50 2.05 2.53 2.99 3.44];

y=[1.0499 1.3274 1.3588 1.2550 1.6322 2.5523 3.4462 3.6684];

z=0:0.01:4;

s=cubicinerp(x,y,z);

figure(1)

plot(x,y,’*r’,z,s,’-b’)

xlabel(‘x’)

ylabel(‘y’)

legend(‘real values’,’interpolation result’)