三次樣條插值

1、三次樣條插值1. 算法原理由于在許多實際問題中，要求函數(shù)的二階導(dǎo)數(shù)連續(xù)，人們便提出了三次樣條插值函數(shù)，三次樣條插值函數(shù)是由分段三次函數(shù)拼接而成的，在連接點處二階導(dǎo)數(shù)連續(xù)。設(shè)S(x)在節(jié)點處的二階導(dǎo)數(shù)，其中為待定參數(shù)。由S（x）是分段三次多項式可知，是分段線性函數(shù)，在子區(qū)間上可以表示為其中，對S（x）兩端積分兩次得其中和為積分常數(shù)。由插值條件得由此解得代入得：求導(dǎo)得：令得在處的左導(dǎo)數(shù) 又令得在處右導(dǎo)數(shù) ，從而有，由在節(jié)點處一階導(dǎo)數(shù)的連續(xù)性知，即兩端同乘得，記，則關(guān)于的方程組寫成。三種邊界條件的三彎矩方程：（1）第一種邊界條件：已知。取，這時方程組減少了兩個未知量，變成只含n-1個未知量的n-1個

2、方程的方程組，用矩陣表示為可用追趕法求解出后，即得三條樣條插值函數(shù)。（2） 第二種邊界條件，已知，記，則有，得，即，其中，得到方程組，用矩陣表示為，該方程組的系數(shù)矩陣是嚴(yán)格三對角占優(yōu)矩陣，可用追趕法求解。（3）第三種邊界條件：周期型邊界條件。已知是以為周期的周期函數(shù)，則由周期性可知，這時將點看成是內(nèi)節(jié)點，則有，也即，其中，方程組第1個方程為：，所以方程組為，用矩陣表示為，顯然系數(shù)矩陣為嚴(yán)格對角占優(yōu)矩陣，可用LU分解法求解。2. 程序框圖3. 源程序function x=mchase(A,d)%追趕法n=length(d);u=zeros(n,1);u(1)=A(1,1);for k=2:n l

