三次样条函数的自动求法摘要：在MATLAB 的The Spline Toolbox 中，没有给出三次样条函数表达的求法，可在教学过程中，或在实际问题中，我们需要知道样条插值函数的分段表达式。

在现行数值分析教材中，一般都是通过解方程确定三次样条插值函数的表达式，但这种方法的工作量很大。

在本文中，我们用MATLAB 语言编制了三个程序，给出在三种边界条件下，三次样条插值函数表达式的自动求法。

关键词：三次样条函数；边界条件；插值0 引言分段低次插值多项式具有计算简单、收敛性有保证、数值稳定等优点，但它不能保证整条曲线的光滑性甚至连续性，从而不能满足一些工程技术的要求。

从20 世纪60 年代开始，由航空、造船等工程设计的需要而发展起来的样条插值方法，既保留了分段低次插值多项式的各种优点，又保证了插值函数的光滑性，已在许多领域里得到越来越广泛的应用。

在教学过程中，或在实际问题中，我们需要知道样条插值函数的分段表达式。

可在MATLAB 的The Spline Toolbox 中，没有直接给出三次样条函数表达式的求法，在现行数值分析教材中，一般都是在给定条件下，通过解方程而确定三次样条插值函数的表达式，尽管在计算过程中可借助数学软件来完成，但这种方法的工作量仍然很大。

本文中，利用数学软件MATLAB ，我们给出了三次样条插值函数表达式的自动求法，这样不但解决了上述问题，而且给出了用数学软件解决实际问题的一个范例。

