|
依据附件图片上的思想写的三次样条插值:
function yy = maspline(x,y,dy0,dyn,xx)
m = length(x);
n = length(y);
if m~=n
error('x与y长度不一致');
end
s = zeros(n+2,1);
a = zeros(n+2,4);
b = zeros(n+2,n-1);
for i =1:n
for j = 1:n-1
if i<j
b(i,j) = 0;
else
b(i,j) = (x(i)-x(j))^3/factorial(3);
end
end
a(i,1) = 1;
a(i,2) = x(i);
a(i,3) = x(i)^2;
a(i,4) = x(i)^3;
s(i) = y(i);
end
a(n+1,1) = 0;
a(n+2,1) = 0;
a(n+1,2:4) =[1 x(1) x(1)^2];
a(n+2,2:4) = [1 x(n) x(n)^2];
s(n+1) = dy0;
s(n+2) = dyn;
for i = 1:n-1
if 1>i
b(n+1,i) = (x(1)-x(i))^2/2;
end
end
for i = 1:n-1
if n>i
b(n+2,i) = (x(n)-x(i))^2/2;
end
end
A = [a b];
e = A\s;
t = ones(1,n-1);
tt = floor(xx);
ti = find(x == tt);
for i = 1:n-1
if i<=ti
t(i) = 0;
else
t(i) = (xx-x(i))^3/factorial(3);
end
end
b = [1 xx xx^2 xx^3];
B = [b t]';
yy = e'*B |
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有帐号?注册
x
|