3、(k)=A(k,k-1)/u(k-1); u(k)=A(k,k)-l(k)*A(k-1,k);endy=zeros(n,1);y(1)=d(1);for i=2:n y(i)=d(i)-l(i)*y(i-1);endx=zeros(n,1);x(n)=y(n)/u(n);for i=n-1:-1:1 x(i)=(y(i)-A(i,i+1)*x(i+1)/u(i);endxendfunction T=mspline1(x0,y0,f21,f22,xx)%三次樣條插值函數(shù)第一種邊界條件%x0、y0分別為節(jié)點的橫坐標(biāo)和縱坐標(biāo)；%f21為左端點的二階導(dǎo)數(shù)值,f22為右端點的二階導(dǎo)數(shù)值；xx為由插值點組

4、成的向量n=length(x0)-1;%計算小區(qū)間數(shù)for i=1:n h(i)=x0(i+1)-x0(i);endfor i=1:n-1 mu(i)=h(i)/(h(i)+h(i+1); lamda(i)=1-mu(i); d(i)=6*(y0(i+2)-y0(i+1)/h(i+1)-(y0(i+1)-y0(i)/h(i)/(h(i)+h(i+1);endA=zeros(n-1);for i=1:n-2 A(i+1,i)=mu(i+1);%次下對角線 A(i,i+1)=lamda(i);%次上對角線 A(i,i)=2;%主對角線endA(n-1,n-1)=2;dd=zeros(n-1,1);

5、%右端列向量for i=2:n-2 dd(i)=d(i);enddd(1)=d(1)-mu(1)*f21;dd(n-1)=d(n-1)-lamda(n-1)*f22;M=mchase(A,dd);%追趕法求解M值hmulamdaAddM=f21,M',f22't=sym('t');a=zeros(n,1);b=zeros(n,1);c=zeros(n,1);e=zeros(n,1);for i=1:n a(i)=M(i)./(6*h(i); b(i)=M(i+1)./(6*h(i); W1(i)=b(i)-a(i); W2(i)=3*(a(i).*x0(i+1)

6、-b(i).*x0(i); c(i)=y0(i)./h(i)-h(i).*M(i)/6; e(i)=y0(i+1)./h(i)-h(i).*M(i+1)/6; W3(i)=3*b(i).*x0(i).2-3*a(i).*x0(i+1).2+e(i)-c(i); W4(i)=a(i).*x0(i+1).3-b(i).*x0(i).3+c(i).*x0(i+1)-e(i).*x0(i); Si(t)=W1(i).*t3+W2(i).*t2+W3(i).*t+W4(i)%每個小區(qū)間的三次樣條差值函數(shù)表達(dá)式endm=length(xx);T=zeros(m,1);for k=1:m for j=1:n

7、 if (xx(k)>=x0(j)&(xx(k)<x0(j+1) T(k)=W1(j).*(xx(k).3)+W2(j).*(xx(k).2)+W3(j).*xx(k)+W4(j); end endendTEndfunction T=mspline2(x0,y0,f11,f12,xx)%三次樣條插值函數(shù)第二種邊界條件%x0、y0分別為節(jié)點的橫坐標(biāo)和縱坐標(biāo)；%f11為左端點的二階導(dǎo)數(shù)值,f12為右端點的二階導(dǎo)數(shù)值；xx為由插值點組成的向量n=length(x0)-1;%計算小區(qū)間數(shù)for i=1:n h(i)=x0(i+1)-x0(i);endfor i=1:n-1 mu(i

8、)=h(i)/(h(i)+h(i+1); lamda(i)=1-mu(i); d(i)=6*(y0(i+2)-y0(i+1)/h(i+1)-(y0(i+1)-y0(i)/h(i)/(h(i)+h(i+1);endA=zeros(n+1);%系數(shù)矩陣dd=zeros(n+1,1);%右端列向量for k=2:n A(k,k)=2;%主對角線元素 A(k,k-1)=mu(k-1);%次下對角線元素 A(k,k+1)=lamda(k-1);%次上對角線元素endA(1,1)=2;A(1,2)=1;A(n+1,n+1)=2;A(n+1,n)=1;dd(1)=6*(y0(2)-y0(1)/h(1)-f1

9、1)/h(1);dd(n+1)=6*(f12-(y0(n+1)-y0(n)/h(n)/h(n);for k=2:n dd(k)=d(k-1);endM=mchase(A,dd);%追趕法求解M值hmulamdaAddMt=sym('t');a=zeros(n,1);b=zeros(n,1);c=zeros(n,1);e=zeros(n,1);for i=1:n a(i)=M(i)./(6*h(i); b(i)=M(i+1)./(6*h(i); W1(i)=b(i)-a(i); W2(i)=3*(a(i).*x0(i+1)-b(i).*x0(i); c(i)=y0(i)./h(i

10、)-h(i).*M(i)/6; e(i)=y0(i+1)./h(i)-h(i).*M(i+1)/6; W3(i)=3*b(i).*x0(i).2-3*a(i).*x0(i+1).2+e(i)-c(i); W4(i)=a(i).*x0(i+1).3-b(i).*x0(i).3+c(i).*x0(i+1)-e(i).*x0(i); Si(t)=W1(i).*t3+W2(i).*t2+W3(i).*t+W4(i)%每個小區(qū)間的三次樣條差值函數(shù)表達(dá)式endm=length(xx);T=zeros(m,1);for k=1:m for j=1:n if (xx(k)>=x0(j)&(xx(

11、k)<x0(j+1) T(k)=W1(j).*(xx(k).3)+W2(j).*(xx(k).2)+W3(j).*xx(k)+W4(j); end endendTend%計算實習(xí)，第一種邊界條件clear;x=-1:0.2:1;%輸入節(jié)點橫坐標(biāo)y=ones(1,11)./(ones(1,11)+25*x.2)%計算節(jié)點縱坐標(biāo)t=sym('t') ;f=1/(1+25*t2);%定義函數(shù)f2=diff(f,2)%函數(shù)的二階導(dǎo)數(shù)式f21=vpa(subs(f2,'t',-1)%計算左端點的二階導(dǎo)數(shù)值f22=vpa(subs(f2,'t',1)%

12、計算右端點的二階導(dǎo)數(shù)值xx=-1:0.1:1;T=mspline1(x,y,f21,f22,xx); T=T'ezplot(f,-1 1);%畫出函數(shù)f的曲線hold onplot(x,y,':',xx,T,'-');%根據(jù)函數(shù)計算和插值計算的結(jié)果畫出的曲線%計算實習(xí)，第二種邊界條件clear;x=-1:0.2:1;%輸入節(jié)點橫坐標(biāo)y=ones(1,11)./(ones(1,11)+25*x.2)%計算節(jié)點縱坐標(biāo)t=sym('t') ;f=1/(1+25*t2);%定義函數(shù)f1=diff(f,1)%函數(shù)的一階導(dǎo)數(shù)式f11=vpa(subs(f1,'t',-1)%計算左端點的一階導(dǎo)數(shù)值f12=vpa(subs(f1,'t',1)%計算右端點的一階導(dǎo)數(shù)值xx=-1:0.1:1;T=mspline2(x,y,f11,f12,xx); T=T'ezplot(f,-1 1);%畫出函數(shù)f的曲線hold onplot(x,y,':',xx,T,'-');%根據(jù)函數(shù)計算和插值計算的結(jié)果畫出的曲線4. 計算結(jié)果第一種邊界條件：x-1-0.8-0.6-0.4-0.200.20.40.60.81函數(shù)計算值0.03850.05880.