## Question

## Transcribed Text

## Solution Preview

This material may consist of step-by-step explanations on how to solve a problem or examples of proper writing, including the use of citations, references, bibliographies, and formatting. This material is made available for the sole purpose of studying and learning - misuse is strictly forbidden.

function [ynew,a,b,c,d] = NatSpline(x,y,xnew)% ynew = NatSpline(x,y)

% Natural spline interpolation.

% x,y are column n-vectors.

% Assume that n >= 4 and x(1) < ... x(n).

% xnew -- is a set of points in the interval x(1) to x(n) for which

% the natural spline needs to be evaluated.

% On the i-th interval [x(i),x(i+1)] the spline has the form:

% S(z) = a(i) + b(i)(z-x(i)) + c(i)(z-x(i))^2 + d(i)(z-x(i))^2(z-x(i+1).

%

% Outputs:

% ynew -- the values of the natural spline evaluated at the points xnew.

% a,b,c,d,-- the vector of coefficients of the splines

% Usage:

% [ynew,a,b,c,d] = NatSpline(x,y,xnew)

% First, set up all but the first and last equations that

% define the vector of interior slopes s(2:n-1).

n = length(x);

Dx = diff(x);

yp = diff(y) ./ Dx;

T = zeros(n-2,n-2);

r = zeros(n-2,1);

for i=2:n-3

T(i,i) = 2*(Dx(i) + Dx(i+1));

T(i,i-1) = Dx(i+1);

T(i,i+1) = Dx(i);

r(i) = 3*(Dx(i+1)*yp(i) + Dx(i)*yp(i+1));

end

T(1,1) = 2*Dx(1) + 1.5*Dx(2);

T(1,2) = Dx(1);

r(1) = 1.5*Dx(2)*yp(1) + 3*Dx(1)*yp(2);

T(n-2,n-2) = 1.5*Dx(n-2)+2*Dx(n-1);

T(n-2,n-3) = Dx(n-1);

r(n-2) = 3*Dx(n-1)*yp(n-2) + 1.5*Dx(n-2)*yp(n-1);

%we now need to solve

%stilde = T\r;

A=diag(T);

B=zeros(size(A));

for i=2:n-2

B(i)=T(i,i-1);

end

C=zeros(size(A));

for i=1:n-3

C(i)=T(i,i+1);

end

%call the tri diagonal solver: (below)

stilde=tridiag( A, B, C, r );

s1 = (3*yp(1) - stilde(1));

sn = (3*yp(n-1) - stilde(n-2));

s = [s1;stilde;sn];

% Compute the a,b,c,d vector of coefficients for the splines.

a = y(1:n-1);

b = s(1:n-1);

c = (yp - s(1:n-1)) ./ Dx;

d = (s(2:n) + s(1:n-1) - 2*yp) ./ (Dx.* Dx);

%next evaluate at the points xnew to obtain the values ynew

for j=1:length(xnew)

x_val=xnew(j);...