三次样条插值自然边界条件

三次样条插值函数MATLAB 编程实现三次样条插值函数为()()[)()[]1011,,,,n n n S x x x x S x S x x x x-⎧∈⎪=⎨⎪∈⎩ 利用三次埃尔米特插值函数表示三次样条插值函数，即()()()()())111111,,j j j j j j j j j j j S x y x y x m x m x x x x ααββ++++++⎡=+++∈⎣（0,1,,1j n =-）基函数满足()()()()()()21112111121121111212jj j j j j j j j j j j j j j j j j j jj j j j x x x x x x x x xx xx x x x x x xx xx x x x xx x x x x x xααββ++++++++++++⎛⎫⎛⎫--=+ ⎪⎪ ⎪⎪--⎝⎭⎝⎭⎛⎫⎛⎫--=+ ⎪⎪ ⎪⎪--⎝⎭⎝⎭⎛⎫-=-⎪ ⎪-⎝⎭⎛⎫-=-⎪ ⎪-⎝⎭由上式易得()()()()()()()()()()()()()()1331111331112211112211612612246246j j j j j j j j j j j j j j j j j j j j j j jj j j j j x x x x x x xx x x x x x x x x x x x x xx xx x x x x x x x x ααββ+++++++++++++++''=---+''=-+--+''=---+''=---则有()()()()()()()()()()()111111113333111111122221111661212242466j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j S x y x y x m x m x x x x x y x y x x x x x x x x x x x x x m x m x x x x x x x x x ααββ+++++++++++++++++++''''''''''=+++⎡⎤⎡⎤++⎢⎥⎢⎥=-+-+⎢⎥⎢⎥----⎣⎦⎣⎦⎡⎤++⎢⎥+-+-⎢⎥----⎣⎦)1,j j x x x +⎡⎤⎢⎥⎡∈⎣⎢⎥⎣⎦（0,1,,1j n =-）同理有()()()()()()()()()()()()()()()11111113333111111122221111661212242466j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j S x y x y x m x m x x x x x y x y x x x x x x x x x x x x x m x m x x x x x x x x x ααββ------------------''''''''''=+++⎡⎤⎡⎤++⎢⎥⎢⎥=-+-+⎢⎥⎢⎥----⎣⎦⎣⎦⎡⎤⎡++⎢⎥+-+-⎢⎥----⎣⎦⎣)1,j j x x x -⎤⎢⎥⎡∈⎣⎢⎥⎦（1,,j n =）根据样条函数二阶导数连续性，即()()100j j j j S x S x +''''+=-（1,,1j n =-）即()()()()()()()()()()()()()()()()111111332211111111113322111166426624j j jj j j j j j j jj j j j j j j j j j j j j j jj j j j j jj j j j j j j j x x y x x y x x x x m m x x x x x x x x x x y x x y x x x x m m x x x x x x x x ++++++++++--------------+++--------=+++----（1,,1j n =-）化简得()()()()()111111111111233j j j j j j j j j j j j j j j j j j jj j xx m x x m x x m x x x x y y y y x x x x +-+--+-++-+--+-+---=-+---（1,,1j n =-）可得线性方程组()()()()()()()()()()0121201023231213121221111110212110211032213221322122233333n n n n n n n n n n n n m m x x x x x x m x x x x x x m x x x x x x m m m x x x x y y y y x x x x x x x x y y y y x x x x y ------⨯+-+⨯⎛⎫ ⎪ ⎪ ⎪---⎛⎫⎪ ⎪--- ⎪ ⎪⎪ ⎪⎪ ⎪ ⎪--- ⎪⎝⎭ ⎪ ⎪ ⎪⎝⎭---+------+---=()()()121112112113n n n n n n n n n n n n n x x x x y y y x x x x ----------⨯⎛⎫⎪ ⎪ ⎪ ⎪⎪⎪ ⎪-- ⎪-+- ⎪--⎝⎭为了使样条插值问题有惟一解，我们在原有方程基础上增加两个边界条件。

1
样条函数的定义 定义4.1 设区间[a,b]上给定一个节点划分
a=x0<x1<……<xn-1<xn=b 如果存在正整数k使得[a,b]上的分段函数s(x)满足 如下两条： (1)在[a,b]上有直到k-1阶连续导数。 (2)在每个小区间[xi,xi+1]上是次数不大于k的多项式。 则称分段函数s(x)是以(2.6)为节点集的k次样条函数。
x xi i i 1 hi
) mi1hi
( ) xi1x
1 hi
x [xi , xi1], hi xi1 xi , i 0,1,, n 1
(x) (2x 1)( x 1)2,1(x) x(x 1)2 13
对Si (x)求二阶导数 ,并整理后得
Si( x)
6( xi
xi 1 hi3
2x)
因为分段三次Hermite插值多项式已经至少是一阶连续 可导了，为了让它成为三次样条函数只需确定节点处 的一阶导数使这些节点处的二阶导数连续即可！
S(xi 0) S(xi 0), i 1,, n 1
S(x)
y ( xxi i 0 hi
)
y ( ) m h ( xi1x i1 0 hi
( yn
yn 1 )
2 hn1
(mn1
2mn )
立即可得下式：
21
其中：
nm1 nmn1 2mn gn
n
h0
h0 hn1
, n
hn1 h0 hn1
1 n
gn
3 n
y1 y0 h0
n
yn
yn1 hn1
联合基本方程得一个广义三对角或周期三对角方程组：
2 1
1
1
2

用Matlab实现了3次样条曲线插值的算法。

边界条件取为自然边界条件，即：两个端点处的2阶导数等于0；共包含3各个函数文件，主函数所在文件（即使用的时候直接调用的函数）为spline3.m,另外两个函数文件是在splin3函数文件中被调用的自定义函数。

一个是GetParam.m,一个是GetM.m。

%GetParam.m文件的内容：%根据给定的离散点的横坐标所构成的向量，计算各个区间段的h值；function GetParam(Vx,Vy)global gh;global gf;global gu;global gr;global gd;global gff;global gM;%global gn;%n=length(Vx);%length()为向量Vx所含元素的个数；%n=legth(Vx);%gn=n;%n=gn;n=length(Vx);gh(1)=Vx(2)-Vx(1);gf(1)=(Vy(2)-Vy(1))/gh(1);for i=2:1:n-1%从区间0到区间n-1; gh(i)=Vx(i+1)-Vx(i);gf(i)=(Vy(i+1)-Vy(i))/gh(i);gu(i)=gh(i-1)/(gh(i-1)+gh(i));gr(i)=1-gu(i);gff(i)=(gf(i-1)-gf(i))/(Vx(i-1)-Vx(i+1)); gd(i)=6*gff(i);end%设置与边界条件有关的参数；gM(1)=0;%起点的2阶导数；gM(n)=0;%终点的2阶导数；end%GetM.m文件的内容：function GetM(Vx)global gh;global gf;global gu;global gr;global gd;global gff;global gM;%global gn;nn=length(Vx);%nn=gn;n=nn-2;b=zeros(n,1);A=zeros(n,n);A(1,1)=2;A(1,2)=gr(2);b(1)=gd(2)-gu(2)*gM(1);for i=2:1:n-1A(i,i)=2;A(i,i-1)=gu(i+1);A(i,i+1)=gr(i+1);b(i)=gd(i+1);endA(n,n-1)=gu(n);A(n,n)=2;b(n)=gd(nn-1)-gr(nn-1)*gM(nn); X=(inv(A))*b;for i=2:1:nn-1gM(i)=X(i-1);end%主函数文件spline3.m的内容：function result=spline3(x,Vx,Vy) global gh;global gf;global gu;global gr;global gd;global gff;global gM;%global gn;GetParam(Vx,Vy);GetM(Vx);%n=length(Vx);%n=gn;n=length(Vx);nn=length(x);y=zeros(1,nn);for j=1:1:nni=1;while(x(j)>Vx(i+1))endsn=i;t1=(Vx(sn+1)-x(j))^3/(6*gh(sn));t1=t1*gM(sn);t2=(x(j)-Vx(sn))^3/(6*gh(sn));t2=t2*gM(sn+1);t3=Vy(sn)-gM(i)*((gh(i))^2)/6;t3=t3*(Vx(sn+1)-x(j))/gh(sn);t4=Vy(sn+1)-gM(sn+1)*((gh(sn))^2)/6;t4=t4*(x(j)-Vx(sn))/gh(sn);y(j)=t1+t2+t3+t4;endresult=y;end函数调用的时候，result=spline3(x,Vx,Vy)，x为代求点的横坐标向量，(Vx,Vy)为已知的点的坐标。

(x
a)
m
m次截断多项式
a
.
7
定理5.5 任意s(x)∈Sm(x1,x2,…,xn)均可唯一地表示为
n
s(x)pm(x) cj(xxj)m , x (4-31) j1
其中pm(x)∈Pm，cj(j=1,2,…,n)为实数。
定理5.6 为使s(x)∈Sm(x1,x2,…,xn)，必须且只须存在pm(x)∈Pm
8
例1 验证分片多项式是三次样条函数。
1 2x
x 3
S ( x) 2825x9x2x3 3x1
2619x3x2x3 1x0
2619x3x2
0 x
解 利用上面的定理(光滑因子)验证．
(x 3)3,
2(x 1)3,
x3,
所以由定理5.5可知该函数为三次样条函数．
例，设
x3x2
0x1
S(x) a3xb2 xc x11x2
信息;

样? ?条?插插值值：:（样条函数—满足一定光滑性的分段多项式）。 局部性好, 满足一定光滑性, 收敛性保证, 只需要函数值
信息。
.
2
样条函数是一个重要的逼近工具，在插值、数值微分、曲 线拟合等方面有着广泛的应用。
定义5.3 对区间(-∞，+∞)的一个分割：
: x 1 x 2 x n ,
n
p n (x )p n 1 (x ) c n (x x n )m p0(x) cj(xxj)m j1
为了便于表示分段信息, 引进截断多项式：
(x a)m
(x a)m , x a,
0, x a,
(5-30)
易见
(x
a)
m
∈Cm-1(-∞,+∞)

(x)
lim
x xk
Sk 1( x)
lim
x
x
k
Sk (x)
lim
x
x
k
Sk1( x)
k 1,2,,n 1
------(4)
lim
x
x
k
Sk( x)
lim
x
x
k
Sk1( x)
共4n 2个条件
5
Sk (x)是[xk , xk 1 ]上的三次样条插值多项式,应有4个待定的系数 即要确定S(x)必须确定4n个待定的系数 少两个条件 并且我们不能只对插值函数在中间节点的状态进行限制 也要对插值多项式在两端点的状态加以要求 也就是所谓的边界条件:
例. 使用不同的插值方法于函数
y
1
1 x2
x [5,5]
最后，介绍一个有用的结论
定理 . 设f (x) C 2[a,b], S(x)是以xk (k 0,1,, n)
为节点, 满足任意边界条件的三次样条插值函数，
设hi
xi 1
xi
,
h
max
0in1
hi
,
min
0in1
hi
,
则当 h
c 时
S(x)和S(x)在[a,b]上一致收敛到f (x)和f (x)
------(6)
13
由(11)式,可知
S0( x0
)
6( x0
x1 h03
2 x0
) ( y1
y0 )
6 x0
2 x0 h02
4 x1
m0
6 x0
4 x0 h02
2 x1
m1
6 h02
(

三次样条插值自然边界条件(一)三次样条插值自然边界条件引言•什么是样条插值？•为什么需要自然边界条件？样条插值基本概念•样条定义和性质•插值问题三次样条插值算法1.数据准备–输入数据的获取和整理2.线性方程组的建立–利用数据点和导数的关系建立线性方程组3.矩阵求解–解线性方程组获取样条插值函数的系数4.插值函数构建–使用线性方程组的解构建三次样条插值函数自然边界条件的引入•什么是自然边界条件？•为什么要使用自然边界条件？自然边界条件的应用1.边界条件的设置–左、右边界条件的定义2.线性方程组建立–利用自然边界条件建立线性方程组3.矩阵求解–解线性方程组获取带有自然边界条件的样条插值函数的系数4.插值函数构建–使用线性方程组的解构建带自然边界条件的三次样条插值函数结论•自然边界条件在三次样条插值中的应用优势•对比其他边界条件的选择参考文献参考文献1.DeBoor, C. (1978). A Practical Guide to Splines. NewYork: Springer-Verlag.2.王援胜. 数值计算方法（第2版）. 中国科学技术大学出版社,2014.3.徐健, 曹洪欣, & 张平原. (2003). 数值分析与算法 C++语言实现. 清华大学出版社.4.Zhang, Q., & Jiao, L. (2001). Piecewise polynomialinterpolation with adaptive knot values. Journal ofcomputational and applied mathematics, , .。

三次样条插值多项式实验的目的及意义：为了取得理想结果：在不增加更多的插值条件下，能够求得一个插值多项式，既有良好的逼近效果，又有好的光滑性，引进三次样条插值 多项式。

如果已知函数y=f(x)在节点a=x0<x1<…<xn=b 处的函数值和导数值：()i i x f y =，i=0,1,2,…,n如果S(x)满足条件：1. S(x)是一个分段的三次多项式且()i i y x S =；2. S(x)在[a,b]具有二阶连续导数。

则称S(x)是三次样条插值函数。

S(x)的具体形式为：()()()()⎪⎪⎩⎪⎪⎨⎧∈∈∈=-]12,121,01,[,...............][,][,n n n x x x x s x x x x s x x x x s x s其中()x S i 在[]i i x x ,1-上是三次多项式()iiiiid x c x b x a x S +++=23由插值条件()ii y x S =，i=0,1,2,…,n ，得n+1个条件。

边界条件一：()()nn y x S y x S '',''00== 边界条件二：()()nn y x S y x S '''',''''00==数学公式：()()2211133[2]()[2()]()i i i i i i i i i i ih x x x x h x x x x H x y y h h ---+-----=++2211122()()()()i i i i i i i ix x x x x x x x m m h h -------+ 算法描述：Step1:输入未知数X 及(xi,yi),i=0,1,…,n ； Step2:计算步长H[i]; Step3:计算[][]()⎪⎪⎪⎩⎪⎪⎪⎨⎧+=-=+=-+++i i i i i i i ii i i i i x x f x x f u g u hh h ,,311111λλλStep4:根据边界条件，求解相应的方程得到m0,m1,…, mn Step5:判断X 属于[]i i x x ,1-，i=1,2,…,n 中的哪一个Step6:计算()x s y i i ≈Step7:输出y. 程序原代码如下： #include "stdio.h" #define N 5 void main() { int i,k; float X,s,y0,yn;float a[N][N+1],h[N],u[N],v[N],g[N],m[N],p[N],q[N],w[N];printf("please input X:"); //X 为未知数的大小scanf("%f",&X);printf("please input x:"); //输入x的大小for(i=0;i<N;i++)scanf("%f",&a[i][0]);printf("please input y:"); //输入y的大小for(i=0;i<N;i++)scanf("%f",&a[i][1]);for(i=1;i<N;i++)h[i]=a[i][0]-a[i-1][0]; //计算步长for(i=1;i<N;i++){v[i]=h[i+1]/(h[i]+h[i+1]);u[i]=1-v[i];g[i]=3*u[i]*(a[i+1][1]-a[i][1])/h[i+1]+3*v[i]*(a[i][1]-a[i-1][1])/h[i]; }printf("\t(1)已知边界条件1\n");printf("\t(2)已知边界条件2\n");printf("请选择边界条件序号:");scanf("%d",&k);if(k==1){printf("请输入y0和yn的一阶导:"); //输入边界条件一scanf("%f%f",&m[0],&m[N-1]);p[0]=0; //用追赶法求解m[N]q[0]=0;g[1]=g[1]-v[1]*m[0];g[N-2]=g[N-2]-u[N-2]*m[N-1];for(i=1;i<N;i++){w[i]=2-u[i]*p[i-1];p[i]=v[i]/w[i];q[i]=(g[i]-u[i]*q[i-1])/w[i];}m[N-2]=q[N-2];for(i=N-3;i>0;i--)m[i]=q[i]-p[i]*m[i+1];printf("输出m[i]的值:\n");for(i=0;i<N;i++)printf("%f\n",m[i]);for(i=1;i<N;i++) //计算最终结果if(X>a[i-1][0]&&X<a[i][0])s=(h[i]+2*(X-a[i-1][0]))*(X-a[i][0])*(X-a[i][0])*a[i-1][1]/(h[i]*h[i]*h[i]) +(h[i]-2*(X-a[i][0]))*(X-a[i-1][0])*(X-a[i-1][0])*a[i][1]/(h[i]*h[i]*h[i])+ (X-a[i-1][0])*(X-a[i][0])*(X-a[i][0])*m[i-1]/(h[i]*h[i])+(X-a[i-1][0])*(X-a[i-1][0])*(X-a[i][0])*m[i]/(h[i]*h[i]);printf("s(%f)=%f\n",X,s);}if(k==2){printf("请输入y0和yn的二阶导:"); //输入边界条件二scanf("%f%f",&y0,&yn);g[0]=3*(a[1][1]-a[0][1])/h[1]-h[1]*y0/2;g[N-1]=3*(a[N-1][1]-a[N-2][1])/h[N-1]+h[N-1]*yn/2;q[0]=g[0];u[0]=1;v[N-1]=1;w[0]=2;for(i=1;i<N;i++){w[i]=2-v[i]*u[i-1]/w[i-1];q[i]=g[i]-v[i]*q[i-1]/w[i-1];}m[N-1]=q[N-1]/w[N-1];for(i=N-2;i>=0;i--)m[i]=(q[i]-u[i]*m[i+1])/w[i];printf("输出m[i]的值:\n");for(i=0;i<N;i++)printf("%f\n",m[i]);for(i=1;i<N;i++)if(X>=a[i-1][0]&&X<=a[i][0])s=(h[i]+2*(X-a[i-1][0]))*(X-a[i][0])*(X-a[i][0])*a[i-1][1]/(h[i]*h[i]*h[i]) +(h[i]-2*(X-a[i][0]))*(X-a[i-1][0])*(X-a[i-1][0])*a[i][1]/(h[i]*h[i]*h[i])+ (X-a[i-1][0])*(X-a[i][0])*(X-a[i][0])*m[i-1]/(h[i]*h[i])+(X-a[i-1][0])*(X-a[i-1][0])*(X-a[i][0])*m[i]/(h[i]*h[i]);printf("s(%f)=%f\n",X,s);}}数值计算：已知y=f(x)的如下数值求三次样条插值函数S（x），满足条件1.s’(0)=0,s’(4)=482.s’’(0)=0,s’’(4)=24Please input X:2.5Please input x:0 1 2 3 4 Please input y:-8 -7 0 19 56(1)已知边界条件1(2)已知边界条件2请选择边界条件的序号:1请输入y0和yn的一阶导:0 48 0.0000003.00000012.00000027.00000048.000000s(2.500000)=7.625000press any key tocontinue请选择边界条件的序号:2请输入y0和yn的二阶导:0 24 -0.0000003.00000011.99999927.00000248.0000007.625000press any key tocontinue s(2.500000)=7.625000 对计算结果进行评价分析：()()443845h M x S x f ≤-三次样条插值函数与三次Hermite 插值函数相比，不仅光滑度有提高，而且要求求解时还不需要增加内节点处的导数值，因此比较实用。

6( xi xi 1 2 x ) ( yi 1 yi ) 3 hi 1
而
2 4 6 S ( xi 0) mi 1 mi 2 ( yi yi 1 ) hi hi hi 4 2 6 S ( xi 0) mi m i 1 2 ( yi 1 yi ) hi 1 hi 1 hi 1
n
当n 时，Ln ( x )只在 | x | 3.63 内收敛，而在该区间外 是发散的。
从图中可以看出，在 0 附近插值效果是好的，即余项较 小，另一种现象是插值多项式随节点增多而振动。这种插值 多项式当节点增加时反而不能更好地接近被插值函数的现象， 称为龙格现象。
上述现象告诉我们用高次插值多项式是不 妥当的，从数值计算上可解释为高次插值多项 式的计算会带来舍入误差的增大，从而引起计 算失真。因此，实践上作插值时一般只用一次、 二次最多用三次插值多项式。
式中x [ xi 1 , xi ] (i 1,2,, n)
第(2)步
为了确定mi，需要用到S ( x )的二阶导数在节点连续 的条件， S ( x )在[ xi 1 , xi ]和[ xi , xi 1 ]上的二阶导数分别为
Si( x ) 6 x 2 xi 1 4 xi 6 x 4 xi 1 2 xi mi 1 mi 2 2 hi hi ( x [ xi 1 , xi ])
若记hi xi xi 1，则上式可写为
( x x i ) 2 hi 2( x x i 1 ) ( x x i 1 ) 2 hi 2( x i x ) Si ( x) y i 1 yi 3 3 hi hi ( x x i ) 2 ( x x i 1 ) ( x x i 1 ) 2 ( x x i ) m i 1 mi 2 2 hi hi

hk hk 1
mk 1 2mk
hk 1 hk hk 1
k 1
3( hk yk1 yk hk1 yk yk1 )
hk hk 1
hk
hk hk 1
hk 1
k mk1 2mk k mk 1 gk
k
hk
hk hk 1
k
hk 1 hk hk 1
gk
3(k
yk yk 1 hk 1
k 0,1,2
小结
1 x3 3 x2 7 x 1
8 8 4
1 x2
S(x)
1 8
x3
3 8
x2
7 4
x
1
3 x3 45 x2 103 x 33
88
4
2x4 4x5
f (3) S(3) 17 4
最后，介绍一个有用的结果
定理 . 设f (x)C2[a,b],S(x)是以xk (k 0,1,,n)
m2 m3
g0 g1 g2
g3
解方程组得：m0
17 8
, m1
Байду номын сангаас
7 4
, m2
5 4
, m3
19 8
将上述结果代入(10)式
S0 ( x)
1 8
x3
3 8
x2
7 4
x
1
S1 ( x)
1 8
x3
3 8
x2
7 4
x
1
1 x 2 2x4
S2(x)
3 8
x3
45 8
x2
103 4
x
33
4x5
注：三次样条与分段 Hermite 插值的根本区别在于S(x)自 身光滑，不需要知道 f 的导数值（除了在2个端点可能需 要）；而Hermite插值依赖于f 在所有插值点的导数值。

确定S ( x )必须确定4n个待定的系数的条件正好也是4n个
即
lim S ( x ) = lim S ( x ) lim S ′ ( x) = lim S ′ ( x) = m lim S ′′( x) = lim S ′′ ( x) S ′( x ) = f ′ S ′( x ) = f ′ 或
将(13)式化为矩阵形式
2 λ2
µ1 λ3
2
µ2 λ4
2
µ3
2 O O O
λn − 2
λn − 1
O 2
m1 g 1 − λ1 f 0′ m2 g2 g3 m3 M = M M M µ n − 2 mn − 2 gn−2 2 mn − 1 g n − 1 − µ n − 1 f n′
共4 n − 2个条件
′ ′ lim S k′( x ) = lim S k′− 1 ( x ) + −
S k ( x )是[ xk , xk + 1 ]上的三次样条插值多项式, 应有4个待定的系数 即要确定S ( x )必须确定4n个待定的系数
少两个条件
并且我们不能只对插值函数在中间节点的状态进行限制 也要对插值多项式在两端点的状态加以要求 也就是所谓的边界条件: 第一类(一阶)边界条件: 第二类(二阶)边界条件 第三类(周期)边界条件
f(x) H(x) S(x)
二、三次样条插值多项式
a ≤ x0 , x1 ,L , xn ≤ b为区间[ a , b ]的一个分割 如果函数f ( x )在节点x0 , x1 ,L , xn处的函数值为
f ( x j ) = y j , j = 0 ,1,L , n 如果S ( x )是f ( x )的三次样条插值函数, 则其必满足

则有: S(-1)= – a1+b1–c1+d1=f(-1)=1, S(1)=a2+b2+c2+d2= f (1) =1,
S(0)=d1= f(0)=0,
S(0-0)= d1=S(0+0)=d2, S'-(0)= c1= S'+(0)= c2, S''-(0)=b1= S''+(0)=b2 由自然边界条件: S''(0)= – 6a1+2b1=0, S'(1)= 6a2+2b2=0 解方程组,得 a1=-a2=1/2, b1=b2=3/2, c1=c2=d1=d2=0
令 S( x j ) M j，（j 0,1,2,, n)
参数
插值函数， 因为S j ( x)是三次样条
所以, S [ x j , x j 1 ]上是一次函数,由两点拉格朗日插值可表示为 j ( x)在 x j 1 x x xj S ( x) Mj M j 1 , x [ x j , x j 1 ], hj x j 1 x j hj hj (2.46) 对上式积分,得 ( x j 1 x)2 ( x x j )2 S ( x) Mj M j 1 c1 , (2.47) 2h j 2h j
( xi , f ( xi )),(i 0,1,, n)且 a x0 x1 xn b ;
f ( x ) 于 [a , b] 存在 ，则 ( a （或（ ) b ）或（ c ）） (2)给定边界条件
) b）或（ c））。 唯一3次样条插值函数 S ( x ),且满足 (a（或（
第二章 插值法
例 2.13 已知 f (–1) = 1, f (0) = 0, f (1) = 1.求[–1，1] 上的三次自然样条(满足自然边界条件). 解设

三次样条插值1. 算法原理由于在许多实际问题中，要求函数的二阶导数连续，人们便提出了三次样条插值函数，三次样条插值函数是由分段三次函数拼接而成的，在连接点处二阶导数连续。

设S(x)在节点i x 处的二阶导数),1,0()(''n i M x S i i ，⋯⋯==，其中i M 为待定参数。

由S （x ）是分段三次多项式可知，)(''x S 是分段线性函数，)(''x S 在子区间[]i i x x ,1-上可以表示为ii i ii i i i i i i i i i i i x x x M h x x M h x x M x x x x M x x x x x S ≤≤-+-=--+--=-------1111111,)(''其中1--=i i i x x h ，对S （x ）两端积分两次得()ii i i i i i i i i i i x x x x x c x x b M h x x M h x x x S ≤≤-+-+-+-=-----111113),(6)(6)()(其中i b 和i c 为积分常数。

由插值条件()i i i i y x S y x S ==--)(,11得i i i i i i i i i i y h c M h y h b M h =+=+--6,62112由此解得i i i i i i i i i i h M h y c h M h y b /6,/62121⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎭⎫ ⎝⎛-=--代入得：()()()ii ii i i i i i i i i i ii i i i x x x h x x M h y h x x M h y Mh x x Mh x x x S ≤≤-⎪⎪⎭⎫ ⎝⎛-+-⎪⎪⎭⎫ ⎝⎛-+-+-=------1121213113,6666 求导得：()()()i i i i i i i i i ii i ii x x x h M M h y y Mh x x Mh x x x S ≤≤---+-+--=-----1112112,622' 令i x x =得()x S 在i x 处的左导数 ()ii i i i i i i i i i i i i i h y y M h M h h M M h y y M h x S 11113662'------++=---+=又令1-=i x x 得()x S 在1-i x 处右导数 ()ii i i i i i i i i i i i i i i h y y M h M h h M M h y y M h x S 1111116362'------+-+--=---+-=， 从而有1111163)('++++++-+--=i ii i i i i h y y M h M h x S ，由()x S 在节点i x 处一阶导数的连续性知1,2,1),(')('-⋯⋯==+-n i x S x S i i ，，即1,2,16361111111-⋯⋯=---=+++-+++++-n i h y y h y y M h M h h M h ii i i i i i i i i i i i ，， 两端同乘16++i i h h 得)(62111111111ii i i i i i i i i i i i i i i i h y y h y y h h M h h h M M h h h -++++++-+---+=++++，记ii i i i i i i i h h h h h h μλμ-=+=+=+++1,111，1,,2,1],,,[66111111-⋯⋯==⎪⎪⎭⎫⎝⎛---+=+--+++n i x x x f h y y h y y h h d i i i i i i i i i i i i ，则关于i M 的方程组写成1,2,1,211-⋯⋯==+++-n i d M M M i i i i i ，λμ。

若边界条件为 S ′′( x 0 ) = S ′′( x n ) = 0 ，即满足自然边界条件,则得两端的方程为
2m0 + m1 = 3 f [ x0 , x1 ] = g 0 ; mn −1 + 2mn = 3 f [ xn −1 , xn ] = g n .
于是，用矩阵形式表为
2 1 0 λ 2 µ 1 1 0 λ2 2 L L L 0 0 0 0 0 0
用
y j +1 − y j 1 1 = f [ x j , x j +1 ] ，上面方程可简化为 + 除全式，并注意 y j = f j , h j −1 h j hj
λ j m j +1 + 2m j + µ j m j +1 = g j
其中
( j = 1, 2, L , n − 1);
λj =
hj h j −1 + h j
( j = 0, 1, L , n) ，共有 4n − 2
，可根据实际 通常可在区间 [a, b] 端点 a = x0 , b = xn 上各加一个条件（称为边界条件） 问题的要求给定。常见的有以下三种： 1° 已知两端的一阶导数值，即
S ′( x0 ) = f 0′, S ′( xn ) = f n′ .
g1 − λ1 f 0′ g2 g3 = . M g n−2 g n −1 − µ n −1 f n′
如果边界条件为 S ′′( x0 ) = f 0 ,
''
S ′′( xn ) = f n′′ ，则得两个方程
h 2m0 + m1 = 3 f [ x0 , x1 ] − 0 f 0′′ = g 0 ; 2 m + 2m = 3 f [ x , x ] + hn −1 f ′′ = g . n −1 n n −1 n n n 2

一、问题提出
为给定的节点， 设 x0 , x1 xn 为给定的节点，yi = f ( xi ) ，i = 0,1, n 为相应的函数值， 为相应的函数值，求一个次数不超过 n 的多项式 Pn (x)， 使其满足
Pn ( xi ) = yi，
i = 0,1, n .
这类问题称为插值问题。 称为被插值函数 P 被插值函数， 这类问题称为插值问题。 f ( x) 称为被插值函数， n ( x) 称 插值问题 插值函数， 称为插值节点 为插值函数， x0 , x1 xn 称为插值节点
六、 分段插值
所谓分段插值，就是将被插值函数逐段多项式化。 所谓分段插值，就是将被插值函数逐段多项式化。在每 个 [ xi , xi +1 ] 子段上构造插值多项式，然后把它们装配在一， 子段上构造插值多项式，然后把它们装配在一， 作为整个区间 [ a, b ] 上的插值函数，即称为分段多项式。如果 上的插值函数，即称为分段多项式。 次式， 函数 Sk ( x ) 在每个子段上都是 k 次式，则称为 一般（低次： 一般（低次：k=1,2,3） ） 次式。 k 次式。
f [ x0 , x1 ] = 5, f [ x0 , x1, x2 , x3 ] = 1,
N n ( x) = 0 5( x 1) + 2( x 1)( x 2)
+ ( x 1)( x 2)( x 3)
= x3 4 x + 3
五、 Hermite插值多项式 插值多项式
给定的是节点上的函数值和导数值 问题： 问题：已知
∑
i=0
y i li ( x )
( x x0 ) ( x xi 1 )( x xi +1 ) ( x xn ) , i = 0,1, n ( xi x0 ) ( xi xi 1 )( xi xi +1 ) ( xi xn )

/****************函数说明*********************///pxpy为已知的数据点，xs为要插值的x坐标，最终会得到xs坐标下的y值using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace spline{class Program{staticvoid Main(string[] args){point[] points = new point[13];double[] px = { 64, 304, 544, 1035, 1502, 2061, 2540, 2939, 3498, 3897, 4456, 4696, 4936 }; double[]py={4.663,4.476,4.969,4.853,4.766,5.2415,4.795,4.7055,5.4565,5.5515,5.411,5.8385,6.7085} ;for (int i = 0; i< 13; i++){points[i] = new point();points[i].x =px[i];points[i].y = py[i];}point.DeSortX(points);double[] xs = {64, 304, 544, 1023, 1502, 1981, 2460, 2939, 3418, 3897, 4376, 4616, 4856}; splineInsertPoint(points, xs, 1);Console.ReadLine();}staticdouble[] splineInsertPoint(point[] points, double[] xs, int chf){int plength = points.Length;double[] h = newdouble[plength];double[] f = newdouble[plength];double[] l = newdouble[plength];double[] v = newdouble[plength];double[] g = newdouble[plength];//三转角法的待定一阶系数法。

S (x0 ) S (xn ), S (x0 ) S (xn )
这样，由上给定的任一种边界条件加上插值条件
和连接条件，就能得出4n个方程，可以惟一确 定4n个系数。从而得到三次样条插值函数S(x) 在各个子区间xi , xi+1上的表达式S(xi） (i=1,2,…,)。但是，这种做法当n较大时，计 算工作很大，不便于实际应用。因此我们希望找
三次样条插值函数S(x)是一个分段三次多项式,要求出 S(x),在每个小区间xi,xi+1上要确定4个待定参数,若用 Si(x)表示它在第i个子区间xi,xi+1上的表达式，则
Si (x) ai0 ai1x ai2 x2 ai3 x3 i 0,1, , n 1
其中四个待定系数为 ai0 , ai1, ai2 , ai3 ,子区间共有n个 所以要确定S(x)需要4n个待定系数。
第二种类型：给定两端点f(x)的二阶导数值：
S (x0 ) f (x0 ), S (xn ) f (xn )
作为特例,S(x0 ) S(xn ) 0 称为自然边界条件。满 足自然边界条件的三次样条插值函数称为自然样 条插值函数。
第三种类型：当f(x)是以为 xn x0 周期的函数时， 则要求S(x)也是周期函数,这时边界条件应满足 当 f (x0 ) f (xn ) 时，
i 1,2, , n 1
上述二式共给出了4n-2个条件,而待定系数有4n个,因此 还需要2个条件才能确定S(x),通常在区间端点上
a x0 ,b xn 各加一个条件,称为边界条件, 常用边 界条件有三种类型。
第一种类型：给定两端点f(x)的一阶导数值：
S (x0 ) f (x0 ), S (xn ) f (xn )
另一方面,要求分段三次多项式S(x)及其导数 S(x) 和 S(x) 在整个插值区间a,b上连续,则要求它们在 各个子区间的连接点 x0 , x1 , , xn1 上连续， 即满足条件

画图：
x1=0:.01:1;y1=polyval(S1(1,:),x1-X(1)); x2=1:.01:2;y2=polyval(S1(2,:),x2-X(2)); x3=2:.01:3;y3=polyval(S1(3,:),x3-X(3)); x4=3:.01:4;y4=polyval(S1(4,:),x4-X(4)); x5=4:.01:5;y5=polyval(S1(5,:),x5-X(5)); x6=5:.01:6;y6=polyval(S1(6,:),x6-X(6)); >> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.') >> hold on >> x1=0:.01:1;y1=polyval(S2(1,:),x1-X(1)); x2=1:.01:2;y2=polyval(S2(2,:),x2-X(2)); x3=2:.01:3;y3=polyval(S2(3,:),x3-X(3)); x4=3:.01:4;y4=polyval(S2(4,:),x4-X(4)); x5=4:.01:5;y5=polyval(S2(5,:),x5-X(5)); x6=5:.01:6;y6=polyval(S2(6,:),x6-X(6)); >> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.') >> hold on >> x1=0:.01:1;y1=polyval(S3(1,:),x1-X(1)); x2=1:.01:2;y2=polyval(S3(2,:),x2-X(2)); x3=2:.01:3;y3=polyval(S3(3,:),x3-X(3)); x4=3:.01:4;y4=polyval(S3(4,:),x4-X(4)); x5=4:.01:5;y5=polyval(S3(5,:),x5-X(5)); x6=5:.01:6;y6=polyval(S3(6,:),x6-X(6)); >> plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,X,Y,'.')