三次样条插值——三转角方程的算法设计

三转角样条插值摘要:对于插值多项式的次数问题，人们一般认为：插值多项式的次数越高，精度越高。

比如线性插值的误差就比抛物线插值的误差要大。

但是在20世纪初龙格（Runge ）现象说明插值多项式并非次数越高精度越高。

人们利用样条插值来拟合曲线，可以利用足够的的数据点，还有公式简单，运算量节省等好处，其中最常用的样条函数是三次样条。

在课本《现代数字数学与计算》中采用的是三弯矩方法推导得出三次样条插值多项式的计算公式，本文将采用三转角法推出三次样条插值多项式的计算公式，并用其拟合几个函数。

在发现等距节点在拟合不连续函数的结果不让人满意时，使用Chebyshev 插值节点构造不均匀节点来减少其拟合的误差，最后用三转角样条插值函数拟合画出手写字母。

三转角法的理论：对于b t t t a n =<<<=...10，考虑],[2b a C s ∈的三次样条插值函数，使得n i t f f t s i i i ,...,1,0),()(===，且满足第一类边界条件00'0')()(β==t f t s ，n n n t f t s β==)()(''。

验证在区间n i t t i i ,...,2,1],,[1=-上满足)()()()()(11t v t u t q f t p f t s i i i i i i i i i ββ+++=--其中1--=i i i t t h ，)](2[)()(132--+-=i i ii i t t h ht t t p ， )](2[)()(321i i ii i t t h h t t t q ---=- ，212)()()(ii i i h t t t t t u ---=， 221)()()(ii i i h t t t t t v --=- ，而 i β可由递推公式)(3)11(31)11(21212112121111ii i i i ii i i i i ii ih f h f h h f h h h h -++++++--+-=+++βββ所确定。

三次样条插值的求解摘要：分段低次插值虽然解决了高次插值的振荡现象和数值不稳定现象，使得插值多项式具有一致收敛性，保证了插值函数整体的连续性，但在函数插值节点处不能很好地保证光滑性要求，这在某些要求光滑性的工程应用中是不能接受的。

如飞机的机翼一般要求使用流线形设计，以减少空气阻力，还有船体放样等的型值线，往往要求有二阶光滑度（即有二阶连续导数）。

因此，在分段插值的基础上，引进了一种新的插值方法，在保证原方法的收敛性和稳定性的同时，又使得函数具有较高的光滑性的样条插值。

关键字:三转角方程 三弯矩阵方程0. 引言1，三次样条函数定义1：若函数2()[,]S x a b C ∈,且在每个小区间上1,j j x x +⎡⎤⎦⎣上是三次多项式，其中01n a x x x b ⋯=<<<= 是给定节点，则称()s x 是节点01,,,n x x x ⋯上的三次样条函数。

若节点j x 上 给定函数值()j j y f x =(0,1,)j n ⋯= ,且()j j s x y = (0,1,)j n ⋯= （1.1）成立，则称 ()s x 为三次样条差值函数。

从定义知，要求出()s x ，在每个应小区间1[,]j j x x + 上确定4个待定系数，共有 n 个小区间，故应确定4n 个参数，根据()s x 在[,]a b 上二阶导数连续，在节点()1,2,3,,1j x j n ⋯=-处应满足连续性条件(0)(0),j j s x s x -=+ ''(0)(0),j j s x s x -=+''''(0)(0)j j s x s x -=+ （1.2） 共有 3n-3个条件，再加上()s x 满足插值条件（1.1），共有4n-2个条件，因此还需要2个条件才能确定()s x 。

通常可在区间[,]a b 端点0,n a x b x ==上各加一个条件（称边界条件），边界条件可根据实际的问题要求给定。

2012-2013（1）专业课程实践论文三次样条插值樊旭，0818180212，R数学08-2班当插值节点很多时，插值多项式的次数就会很高，这不仅增大了计算量，还会影响结果的精确度。

虽然可以采用原始的插值法，但是主要缺点就是分段接头处不光滑，导数不连续。

因此构想了这样的函数，既能分段的低次插值，又能保证接头处光滑，就产生了三次样条插值。

三次样条插值应该满足的条件：1：满足定义2：S(f x -0)=S(f x +0), S '(f x -0)= S '(f x +0),S ''( f x -0)= S ''(f x +0) 3：边界条件：(1) S '(n x )=0f 'S '(n x )=n f '(2) S ''(0x )=0f 'S ''(n x )=n f '(3) S ''(0x )=S ''(n x )=04：周期函数时：S(0x +0)=S(n x -0),S '(0x +0)=S '(n x -0),S ''(0x +0)=S ''(n x -0),计算方程：1：三转角方程：S(x)=[][][][]11221112212233)()()(2)()(2)(++++++--+--+-+-+-+-j j j j j j j j j j j j j j j j j j m h x x x x m h x x x x y h x x h x x y h x x h x x2：三弯矩方程：S(x)=()j j j j i j j j j i j j J j j j h x x h M y h x x h M y h x x M h x x M -⎪⎭⎫ ⎝⎛-+-⎪⎭⎫ ⎝⎛-+-+-+++++666)(6221113131 以1为例用VC++程序编程就求解NN输入n,xi,yi,和待估点xx ，边界条件。

数值逼近实习报告题目：三次样条插值II班级：计算111班学号：3110811009姓名：刘艳平指导老师：秦新强2013-3题目：三次样条插值（三转角算法） 姓名：刘艳平 学号：3110811009 一、目的意义以上我们讨论的分段线性插值，逼近程度虽然好，但光滑性差，分段三次Hermite 插值，逼近程度好，光滑性也有所提高，但也增加了更多的条件，不太实用。

接下来我们研究的三次样条插值多项式，这就结合二者的优点，即逼近程度虽然好，光滑性强，不需要增加太多的条件，很实用。

当已知区间两个端点的导数值时，既满足边界条件1，用三弯矩算法来求解，相对简单，且很准确。

二、数学公式⎪⎪⎩⎪⎪⎨⎧∈∈∈-],[),(............],[),(],[),(1212101n n nx x x x s x x x x s x x x x s 具体：iii i i ii i iii i i i ii i i i m h x x x x m h x x x x y h x x x x h y h x x x x h x s 22112213211321))(())(())]((2[))]((2[)(--------+--+---+--+=其中：n i x x x x x h i i i i i ,...,2,1],,[,11=∈-=-- i i i i y x s y x s ==--)(,)(11需要求出：ii i i m x s m xs ='='--)(,)(11三、计算流程 Step 1:输入节点n x x x ,...,,10，函数值ny y y ,...,,10,边界条件及x.Step 2:计算⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+=-=+=-++++]),[],[(3111111i i i i i i i ii i ii i i i i x x f x x f g h h h h h h λμλμλStep 3：根据边界条件，求解相应的方程得到nm m m ,...,,10.Step 4：判断x 属于],[1i i x x -，i=1,2,...,n中的那一个.Step 5：计算)(x s y i ≈.Step 6：输出y. 四、代码#include <iostream.h> //三次样条三转角边界条件1 #include "process.h" #define N 2 #define M 4double a[N],d[N],c[N],p[N],q[N],b[N],xx[N]; void p_q() { int k; if(d[0]==0) {cout<<"Method failed"<<endl;exit(0);}p[0]=d[0];q[0]=c[0]/p[0];for(k=1;k<N-1;k++){p[k]=d[k]-a[k]*q[k-1];if(p[k]==0){cout<<"Method failed"<<endl; exit(0);}q[k]=c[k]/p[k];}p[N-1]=d[N-1]-a[N-1]*q[N-2]; if(p[N-1]==0){cout<<"Method failed"<<endl; exit(0);}}void jisuan(){int k;double y[N];y[0]=b[0]/p[0];for(k=1;k<=N-1;k++)y[k]=(b[k]-a[k]*y[k-1])/p[k];xx[N-1]=y[N-1];for(k=N-2;k>=0;k--)xx[k]=y[k]-q[k]*xx[k+1];for(k=0;k<N;k++)cout<<"m的值为:"<<xx[k]<<'\t'<<endl;}void main(){double X[M],Y[M],h[M],v[M],u[M],m[M],g[M],mm[M]; int i;double x,y1,y2,w1,w2,w3,w4,w,hh;cout<<"输入已知点的X坐标:"<<endl;for(i=0;i<M;i++)cin>>X[i];cout<<"输入已知点的Y坐标:"<<endl;for(i=0;i<M;i++)cin>>Y[i];cout<<"输入边界条件:"<<endl;cin>>y1>>y2;cout<<"输入需要计算的x值:"<<endl;cin>>x;for(i=1;i<M;i++){h[i]=X[i]-X[i-1];}for(i=1;i<M-1;i++){v[i]=h[i+1]/(h[i]+h[i+1]);u[i]=h[i]/(h[i]+h[i+1]);g[i]=3*((v[i]*(Y[i+1]-Y[i])/(X[i+1]-X[i]))+(u[i]*(Y[i]-Y[i-1])/(X[i]-X[i-1])));}a[0]=0;b[0]=g[1]-v[1]*y1;for(i=1;i<N;i++)a[i]=v[i+1];for(i=0;i<N;i++)d[i]=2;for(i=0;i<N;i++)c[i]=u[i+1];for(i=1;i<N;i++)b[i]=g[i+1];p_q();jisuan();mm[0]=y1;mm[1]=xx[0];mm[2]=xx[1];mm[3]=y2;for(i=1;i<M;i++){hh=X[i]-X[i-1];if(x>=X[i-1]){w1=((hh+(x-X[i-1])*2)*(x-X[i])*(x-X[i])*Y[i-1])/(hh*hh*hh);w2=((hh-(x-X[i])*2)*(x-X[i-1])*(x-X[i-1])*Y[i])/(hh*hh*hh);w3=((x-X[i-1])*(x-X[i])*(x-X[i])*mm[i-1])/(hh*hh);w4=((x-X[i])*(x-X[i-1])*(x-X[i-1])*mm[i])/(hh*hh);w=w1+w2+w3+w4;}}cout<<"x 点的结果为"<<w<<endl; }五、数值实例已知函数y=f(x),的如下数据，试求其在区间]3,0[上的三次样条插值函数S(x). X 0 1 2 3 Y0 1 0 1 y '1解：在此题中，解得：151,15421=-=m m其中：0,130==m m ，一起带入上边公式即可得到所求的三次样条函数如下⎪⎪⎪⎩⎪⎪⎪⎨⎧∈--+---=∈--+------=∈---+--=]3,2[,)3)(2(151)2)](3(21[)(]2,0[),2()1(151)2)(1(154)2)](1(21[)(]1,0[),1(154)1()]1(21[)(22322222221x x x x x x s x x x x x x x x s x x x x x x x x s 假设，求出x=1.5处的函数值为：0.458333 程序运行结果截图如下：六、对计算机结果进行分析评价三次样条插值函数与三次Hermite插值函数相比，不仅光滑度有提高，而且要求求解时还不需要增加内节点处的导数值，因此比较实用。

非常类似前面的三弯矩法，这里的sanzhj函数和intersanzhj作用相当于前面的sanwanj和intersanwj，追赶法程序通用，代码如下。

%%%%%%%%%%%%%%%%%%%function [newu,w,newv,d]=sanzhj(x,y,x0,y0,y1a,y1b)% 三转角样条插值% 将插值点分两次输入，x0 y0 单独输入% 边值条件a的一阶导数 y1a 和b的一阶导数 y1bn=length(x);m=length(y);if m~=nerror('x or y 输入有误，再来');endv=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1);w=2*ones(n-1,1);h0=x(1)-x0;h=zeros(n-1,1);for k=1:n-1h(k)=x(k+1)-x(k);endv(1)=h0/(h0+h(1));u(1)=1-v(1);d(1)=3*(v(1)*(y(2)-y(1))/h(1)+u(1)*((y(1)-y0))/h0);%for k=2:n-1v(k)=h(k-1)/(h(k-1)+h(k));u(k)=1-v(k);d(k)=3*(v(k)*(y(k+1)-y(k))/h(k)+u(k)*(y(k)-y(k-1))/h(k-1));endd(1)=d(1)-u(1)*y1a;d(n-1)=d(n-1)-v(n-1)*y1b;newv=v(1:n-2,:);newu=u(2:n-1,:);%%%%%%%%%%%%function intersanzhj(x,y,x0,y0,y1a,y1b)% 三转角样条插值%第一部分n=length(x);m=length(y);if m~=nerror('x or y 输入有误，再来');end%重新定义hh=zeros(n,1);h(1)=x(1)-x0;for k=2:nh(k)=x(k)-x(k-1);end% 调用三转角函数[a,b,c,d]=sanzhj(x,y,x0,y0,y1a,y1b);% 三对角方程m=chase(a,b,c,d);% 求MM=[1;m;0];% 求插值函数fprintf('三次样条（三转角）插值的函数表达式\n');syms X ;fprintf('S0--1:\n');S(1)=collect((y0/h(1).^3)*(X-x(1)).^2*(h(1)+2*(X-x0))+(y(1)/h(1).^3)*(X-x0).^2*(h(1)+2*(x(1)-X))+(M(1)/h(1).^2)*(X-x0)*(X-x(1)).^2+(M(2)/h(1).^2)*(X-x(1))*(X-x0).^2);for k=2:nfprintf('S%d--%d:\n',k-1,k);S(k)=collect((y(k-1)/h(k).^3)*(X-x(k)).^2*(h(k)+2*(X-x(k-1)))+(y(k)/h(k).^3)*(X-x(k-1)).^2*(h(k)+2*(x(k)-X))+(M(k)/h(k).^2)*(X-x(k-1))*(X-x(k)).^2+(M(k+1)/h(k).^2)*(X-x(k))*(X-x(k-1)).^2);endS=S.';disp(S);fprintf('以上为样条函数（三转角）解析式，显示为手写习惯如下：\n');pretty(S);%第二部分%是否继续运行程序myloop=input('继续运行程序输入“1”，否则输入“0”\n');if myloopwhile myloopxi=input('输入需要计算的点的值，并按回车键\n');if xi>x0|xi<x(n)fprintf('现在开始计算输入点的插值函数值……\n');elsefprintf('输入数值不在插值范围内，请重新输入\n');xi=input('输入需要计算的点的值，并按回车键……\n'); end% 确定输入的数值应该使用哪个解析式newx=[x0;x];[r,suoy]=min(abs(newx-xi));fprintf('输入点的插值函数值为:\n\n');fprintf('\t');if xi<=newx(suoy)f=subs(S(suoy-1),X,xi);elsef=subs(S(suoy),X,xi);enddisp(f);myloop=input('继续计算输入“1”，终止计算输入“0”\n'); endelsereturn;end。

三次样条插值算法详解下面详细介绍三次样条插值方法的具体步骤：1.数据准备：首先，需要获得一组数据点，这些数据点包含了所需插值曲线的关键信息。

通常情况下，数据点是从实际观测中获得的。

2.区间划分：将插值区间划分为若干个小区间，每个小区间对应一个三次函数。

3. 函数构建：对于每个小区间，在该区间内构建一个三次函数。

这里使用三次多项式进行构建，形如 f(x) = ax^3 + bx^2 + cx + d。

每个小区间内的函数有四个待定系数：a、b、c、d。

4.条件设置：为了确定每个小区间内的函数，需要设置相应的条件。

一般来说，需要满足以下两个条件：(a) 函数值条件：保证每条小区间内的函数值通过对应的数据点。

即，对于每个小区间，函数值满足 f(xi) = yi，其中(xi, yi)表示第i个数据点的横纵坐标。

(b)导数条件：保证每个小区间内函数的导数连续。

这可以通过限制每个小区间内的函数的一阶导数（即斜率）相等来实现。

5.矩阵方程求解：根据上述条件设置，可以得到一个线性方程组，其中待求的系数为未知数。

将上述条件代入方程组中，然后求解该方程组以获得每个小区间内的函数系数。

6.曲线绘制：通过得到的函数系数，可以计算每个小区间内的函数值，并连接这些函数值，最终得到整个插值曲线。

三次样条插值方法是一种非常强大和灵活的插值方法，适用于各种类型的数据点，包括均匀和非均匀间距的数据。

通过调整划分区间的个数，可以控制插值曲线的光滑程度。

一般来说，插值区间越多，插值曲线越平滑，但对输入数据的噪声更敏感。

总结起来，三次样条插值是一种高级的插值方法，通过构建三次函数来逼近数据点，可以产生平滑的插值曲线。

它的基本思想是将插值区间划分为若干个小区间，并在每个小区间内构建一个三次函数。

通过设置函数值条件和导数条件，可以得到一个线性方程组，从而求解出每个小区间内的函数系数。

最终连接每个小区间内的函数值，得到整个插值曲线。