Matlab Natural Spline interpolation find coefficients

Question

zach on 24 Oct 2014

0
Link

Direct link to this question

https://au.mathworks.com/matlabcentral/answers/159899-matlab-natural-spline-interpolation-find-coefficients

Commented: Image Analyst on 24 Oct 2014

I have a problem for a natural cubic spline.

I was given these three matlab (evalspline.m, trisolve.m, and spline skeleton) files to fill in and get the proper responses:

%evalspline.m % function S = evalspline(a,b,c,d,x,y) % % Inputs: % a-b-c-d the coefficients that define the spline % x are the node-points % y the point(s) where you want to evaluate the spline %

function S = evalspline(a,b,c,d,x,y)

if( length(a) ~= length(x) ) error(['The first argument (a) must be the same length as the 5th' ... ' argument (x)']); end

if( length(b) < length(x)-1 ) error(['The second argument (b) cannot be more then one element' ... ' shorter than the 5th argument (x)']); end

if( length(c) ~= length(x) ) error(['The third argument (c) must be the same length as the 5th' ... ' argument (x)']); end

if( length(d) < length(x)-1 ) error(['The second argument (d) cannot be more then one element' ... ' shorter than the 5th argument (x)']); end

if( min(y) < min(x) ) error(['The values in the 6th argument (y) should not be less than' ... ' the values in the 5th argument (x)']); end

if( max(y) > max(x) ) error(['The values in the 6th argument (y) should not be greather' ... ' than the values in the 5th argument (x)']); end

S = zeros(size(y));

for k = 1:length(y) n = min([length(b) max(find(x<=y(k)))]); S(k) = polyval([d(n) c(n) b(n) a(n)],(y(k)-x(n))); end

%Trisolve.m % Solver for tri-diagonal matrices. % % T must be in compact form: the sub-diagonal elements in the first % column, the diagonal in the second column, and the super-diagonal % in the third column. % % Note that T(1,1) and T(n,3) are never accessed, i.e. the % subdiagonal entries start on the second row, and the % superdiagonal elements end on the (n-1)st row. % function [x]=trisolve(T,b)

if nargin ~= 2 error('trisolve:: two input arguments required.') end if nargout ~= 1 error('trisolve: one output argument required.') end

[n,m] = size(T); if m ~= 3 error('trisolve: matrix must have three columns.') end

[N,M] = size(b);

if N ~= n error('trisolve: rhs and matrix must have same number of rows.'); end

work = zeros(n,1);

work(1) = T(1,2); x(1,:) = b(1,:);

% Forward sweep. for i=2:n beta = T(i,1)/work(i-1); x(i,:) = b(i,:) - beta*x(i-1,:); work(i) = T(i,2) - beta*T(i-1,3); end

x(n,:) = x(n,:)/work(n);

% Backward sweep. for i=n-1:-1:1 x(i,:) = (x(i,:) - T(i,3)*x(i+1,:)) / work(i); end

%%% %%% It gives you the structure of the program, but you have %%% to fill in the missing details!!! %%%

clear format long load duckdata f = f'; % reset f as a column vector x = x'; % reset x as a column vector n = length(x); % set n equal to the number of x-values a = f; % set a equal to the column vector 'f' h = diff(x);

% Matrix 'T' is the compact form of an nxn matrix 'A' where T(:,1) is % the sub-diag of A, T(:,2) is the diag of A, and T(:,3) is the % super-diag of A % % Here "A" is the matrix corresponding the the linear system for % getting the c-coefficients for a natural spline % % 'rhs' is the right-hand-side...

%%%(YOUR CODE HERE)

%%%(END YOUR CODE, PART 1)

c = trisolve(T,rhs); % function that solves for c where T*c=rhs

% % Next, build vectors 'b' and 'd' %

%%%(YOUR CODE HERE)

%%%(END YOUR CODE, PART 2)

y = min(x):0.01:max(x); % set y to arbitary pts to evaluate % the spline at y = y'; % reset y as a column vector

S = evalspline(a,b,c,d,x,y); % function that uses the coeff % defining the spline to evaluate the % spline at each pt in 'y'

% graph pts (x,f) and the spline that fits these pts

plot(y,S,'b-',x,f,'mo','markersize',7) grid on disp('This should look like the duck in the book...') disp('...except maybe for some scaling of the axes.') disp(' ')

%%% Make sure these are column vectors... a = reshape(a,length(a),1); b = reshape(b,length(b),1); c = reshape(c,length(c),1); d = reshape(d,length(d),1);

disp('Coefficients [a b c d]') disp([a [b;0] c [d;0]])

%%%(NOTE) In order to answer all questions, you may have to dig %%%(NOTE) through the code/results carefully to get all the %%%(NOTE) information you need.

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]

For S(x)=1.22 i need to find xk,ck,dk,bk,ak

I was able to figure out S(1.22)=1.4646 but i could not retrieve the other values.

For S(x)=5.5 i need to find xl,cl,dl,bl,al

I was able to figure out S(5.5)=2.1977 but i could not retrieve the other values.

For S(x)=2.25 i need to find xm, cm,dm,bm,am

I was able to figure out S(2.25)=2.1977 but i could not retrieve the other values.