三次样条插值自然边界条件
三次样条插值函数
三次样条插值函数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 ----------⨯⎛⎫⎪ ⎪ ⎪ ⎪⎪⎪ ⎪-- ⎪-+- ⎪--⎝⎭为了使样条插值问题有惟一解，我们在原有方程基础上增加两个边界条件。
详细讲解三次样条插值法及其实现方法
样条函数的定义 定义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次样条曲线插值的算法边界条件取为自然
用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)为已知的点的坐标。
第5章-3-三次样条插值PPT课件
(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(-∞,+∞)
样条函数及三次样条插值PPT课件
(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 插值函数相比，不仅光滑度有提高，而且要求求解时还不需要增加内节点处的导数值，因此比较实用。
第五章(3)三次样条插值
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
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
例：已知一组数据点，编写一程序求解三次样条插值函数满足
并针对下面一组具体实验数据
求解，其中边界条件为.
1）三次样条插值自然边界条件源程序：
function s=spline3(x,y,dy1,dyn)
%x为节点，y为节点函数值，dy1，dyn分别为x=0.25,0.53处的二阶导
m=length(x);n=length(y);
if m~=n
error('x or y输入有误')
return
end
h=zeros(1,n-1);
h(n-1)=x(n)-x(n-1);
for k=1:n-2
h(k)=x(k+1)-x(k);
v(k)=h(k+1)/(h(k+1)+h(k));
u(k)=1-v(k);
end
g(1)=3*(y(2)-y(1))/h(1)-h(1)/2*dy1;
g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)/2*dyn;
for i=2:n-1
g(i)=3*(u(i-1)*(y(i+1)-y(i))/h(i)+v(i-1)*(y(i)-y(i-1))/h(i-1)); end
for i=2:n-1;
A(i,i-1)=v(i-1);
A(i,i+1)=u(i-1);
end
A(n,n-1)=1;
A(1,2)=1;
A=A+2*eye(n);
M=zhuigf(A,g); %调用函数，追赶法求M
fprintf('三次样条（三对角）插值的函数表达式\n');
syms X;
for k=1:n-1
fprintf('S%d--%d:\n',k,k+1);
s(k)=(h(k)+2*(X-x(k)))./h(k).^3.*(X-x(k+1)).^2.*y(k)...
+(h(k)-2*(X-x(k+1)))./h(k).^3.*(X-x(k)).^2.*y(k+1)... +(X-x(k)).*(X-x(k+1)).^2./h(k).^2*M(k)+(X-x(k+1)).*... (X-x(k)).^2./h(k).^2*M(k+1);
end
s=s.';
s=vpa(s,4);
%画三次样条插值函数图像
for i=1:n-1
X=x(i):0.01:x(i+1);
st=(h(i)+2*(X-x(i)))./(h(i)^3).*(X-x(i+1)).^2.*y(i)...
+(h(i)-2.*(X-x(i+1)))./(h(i)^3).*(X-x(i)).^2.*y(i+1)...
+(X-x(i)).*(X-x(i+1)).^2./h(i)^2*M(i)+(X-x(i+1)).*...
(X-x(i)).^2./h(i)^2*M(i+1);
plot(x,y,'o',X,st);
hold on
End
plot(x,y);
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%调用的函数：
%追赶法
function M=zhuigf(A,g)
n=length(A);
L=eye(n);
U=zeros(n);
for i=1:n-1
U(i,i+1)=A(i,i+1);
end
U(1,1)=A(1,1);
for i=2:n
L(i,i-1)=A(i,i-1)/U(i-1,i-1);
U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);
end
Y(1)=g(1);
for i=2:n
Y(i)=g(i)-L(i,i-1)*Y(i-1);
end
M(n)=Y(n)/U(n,n);
for i=n-1:-1:1
M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i);
end
2)在命令窗口输入x，y,dy1,dyn,得到三次样条函数：
>>x=[0.25,0.3,0.39,0.45,0.53];
>>y=[0.5,0.5477,0.6245,0.6708,0.7280];
>>dy1=0;
>>dyn=0;
>>s=spline3(x,y,dy1,dyn)
运行结果：
三次样条（三对角）插值的函数表达式
S1--2:
S2--3:
S3--4:
S4--5:
.5000*(-3600.+.1600e5*X)*(X-.3000)^2+.5477*(5200.-.1600e5*X)*(X-.2500)^2+394.8*(X-.2500)*(X-.3000)^2+355.2*(X-.3000)*(X-.2500)^2
.5477*(-699.6+2743.*X)*(X-.3900)^2+.6245*(1193.-2743.*X)*(X-.3000)^2+109.6*(X-.300 0)*(X-.3900)^2+96.77*(X-.3900)*(X-.3000)^2
.6245*(-3333.+9259.*X)*(X-.4500)^2+.6708*(4444.-9259.*X)*(X-.3900)^2+217.7*(X-.39 00)*(X-.4500)^2+207.6*(X-.4500)*(X-.3900)^2
.6708*(-1602.+3906.*X)*(X-.5300)^2+.7280*(2227.-3906.*X)*(X-.4500)^2+116.8*(X-.45 00)*(X-.5300)^2+109.2*(X-.5300)*(X-.4500)^2
如将三次样条函数加以整理，可用如下程序：
s=collect(s);
则输出结果为：
s=
.4595000000000-13.200*X^3+9.9000000*X^2-1.48800000000*X
.37459306800000-4.2924*X^3+3.66332200*X^2-.137976090000*X
.5180681700000-3.3917*X^3+3.90576600*X^2-.733130370000*X
-.408614400000e-1+2.5768*X^3-4.03188800*X^2+2.868770320000*X。