三次样条曲线

一、样条函数的定义

样条函数属于分段光滑插值，他的基本思想是，在由两相邻节点所构成的每一个小区间内用低次多项式来逼近，并且在各结点的连接处又保证是光滑的(即导数连续)。

设在区间[a,b]上给定一组结点X:

,和一组对应的函数值。若函数S(x)满足下列条件：

(1)在每一个子区间(k=1,2,...n)上，S(x)是一个不超过三次的多项式。

(2)在每一个结点上满足S(xi)=yi,i=0,1,...,n。

(3)S(x)在区间[a,b]上为二次连续可微函数。

则称S(x)为结点X上插值与Y的三次样条插值。

二、三次样条函数的构造

在工程上，构造三次样条插值函数通常有两种方法：一是以给定插值结点处得二阶导数值作为未知数来求解，而工程上称二阶导数为弯矩，因此，这种方法成为三弯矩插值。二是以给定插值结点处得一阶导数作为未知数来求解，而一阶导数右称为斜率，因此，这种方法称为三斜率插值。

三斜率插值法

根据定义，三次样条函数在插值结点处一阶导数应存在。因此设各结点处的一阶导数为：

。利用两点埃尔米特插值公式，就可以得到样条插值函数S(x)在子区间上的表达式为：

 (1)

其中：。由此可知只要确定各结点的一阶导数值mk(k=0,1,2....n)，则各子区间上的三次样条插值函数S(x)也就确定了。

由于S(x)在区间[a,b]上的二阶导数是连续的，即在各结点的左右两子区间上的S(x)虽然不同，但在连接点的二阶导数存在，即在连接点处的二阶左导数与二阶右导数相等：。分别求子区间[xk-1,xk]右端点xk上的二阶左导数，以及子区间[xk,xk+1]左端点xk上的右导数，即可得到：

，(2)与 ，(3) 整理后可得：，k=1,2,...n-1,(4)。

其中，并且ak+bk=1。由于有n+1个插值点，因此需要确定n+1个m(m0,m1,...,mn)。而现在方程组只有n-1个方程(公式(4))，因此还需要另外补充两个条件财能唯一确定n+1个m。在实际应用中，这两个条件可以根据实际的物理意义所给出的边界条件来得到。实际问题中常用的三种边界条件如下：

(1)给定区间两个端点处的一阶导数值，即

则根据方程组(4)可得到新的方程组：

方程组的系数矩阵为，可知此矩阵为三对角矩阵。因此方程组可用追赶法求解。

解出mk后，即得到各结点的一阶导数值，将mk带入各结点的二阶导数值得表达式公式(2)或(3)可求得各结点的二阶导数值，将mk带入各区间上S(x)的表达式公式(1)即可得到各子区间上的三次样条插值函数S(x)。

(2)给定区间两个端点处得二阶到数值，即

，由此可得两个补充方程，其中。与公式(4)联立可得如下方程组:

(3)插值函数为周期函数

yn=y0，且，由此可得两个补充方程为，其中

，并且an+bn=1。最后可得方程组：

在此给出第一种边界条件下三次样条插值的算法实现：

[cpp]  view plain copy

