24 views (last 30 days)

Show older comments

I'm learning about cubic splines, and coding them in matlab. We wrote code for cubic splines, and also did some work by hand, The question on my homework that's stumping me is "In the cubic spline algorithm, explain why we don’t solve the full system of equations, but instead we only solve the system for the c-coefficients and then later solve for the a, b, and d separately." I've been doing a bit of googling today but still don't really understand what the benefit of this process is? I'll include my code for the cubic splines below, if anyone could help guide me on this that would be great!!

function s = cubicInterp(x,y,z)

% x,y are arrays containing the data to interpolate

% z is a scalar or an array containing values where we would like to evaluate the interpolant

% output s is the value(s) of the interpolant evaluated at z

%convert x,y,z to column vectors

x = x(:); y = y(:); z = z(:);

%we developed cubic spline interpolation with x being length n+1

n = length(x)-1;

%preallocate space for all arrays

h = zeros(n,1);

A = zeros(n+1,n+1);

RHS = zeros(n+1,1);

b = zeros(n,1);

d = zeros(n,1);

%fill in the h array

h=diff(x);

%set A(1,1) and A(n+1,n+1)

A(1,1)=1;

A(n+1,n+1)=1;

%fill in the A matrix and the RHS vector (can use one for loop i=2:n)

for i=2:n

A(i,(i-1))=h(i-1);

A(i,i)=2*(h(i-1)+h(i));

A(i,(i+1))=h(i);

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

end

%solve the linear system for the c coefficients

RHS=RHS(:);

c=A\RHS;

%use the c coefficients to solve for the b and d coefficients

for i=1:n

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

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

end

%now that we have all the coefficients for all the splines, we can choose and evaluate the right spline

%loop over all values in the z array

for i = 1:length(z)

%for each, decide which spline to evaluate (which bin does z lie in?)

bin = binFind(x,z(i));

%use the right coefficients to evaluate the correct spline for each z value

s(i)=y(bin)+b(bin).*(z(i)-x(bin))+c(bin).*(z(i)-x(bin)).^2+d(bin).*(z(i)-x(bin)).^3;

end

end

%Modified bin code from FUNCTIONS homework

function bin = binFind(x,z0)

binArray = [x];

%x input an array

%z0 input a scalar

%bin output is the bin in x that z0 belongs to

bin = 1;

while z0 >= binArray(bin) &&bin < length(x)

bin = bin+1;

end

bin = bin-1;

end

function x = spline(varargin)

x = fprintf('This program uses the function spline, which is not allowed to be used for this homework assignment.\n');

end

darova
on 20 Mar 2021

Here is a usefull example

clc,clear

n = 6; % num of points

x = 1:n; % x range

y = rand(1,n); % random y data

pp = spline(x,y); % create coefficient

plot(x,y,'.r') % plot original data

for i = 1:size(pp.coefs,1)

x1 = linspace(x(i),x(i+1),20)-x(i); % range between points

x1 = x1*2-0.5; % larger range

y1 = polyval(pp.coefs(i,:),x1); % calculate in range

line(x1+x(i),y1)

line(x1+x(i),y1+i/5)

end

grid on

darova
on 20 Mar 2021

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!