三次樣條插值多項式matlab

1、三次樣條插值多項式計算物理實驗作業(yè)四陳萬 物理學(xué)2013級主程序：clear,clc;format ratx = 1,4,9,16,25,36,49,64;y = 1,2,3,4,5,6,7,8;fl =;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)%三次樣條插值函數(shù)% x是插值節(jié)點的橫坐標% y是插值節(jié)點的縱坐標% u是插值點的橫坐標% f1是左端點的一階導(dǎo)數(shù)% fn是右端點的一階導(dǎo)數(shù)%a是三對角矩陣對角線下邊一行% b是三對角矩陣對角線% c是三對角矩陣對角線上邊一

2、行% 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) = 6*( deltay(j)/h(j)-(delt

3、ay(j-1)/h(j-1)/sumh;end% 根據(jù)第一類邊界條件，作如下規(guī)定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);%I用 chase 函數(shù)得到 Msym u;for j = 1:n-1u = x(j):x(j+1);v = ones(size(u);S = (M(j)*(x(j+1)*v-u).A3/(

4、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(u,S,'-k');hold onendplot(x,y,'-.*r');xlabel('x'),ylabel('y'),title('cubic spline interp');end子程序2:function M = chase(a,b,c,f)%追趕法求解三對角矩陣

5、方程，Ax=f% a是對角線下邊一行的元素%b是對角線元素%c是對角線上邊一行的元素n = length(b);beta = 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對應(yīng)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(i-1)/(b(i)-a(i)*beta(i-1);en

6、dM(n) = y(n)；for i = (n-1):(-1):1M(i) = y(i)-beta(i)*M(i+1);end end三次樣條插值結(jié)果： Figure 11 Xfife Ed；t View Insert 工 gk叩 我 indw Help心巴 u -1 琳 71國正/, af zi a " ncubic sfpiiTie interp10203040505070與拉格朗日插值作對比:Q Figure 1-OXEditInipd Tools Delktcp Window HelpJtJld必I信 工，有識/ - 以 回值10 EO 304050 GO 即分析圖一知，三次樣條插值結(jié)果與預(yù)期結(jié)果吻合得很好，曲線平滑連續(xù)性好， 在左右短點處的小區(qū)間也吻合得很好，可以延伸到區(qū)間a,b外的一小段，用最鄰近的小區(qū)間插信函數(shù)可以近似求得a,b區(qū)