当前位置：文档之家 › 三次样条插值自然边界条件

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

  • 格式：doc
  • 大小：60.50 KB
  • 文档页数：3

下载文档原格式

  / 3

合集下载

三次样条插值函数

三次样条插值函数

三次样条插值函数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次样条曲线插值的算法边界条件取为自然

用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课件

第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课件

样条函数及三次样条插值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)三次样条插值

第五章(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

第二章三次样条插值

第二章三次样条插值
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 在所有插值点的导数值。

2[1].7三次样条插值

2[1].7三次样条插值
确定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 )

C#三次样条插值函数(自然边界)

C#三次样条插值函数(自然边界)

/****************函数说明*********************///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];//三转角法的待定一阶系数法。

【VIP专享】_三次样条插值课件

【VIP专享】_三次样条插值课件
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,'.')

三次样条插值 python 边界条件

三次样条插值 python 边界条件

三次样条插值 python 边界条件三次样条插值是一种常见的数据拟合方法，适用于任意分段连续函数的插值。

在 Python 中，可以使用 SciPy 库中的 interpolate 模块来进行三次样条插值。

不过，在进行样条插值时，需要指定边界条件，以确定插值的唯一性。

常见的边界条件包括自然边界条件、固定边界条件和周期边界条件。

本文将详细介绍三次样条插值在Python 中的边界条件设置方法。

一、自然边界条件自然边界条件是指插值函数的二阶导数在插值区间两端均为零。

在 SciPy 中，可以使用 interp1d 函数来实现自然边界条件的三次样条插值。

具体实现方法如下：```pythonfrom scipy.interpolate import interp1d# x和y为插值区间的数据点f = interp1d(x, y, kind='cubic',fill_value='extrapolate')```其中，kind 参数设置为 cubic 表示三次样条插值；fill_value 参数设置为 extrapolate 表示在插值区间之外的点依然可以进行插值。

二、固定边界条件固定边界条件是指插值函数在插值区间两端的值与一定值相等。

在 SciPy 中，可以使用 CubicSpline 类来实现固定边界条件的三次样条插值。

具体实现方法如下：```pythonfrom scipy.interpolate import CubicSpline# x和y为插值区间的数据点，left和right为左右端点的值f = CubicSpline(x, y, bc_type='clamped', left=left, right=right)```其中，bc_type 参数设置为 clamped 表示固定边界条件；left 和 right 分别为左右端点的值。

三、周期边界条件周期边界条件是指插值函数在插值区间两端的值相等。

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

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

例：已知一组数据点，编写一程序求解三次样条插值函数满足并针对下面一组具体实验数据求解，其中边界条件为.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~=nerror('x or y输入有误')returnendh=zeros(1,n-1);h(n-1)=x(n)-x(n-1);for k=1:n-2h(k)=x(k+1)-x(k);v(k)=h(k+1)/(h(k+1)+h(k));u(k)=1-v(k);endg(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-1g(i)=3*(u(i-1)*(y(i+1)-y(i))/h(i)+v(i-1)*(y(i)-y(i-1))/h(i-1)); endfor i=2:n-1;A(i,i-1)=v(i-1);A(i,i+1)=u(i-1);endA(n,n-1)=1;A(1,2)=1;A=A+2*eye(n);M=zhuigf(A,g); %调用函数，追赶法求Mfprintf('三次样条（三对角）插值的函数表达式\n');syms X;for k=1:n-1fprintf('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);ends=s.';s=vpa(s,4);%画三次样条插值函数图像for i=1:n-1X=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 onEndplot(x,y);grid on%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%调用的函数：%追赶法function M=zhuigf(A,g)n=length(A);L=eye(n);U=zeros(n);for i=1:n-1U(i,i+1)=A(i,i+1);endU(1,1)=A(1,1);for i=2:nL(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);endY(1)=g(1);for i=2:nY(i)=g(i)-L(i,i-1)*Y(i-1);endM(n)=Y(n)/U(n,n);for i=n-1:-1:1M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i);end2)在命令窗口输入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。

最新三次样条插值的Matlab实现(自然边界和第一边界条件)

最新三次样条插值的Matlab实现(自然边界和第一边界条件)

精品资料三次样条插值的M a t l a b实现(自然边界和第一边界条件)........................................（第一边界条件）源代码：function y=yt1(x0,y0,f_0,f_n,x) _____________(1)%第一类边界条件下三次样条插值；%xi 所求点；%yi 所求点函数值；%x 已知插值点；%y 已知插值点函数值；%f_0左端点一次导数值；%f_n右端点一次导数值；n = length(x0);z = length(y0);h = zeros(n-1,1);k=zeros(n-2,1);l=zeros(n-2,1);S=2*eye(n);for i=1:n-1h(i)= x0(i+1)-x0(i);endfor i=1:n-2k(i)= h(i+1)/(h(i+1)+h(i)); l(i)= 1-k(i);end%对于第一种边界条件：k = [1;k];_______________________(2)l = [l;1];_______________________(3)%构建系数矩阵S：for i = 1:n-1S(i,i+1) = k(i);S(i+1,i) = l(i);end%建立均差表：F=zeros(n-1,2);for i = 1:n-1F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i)); endD = zeros(n-2,1);for i = 1:n-2F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i)); D(i,1) = 6 * F(i,2);end%构建函数D：d0 = 6*(F(1,2)-f_0)/h(1);___________(4)dn = 6*(f_n-F(n-1,2))/h(n-1);___________(5)D = [d0;D;dn];______________(6)m= S\D;%寻找x所在位置，并求出对应插值：for i = 1:length(x)for j = 1:n-1if (x(i)<=x0(j+1))&(x(i)>=x0(j)) y(i) =( m(j)*(x0(j+1)-x(i))^3)/(6*h(j))+...(m(j+1)*(x(i)-x0(j))^3)/(6*h(j))+...(y0(j)-(m(j)*h(j)^2)/6)*(x0(j+1)-x(i))/h(j)+... (y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j) ; break;else continue;endendend（2）（自然边界条件）源代码：仅仅需要对上面部分标注的位置做如下修改: __(1):function y=yt2(x0,y0,x)__(2):k=[0;k]__(3):l=[l;0]__(4)+(5):删除—(6):D=[0:D:0]。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 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。