PAGE
PAGE #/ 4
第一邊界條件)源代碼：
function y=yt1(x0,y0,f_0,f_n,x)(1)
%第一類邊界條件下三次樣條插值；
%xi 所求點；
%yi 所求點函數值；
%x已知插值點；
%y已知插值點函數值；
%f_0 左端點一次導數值；
%f_n 右端點一次導數值；
n = length(x0);
z = length(y0);
h = zeros(n-1,1);
k=zeros(n-2,1);
l=zeros(n-2,1);
S=2*eye(n);
fori=1:n-1
h(i)= x0(i+1)-x0(i);
end
fori=1:n-2
k(i)= h(i+1)/(h(i+1)+h(i));
l(i)= 1-k(i);
end %對于第一種邊界條件：
k = [1;k];(2)
l = [l;1];(3)
%建構系數矩陣 S：
fori = 1:n-1
S(i,i+1) = k(i);
S(i+1,i) = l(i);
end
%建立均差表：
F=zeros(n-1,2);
fori = 1:n-1
F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i)); end
D = zeros(n-2,1);
fori = 1:n-2
F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i));
D(i,1) = 6 * F(i,2);
end
%建構函數 D：
d0 = 6*(F(1,2)-f_0)/h(1);(4)
dn = 6*(f_n-F(n-1,2))/h(n-1);(5)
D = [d0;D;dn];(6)
m= S\D;
%尋找 x 所在位置，并求出對應插值：
fori = 1:length(x)
for j = 1:n-1
if (x(i)<=x0(j+1))&(x(i)>=x0(j))
y(i) =( m(j)*(xO(j+1)-x(i))八3)/(6*h(j))+…
(m(j+1)*(x(i)-xO(j))八3)/(6*h(j))+… (yO(j)-(m(j)*h(jF2)/6)*(xO(j+1)-x(i))/h(j)+… (y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j); break;
else continue; end end end
( 2 )(自然邊界條件)源代碼： 僅僅需要對上面部分标注的位置做如下修改
__(1):function y=yt2(x0,y0,x)
__(2):k=[0;k]
__(3):l=[l;0]
__(4)+(5):删除
—(6):D=[0:D:0]