1. <pre name="code" class="cpp"><pre name="code" class="cpp"><pre name="code" class="cpp">//三次样条插值  
2. //x,y:[in] 已知点对应的函数值  
3. //N:[in] 已知点的个数  
4. //dy0:[in] x[0]的一阶导数  
5. //dyn:[in] x[n-1]的一阶导数  
6. //t:[in] 待求点的值  
7. //z:[out] t中的点所对应的函数值  
8. //nLen:[in] t与z的长度  
9. void Spline3(double *x,double *y,int N,double dy0,double dyn,double *t,double *z,int nLen)  
10. {  
11.     double h0,h1;  
12.     double ak,bk,ck;  
13.     double *pBuffer;  
14.     double *dy;  
15.     CBaseFunc::Allocate(pBuffer,(N+1)*N);  
16.     CBaseFunc::Allocate(dy,N);  
17.   
18.     memset(pBuffer,0,sizeof(double)*(N+1)*N);  
19.     memset(dy,0,sizeof(double)*N);  
20.     h0=x[1]-x[0];  
21.     h1=x[2]-x[1];  
22.     ak=h1/(h0+h1);  
23.     bk=h0/(h0+h1);  
24.     ck=3*(h0*(y[2]-y[1])/h1+h1*(y[1]-y[0])/h0)/(h0+h1);  
25.   
26.     //求方程组的系数矩阵  
27.     pBuffer[0*(N+1)+0]=1;  
28.     pBuffer[0*(N+1)+N]=dy0;  
29.   
30.     pBuffer[1*(N+1)+1]=2;  
31.     pBuffer[1*(N+1)+2]=bk;  
32.     pBuffer[1*(N+1)+N]=ck-ak*dy0;  
33.   
34.     for(int i=2;i<N-2;i++)  
35.     {  
36.         h0=h1;  
37.         h1=x[i+1]-x[i];  
38.         ak=h1/(h0+h1);  
39.         bk=h0/(h0+h1);  
40.         ck=3*(h0*(y[i+1]-y[i])/h1+h1*(y[i]-y[i-1])/h0)/(h0+h1);  
41.   
42.         pBuffer[i*(N+1)+i-1]=ak;  
43.         pBuffer[i*(N+1)+i]=2;  
44.         pBuffer[i*(N+1)+i+1]=bk;  
45.         pBuffer[i*(N+1)+N]=ck;  
46.     }  
47.     h0=h1;  
48.     h1=x[N]-x[N-1];  
49.     ak=h1/(h0+h1);  
50.     bk=h0/(h0+h1);  
51.     ck=3*(h0*(y[N]-y[N-1])/h1+h1*(y[N-1]-y[N-2])/h0)/(h0+h1);  
52.     pBuffer[(N-2)*(N+1)+N-3]=ak;  
53.     pBuffer[(N-2)*(N+1)+N-2]=bk;  
54.     pBuffer[(N-2)*(N+1)+N]=ck-bk*dyn;  
55.   
56.     pBuffer[(N-1)*(N+1)+N-1]=1;  
57.     pBuffer[(N-1)*(N+1)+N]=dyn;  
58.     //追赶法解三对角矩阵,解为各插值点的一阶导数值  
59.     TridiagonalMatrix(pBuffer,N,dy);  
60.   
61.     int nPos=0;  
62.     memset(z,0,sizeof(double)*nLen);  
63.     for(int k=0;k<nLen;k++)  
64.     {  
65.         //寻找各待插值点的位置  
66.         if(t[k]>=x[N-1])  
67.         {  
68.             nPos=N-2;  
69.         }  
70.         else  
71.         {  
72.             nPos=0;  
73.             while(nPos+1<N)  
74.             {  
75.                 if(t[k]>=x[nPos] && t[k]<=x[nPos+1])   
76.                 {  
77.                     break;  
78.                 }  
79.                 nPos++;  
80.             }  
81.         }  
82.         //根据公式计算待插值点的值  
83.         if(nPos+1<N)  
84.         {  
85.             h0=x[nPos+1]-x[nPos];  
86.             z[k]=dy[nPos]*pow(t[k]-x[nPos+1],2.0)*(t[k]-x[nPos])/pow(h0,2.0)+  
87.                 dy[nPos+1]*pow(t[k]-x[nPos],2.0)*(t[k]-x[nPos+1])/pow(h0,2.0)+  
88.                 y[nPos]*pow(t[k]-x[nPos+1],2.0)*(2*(t[k]-x[nPos])+h0)/pow(h0,3.0)+  
89.                 y[nPos+1]*pow(t[k]-x[nPos],2.0)*(-2*(t[k]-x[nPos+1])+h0)/pow(h0,3.0);  
90.         }  
91.     }  
92.     CBaseFunc::Free(pBuffer);  
93.     CBaseFunc::Free(dy);  
94. }  
95.   
96. //追赶法解三对角方程  
97. //a[k][k+1]=a[k][k+1]/a[k][k]  
98. //a[k+1][k+1]-=a[k+1][k]*a[k][k+1]  
99. bool TridiagonalMatrix(double *a,int N,double *result)  
100. {  
101.     for(int k=0;k<N-1;k++)  
102.     {  
103.         a[k*(N+1)+k+1]/=a[k*(N+1)+k];  
104.         a[k*(N+1)+N]/=a[k*(N+1)+k];  
105.   
106.         a[(k+1)*(N+1)+k+1]-=a[(k+1)*(N+1)+k]*a[k*(N+1)+k+1];  
107.         a[(k+1)*(N+1)+N]-=a[(k+1)*(N+1)+k]*a[k*(N+1)+N];  
108.     }  
109.   
110.     //回带求解  
111.     result[N-1]=a[(N-1)*(N+1)+N]/a[(N-1)*(N+1)+N-1];  
112.     for(int k=N-1-1;k>=0;k--)  
113.     {  
114.         result[k]=(a[k*(N+1)+N]-a[k*(N+1)+k+1]*result[k+1]);  
115.     }  
116.     return true;  
117. }  
118.   
119.   
120. //测试函数  
121. void Spline3()  
122. {  
123.     const int N=12;  
124.     double x[N]={0.52,8,17.95,28.65,50.65,104.6,156.6,260.7,364.4,468,507,520};  
125.     double y[N]={5.28794,13.84,20.2,24.9,31.1,36.5,36.6,31,20.9,7.8,1.5,0.2};  
126.     double t[8]={4,14,30,60,130,230,450,515};  
127.     double z[8];  
128.   
129.     double dy0=1.86548;  
130.     double dyn=-0.046115;  
131.   
132.     Spline3(x,y,N,dy0,dyn,t,z,8);  
133. }  
134. </pre><pre name="code" class="cpp"></pre><pre name="code" class="cpp"></pre><pre name="code" class="cpp"></pre><pre></pre>  
135. <pre></pre>  
136. <pre></pre>  
137. <pre></pre>  
138. <pre></pre>  
139. <pre></pre>  
140. <pre></pre>  
141. <pre></pre>  
142. <pre></pre>  
143.   
144. </pre></pre>