1、三次樣條插值多項式計算物理實驗作業四物理學 2013 級程序：clear,clc;format ratx = 1,4,9,16,25,36,49,64;y = 1,2,3,4,5,6,7,8;f1 = 0.5;fn = 1/16;a,b,c,d,M,S = spline(x,y,f1,fn);子程序1：function a,b,c,d,M,S=spline(x,y,f1,fn)% 三次樣條插值函數% x 是插值節點的橫坐標% y 是插值節點的縱坐標% u 是插值點的橫坐標% f1 是左端點的一階導數% fn 是右端點的一階導數% a 是三對角矩陣對角線下邊一行% b

2、是三對角矩陣對角線% c 是三對角矩陣對角線上邊一行% S是插值點的縱坐標n = length(x);h = zeros(1,n-1);deltay = zeros(1,n);miu = zeros(1,n-1);lamda = zeros(1,n-1);d = zeros(1,n-1);for j = 1:n-1h(j) = x(j+1)-x(j);deltay(j) = y(j+1)-y(j);end % 得到 h 矩陣for j = 2:n-1sumh = h(j-1) + h(j);miu(j) = h(j-1) / sumh;lamda(j) = h(j) / sumh;d(j) =

3、 6*( deltay(j)/h(j)-(deltay(j-1)/h(j-1)/sumh;end%根據第一類邊界條件，作如下規定lamda(1) = 1;d(1) = 6*(deltay(1)/h(1)-f1)/h(1);miu(1) = 1;d(n) = 6*(fn-deltay(n-1)/h(n-1)/h(n-1);%輸出三對角矩陣的a,b,ca = miu;b = 2*ones(1,n);c = lamda;M = chase(a,b,c,d); % 調用 chase 函數得到 Msym u;for j = 1:n-1u = x(j):0.1:x(j+1);v = ones(size(u

4、);S = (M(j)*(x(j+1)*v-u)A3/(6*h(j)+M(j+1)*(u-x(j)*v)A3/(6*h(j) +(y(j)-M(j)*h(j)A2/6)*(x(j+1)*v-u)/h(j)+(y(j+1).-M(j+1)*h(j)A2/6)*(u-x(j)*v)/h(j);plot(uS'-k');hold onendplot(x,y,'-.*r');xlabel('x'),ylabel('y'),title('cubic spline interp');end子程序2:function M = c

5、hase(a,b,c,f)%追趕法求解三對角矩陣方程，Ax=f% a是對角線下邊一行的元素%b是對角線元素%c是對角線上邊一行的元素n = length(b);歡迎下載4beta = ones(1,n-1);y = ones(1,n);M = ones(1,n);for i = (n-1):(-1):1a(i+1) = a(i);end%將a矩陣和n對應beta(1) = c(1)/b(1);for i = 2:(n-1)beta(i) = c(i)/( b(i)-a(i)*beta(i-1);endy(1) = f(1)/b(1);for i = 2:ny(i) = (f(i)-a(i)*y

6、(i-1)/(b(i)-a(i)*beta(i-1); endM(n) = y(n);for i = (n-1):(-1):1M(i) = y(i)-beta(i)*M(i+1);endend三次樣條插值結果：10203。40505070Vew Inicrrt Tools Deaktep Window Help住勢賽露）目口回 與拉格朗日插值作對比：Q Figure 1-LXFi c Edit View niert oslsWindow Help，臼4 ：瑞| b W勢昵* - a 西值100豺60j 4020102。304050 GO 70分析圖一知，三次樣條插值結果與預期結果吻合得很好，曲線平滑連續性好，在左右短點處的小區間也吻合得很好，可以延伸到區間a,b外的一小段，用最鄰近的小區間插值函數可以近似求得a,b區間外