1 计算方法定义对于给定的函数值),,1,0)((n k x f y k k ==其中b x x x a n =<<<= 10，如果函数)(x S 满足条件：（1）)(x S 在每个子区间[k k x x ,1-](k=1,2,n , )上都是不高于三次的多项式；（2）)(x S 、)(x S '、)(x S ''在[a,b]上都连续；（3）),2,1,0()(n k y x S k ==。

则称)(x S 为函数)(x f 关于节点n x x x ,,,10 的三次样条插值函数。

要求三次样条插值函数)(x S ，只需在每个子区间[k k x x ,1-]上确定一个三次多项式k k k k k d x c x b x a x S +++=23)()(x S k 共有4个系数，确定它们需要4个条件，因此要完全确定)(x S 共需4n 个条件。

由)(x S 所满足的条件(1)、(2)、(3），可确定4n-2个条件，还缺少两个条件。

这两个条件通常由实际问题对三次样条插值函数在端点的状态要求给出，称之为边界条件，常用的边界条件有以下三类。

第一类边界条件：给定端点处的一阶导数值)(y x S '='，n n y x S '=')( 第二类边界条件：给定端点处的二阶导数值)(y x S ''=''，n n y x S ''='')( 当00=''y ，0=''n y 时，)(x S 在端点处不受力，呈自然状态，故称之为自然边界条件。

第三类边界条件是周期性条件：如果)(x f y =是以b-a 为周期的函数，)(x S 也应是以b-a 为周期的函数，于是)(x S 在端点处满足条件)0()0(00-'=+'x S x S)0()0(-''=+''n n x S x S设),,1,0()(n k M x S k k ==''，记),,2,1(1n k x x h k k k =-=-，则可用k M 表示)(x S k ：kk k k k k k k k k k k k k k k k h x x h M y h x x h M y h x x M h x x M x S 122113131)6()6(6)(6)()(-------+--+-+-=),,2,1],,[(1n k x x x k k =∈- （1）若记⎪⎪⎪⎩⎪⎪⎪⎨⎧---+=-=+=+=-++++++)(611111111k k k k k k K k k k K k k k k k k k h y y h y y h h g h h h h h h μλμ （2） )1,,2,1(-=n k在第一类边界条件下，可得确定0M 、1M 、……、n M 的线性方程⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----n n n n n n g g g g M M M M 1101101111......212............212λμλμ （3） 其中⎪⎪⎩⎪⎪⎨⎧--'='--=-)(6)(61010110n n n n n n h y y y h g y h y y h g 在第二类边界条件下，可得确定0M 、1M 、……、n M 的线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡''-''-=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------n n n n n n n n n y g g g y g M M M M 11220111221122221......22............22λμμλμλμλ （4） 其中：00y M ''=，n n y M ''=。

在第三类边界条件下，可得确定1M 、……、0M M n =的线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----n n n n n n n n g g g g M M M M 121121112211......22............22μλλμλμμλ （5）其中⎪⎪⎪⎩⎪⎪⎪⎨⎧---+=-=+=+=-)(6111011111n n n n n n n n n n n h y y h y y h h g h h h h h h μλμ 2 主要结果在MATLAB 的The Spline Toolbox 中，没有直接给出三次样条函数表达式的求法，只给出其向量表示形式。

根据（1）~（5）式，我们用MATLAB 语言编制了三个小程序，给出了在常用的三种边界条件下的三次样条插值函数表达式的自动求法。

命题 三种边界条件下的三次样条插值函数可分别由以下程序自动求得。

源程序 1 function S=spf1(x,y,xx,m)%1 此函数为三次样条插值函数，计算满足边界条件一时的插值函数（此时需要定义一个系统变量’xx ’，即 syms xx ），函数的系数输出形式为分数，如果系数要以小数形式输出，请给出精度m 。

%2 如果给出的变量xx 不是系统变量，而是一组数值，则计算满足边界条件一时在xx 处的插值。

%3 输入时请将边界条件（即插值函数在两端点的一阶导数）%放在向量y 的第一个分量和最后一个分量。

%4 如果不输入边界条件，默认插值函数在两端点的一阶导数均为零。

n=length(x);If n<2,error("There should be at least two data points."),end If any(diff(x)<0),[x,ind]=sort(x);else,ind=1:n;endx=x(:);dx=diff(x);If all(dx)==0,error("The data abscissae should be distinct."),end[yd yn]=size(y);%if Y happen to be a column % matrix,change it to the expected row matrix.If yn==1,yn=yd;y=reshape(y,1,yn);yd=1;endIf yn==nnotaknot=1;elseif yn==n+2notaknot=0;endslopes=y(:,[1 n+2]).':y(:,[1 n+2])=[];ElseError('Abscissa and ordinate vector should be of the same length.') Endyi=y(:.ind) ;dd=ones(1,yd);Dx=diff(x);divdif=diff(yi)/dx(:,dd);a=zeros(1,n-1);c=a;c=zeros(1,n-1);a(1:n-2)=dx(1:n-2)/(dx(1:n-2)+dx(2:n-1));c(2:n-1)=dx(2:n-1)/(dx(1:n-2)+dx(2:n-1));d(2:n-1)=6*(divdif(2:n-1)-divdif(1:n-2))/(dx(1:n-2)+dx(2:n-1));a(n-1)=1;c(1)=1;If notaknotd(1)=6*(yi(2)-yi(1))/(dx(1)^2);d(n)=6*(yi(n-1)-yi(n))/(dx(n-1)^2);Elsed(1)=6*((yi(2)-yi(1))/dx(1)-endslopes(1))/dx(1);d(n)=6*(endslopes(2)-(yi(n)-yi(n-1))/dx(n-1))/dx(n-1);Endb(1:n)=2;d=tridi(a,b,c,d);If isnumeric(xx)==1m=length(xx);pp=0;For k=1:mFor i=1:n-1If(xx(k)<=x(i+1))&(xx(k)>=x(i))pp(k)=d(i)*(x(i+1)-xx(k))^3/(6*dx(i))+d(i+1)*(xx(k)-x(i))^3/(6*dx(i)) +(y(i)-d(i)*dx(i)^2/6)*(x(i+1)-xx(k))/dx(i)+(y(i+1)-d(i+1)*dx(i)^2/6) *（xx(k)-x(i)）/dx(i);EndEndS=ppEndElseif(isnumeric(xx)+1==1)&(nargin==3)For i=1:n-1Pp(i)=d(i)*(x(i+1)-xx)^3/(6*dx(i))+d(i+1)*(xx(k)-x(i))^3/(6*dx(i))+(y (i)-(d(i)^2)/6)*(x(i+1)-xx(k))/dx(i)+(y(i+1)-(d(i+1)*dx(i)^2)/6)*(xx( k)-x(i))/dx(i);Pp(i)=simplify(pp(i));Fprintf(‘in[%g,%g]\n’,x(i),x(i+1));Fprintf(‘S(%d)=’,i);Pretty(pp(i));EndElseDigits(m);For i=1:n-1Pp(i)=d(i)*(x(i+1)-xx)^3/(6*dx(i))+d(i+1)*(xx-x(i))^3/(6*dx(i))+(y(i) -(d(i)^2)/6)*(x(i+1)-xx)/dx(i)+(y(i+1)-(d(i+1)*dx(i)^2)/6)*(xx-x(i))/ dx(i);Sym2poly(pp(i));Ans=sym(ans,’d’);Poly2sym(ans,xx);Fprintf(‘in[%g,%g]\n’,x(i),x（i+1）);Fprintf(‘S(%d)=’,i);Pretty(ans);EndEnd3 应用实例例 给定函数值k x =1,2,4,5,k y =f(k x )=1,3,4,2,求满足边界条件'S （1）=0，'S （5）=1的三次样条插值函数。