Question about how this cubic spline is solved and why?

조회 수: 11 (최근 30일)
Evan Kardos
Evan Kardos 2019년 7월 12일
댓글: Rena Berman 2021년 5월 6일
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
  댓글 수: 2
Stephen23
Stephen23 2021년 3월 27일
편집: Stephen23 2021년 3월 27일
Original question by Evan Kardos retrieved from Google Cache:
"Question about how this cubic spline is solved and why?"
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
Rena Berman
Rena Berman 2021년 5월 6일
(Answers Dev) Restored edit

댓글을 달려면 로그인하십시오.

답변 (1개)

darova
darova 2021년 3월 20일
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
  댓글 수: 1
darova
darova 2021년 3월 20일
As you can see: spline is a set of polynoms of 3d order. 4 points needed to create polynom, each polynom is used only between two points (except 1st and last, they are used between 3 points)

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Splines에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by