Question

Modify this Matlab code for natural spline so it computes clamp Spline compclampspline(x,y,alpha,beta)

alpha and beta are the derivatives at the ends.

Code for Natural Spline

function [a,b,c,d]=compnatspline(x,y)

function [a,b,c,d]=compnatspline(x,y)

n=length(x);
delx=(x(2:n)-x(1:n-1))';
dely=(y(2:n)-y(1:n-1))';
a=zeros(n,n);
a(1,1)=1; a(n,n)=1;
r(1)=0; r(n)=0;
for row=2:n-1
r(row)=3*(dely(row)/delx(row)-dely(row-1)/delx(row-1));
a(row,row-1)=delx(row-1);
a(row,row)=2*(delx(row-1)+delx(row));
a(row,row+1)=delx(row);
end
c=a\r';
clear a
b(1:n-1)=dely./delx -delx/3.*(2*c(1:n-1)+c(2:n));
d(1:n-1)=(c(2:n)-c(1:n-1))./(3*delx(1:n-1));
a(1:n-1)=y(1:n-1);
c=c(1:n-1);

Answer 1

The modified Matlab script compClampSpline.m is posted below:

function [a, b, c, d] = compClampSpline(x, y, alpha, beta)

n = length(x);
delx = ( x(2:n)-x(1:n-1) )';
dely = ( y(2:n)-y(1:n-1) )';

a = zeros(n,n);
a(1, 1:2) = [2*delx(1) delx(1)]; r(1) = 3*(dely(1)/delx(1) - alpha);
a(n, n-1:n) = [delx(n-1) 2*delx(n-1)]; r(n) = 3*(beta - dely(n-1)/delx(n-1));

for row = 2:n-1
r(row) = 3*(dely(row)/delx(row)-dely(row-1)/delx(row-1));
a(row, row-1) = delx(row-1);
a(row, row) = 2*(delx(row-1) + delx(row));
a(row, row+1) = delx(row);
end

c = a\r';
clear a

b(1:n-1) = dely./delx -delx/3.*(2*c(1:n-1)+c(2:n));
d(1:n-1) = (c(2:n)-c(1:n-1))./(3*delx(1:n-1));
a(1:n-1) = y(1:n-1);
c = c(1:n-1);