该楼层疑似违规已被系统折叠隐藏此楼查看此楼
function
S=Threch1(X,Y,dy0,dyn,xi)
% X
为已知数据的横坐标
%Y
为已知数据的纵坐标
%xi
插值点处的横坐标
%S
求得的三次样条插值函数的值
%dy0
左端点处的一阶导数
% dyn
右端点处的一阶导数
n=length(X)-1;
d=zeros(n+1,1);
h=zeros(1,n-1);
f1=zeros(1,n-1);
f2=zeros(1,n-2);
for
i=1:n
%
求函数的一阶差商
h(i)=X(i+1)-X(i);
f1(i)=(Y(i+1)-Y(i))/h(i);
end
for
i=2:n
%
求函数的二阶差商
f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));
d(i)=6*f2(i);
end
d(1)=6*(f1(1)-dy0)/h(1);
d(n+1)=6*(dyn-f1(n-1))/h(n-1);
%¸
赋初值
A=zeros(n+1,n+1);
B=zeros(1,n-1);
C=zeros(1,n-1);
for
i=1:n-1
B(i)=h(i)/(h(i)+h(i+1));
C(i)=1-B(i);
end
A(1,2)=1;
A(n+1,n)=1;
for
i=1:n+1
A(i,i)=2;
end
for
i=2:n
A(i,i-1)=B(i-1);
A(i,i+1)=C(i-1);
end
M=A\d;
syms
x
;
for
i=1:n
Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))
...
+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);
digits(4);
Sx(i)=vpa(Sx(i));%
三样条插值函数表达式
end
for
i=1:n
disp(
'S(x)='
);
fprintf(
'%s (%d,%d)\n'
,char(Sx(i)),X(i),X(i+1));
end
for
i=1:n
if
xi>=X(i)&&xi<=X(i+1)
S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M
(i+1)-M(i))/(6*h(i))*(xi-X(i))^3;
end
end
disp(
'xi S'
);
fprintf(
'%d,%d\n'
,xi,S);
return
