三次样条插值

《三次样条插值》由会员分享，可在线阅读，更多相关《三次样条插值（9页珍藏版）》请在装配图网上搜索。

1、三次样条插值1. 算法原理由于在许多实际问题中，要求函数的二阶导数连续，人们便提出了三次样条插值函数，三次样条插 值函数是由分段三次函数拼接而成的，在连接点处二阶导数连续。设S(x)在节点Xi处的二阶导数S”(叮二MN二丄，n)，其中Mi为待定参数。由S (x)是分 段三次多项式可知，S(X)是分段线性函数，S(X)在子区间J上可以表示为x xx xx xM =t M +x x i hi1ii 1iS(x) = x_xi M +x xi1ii1M , x x x hii1ii其中hi = xi 一备，对S(x)两端积分两次得(x x)+c (xx ),x x x iii 1i 1i(x x)3

2、(x x )S (x) = i M +t-1 M + b6hi 16hi iii 1bc其中bi和i为积分常数。由插值条件S(x )= y ., S(x )二 y 得i-1i-1h2+ M + b h6i 1iih2=y，i-M + c hi 16 i i i6hix)3(x x )3M +iM +i 16hiiyi-1h m6 T-1 丿xxt1, x x xhi 1ii由此解得rh 2rh 2)b=y - M/ h ,c =iT-16i-1ii、i 61 ,代入得：求导得： x )2(x xM + 2hi1i令x = x得S(x)在x处的左导数iiS (x )=i-12hiM Mii1 h

3、 , x6T i x xi 1i丄h丄y yi-1+ t M +i-13 i hi又令x = Ji得SC)在hi处右导数S (x )M + i12 i1hiM M ii h6 1=Zm3 i-1h “ 丄 y - y J M +i，6 i hihhy -y从 而 有 S (X) i+1 M i+1 M+ i+1i+3 i 6i +1hi+1，由S (x)在节点x 处一阶导数的连续性 知 iS (x ) S (x ),i 1,2,n 1，即 - i+ ih皿丄h + h 皿丄h y - y y - yd+ ii+1+ i+1i+1i ii1，6i-13 i 6i+1hhi+1ii 1,2, .，

4、两端同乘h+hii6h得7 h +hiii+1hM+ 2 M +i+ Mi -1i h + hi +1i +1ii +1ii+1一 yhi+1y yii1hihhH i,九 屮1 Hi h + h i h + hiii+1ii+16h + hii+1(y - yi+1il hi+16fxi-1,xi,xi+1,i1,2,.,n-1则关于 M 的方程组写成iH M + 2M + XM d , i 1,2,n 1。ii -1ii+1i三种边界条件的三弯矩方程：第一种边界条件：已知广S广(b )。取气-广心)，叫广(b)，这时方程组减少了两个未知量，变成只含n-1个未知量M ( 1,2,., n-1

5、)的n-1个方程的方程组i2 M + 九 M d-H Mi i 2i 0H M+ 2 M + 九M d , i 2,3,，ii -1ii +1iH M + 2 M d 九 Mn-1n-2n-2n 2 ，用矩阵表示为n-1n-1n -1nHn-2n -1即得三条样条插值函数。o人M2)第 二 种 边 界 条 件 ，已 知 f (a)，f (b) ， 记y0 f(a),yn f(b) ， 则 有九n-22可用追赶法求解出M (i二1,2,., n-1)后,iS+(x )y ,S-(x )y ，得+00- n nhy -y+ n M + n-1 y ，n -13 n hnnhhy -yh-1M -M

6、 + 厶 yo y ,fM3061 h 0 612M + M = d , M+ 2M = d ， 其 中 d = 6 fx , x , x , d010n -1n n00nn12二6 f x, x , x ,得到方程组n -1 n n2 M + M = d0 1 0=d , i = 1,2,i=x0(j)&(xx(k)=x0(j)&(xx(k)x0(j+1) T(k)=Wl(j).*(xx(k).A3)+W2(j).*(xx(k).A2)+W3(j).*xx(k)+W4(j);endendendTend%计算实习，第一种边界条件clear;x=-l:0.2:l; %输入节点横坐标y=ones(

7、l,ll)./(ones(l,ll)+25*x.2) %计算节点纵坐标 t=sym(t) ;f=l/(l+25* t2);% 定义函数f2=diff(f,2) %函数的二阶导数式f21=vpa(subs(f2, t ,T)%计算左端点的二阶导数值f22=vpa(subs(f2, t ,1)%计算右端点的二阶导数值 xx=-1:0.1:1;T=mspline1(x,y,f21,f22,xx);T=Tezplo t(f,T 1);%画出函数f的曲线hold onplot(x,y,:,xx,T,-);%根据函数计算和插值计算的结果画出的曲线%计算实习，第二种边界条件clear;x=-l:0.2:l;

8、 %输入节点横坐标y=ones(l,ll)./(ones(l,ll)+25*x.2) %计算节点纵坐标 t=sym(t) ;f=1/(1+25* t2);% 定义函数fl=diff(f,l)%函数的一阶导数式fll=vpa(subs(fl, t ,-1)%计算左端点的一阶导数值fl2=vpa(subs(fl, t ,1)%计算右端点的一阶导数值xx=-1:0.1:1;T=mspline2(x,y,f11,f12,xx);T=Tezplo t(f,T 1);%画出函数f的曲线hold onplot (x, y, :, xx, T, -) ; %根据函数计算和插值计算的结果画出的曲线4. 计算结果第一种边界条件：x-1-0.8-0.6-0.4-0.200.20.40.60.81函数计算 值0.03850.05880.10.20.510.50.20.10.05880.0385插值计算 值0.03850.05880.10.20.510.50.20.10.05880第二种边界条件：X-1-0.8-0.6-0.4-0.200.20.40.60.81函数计算 值0.03850.05880.10.20.510.50.20.10.05880.0385插值计算 值0.03850.05880.10.20.510.50.20.10.05880