数值分析课程设计报告书三次样条插值的三弯矩法.docx

数值分析课程设计报告书院系名称:学生姓名:专业名称:班级:时间:实验一三次样条插值的三弯矩法一、实验目的已知数据X , y = f (X ), i = 0, , n及边界条件y G)(x ), j = 0,1 n，求f (x)的i i ij j三次样条插值函数S(x).要求输出用追赶法解出的弯矩向量M = M0, ,Mn及 S ( )(t), i = 0, , m, k = 0,1,2的值.画出y = S (x)的图形，图形中描出插值点(x, y.) 及(t,S(t)分别用o和* 标记.ii 二、实验原理1.用追赶法求解第二类边界条件的三弯矩方程:一 2R112人1 M 一 M:=6.R.2 n-1.12MMn 1nhhj-1j-+ hj,人jh=jh + h,h = x',一x.其中日= jf x0, x0,七f x , x , x 012f x :, x , x n 一 2， n-1 nf x , x , x n-1n2.得出样条函数表达式:(x- x)3(x - x )3 M h2 x - xS(x) = M + M二 + (y -jj)+ (yj6hj+1 6h j 6 hj+13.计算 S(k)(t), i = 0, , m, k = 0,1,2 .i三、实验结果所用数据: x=-2223,-1987,-18465,-1292,-12266,-11056,-08662,-06594,-02671,-00452,0.5385,12564,14398,15415,17646,19678,2236;y=0.83995,1.1696,1.3141,1.6992,1.7312,1.7847,1.8708,1.9262,1.9881,1.9997,1.9511,1.7169,1.618,1.5543,1.3871,1.191,0.81662;d2s1= -4.5000;d2sn= -4.8967; %第二种边界条件t=-2.223,-1.9443,-1.6656,-1.3869,-1.1083,-0.82956,-0.55088,-0.27219,0.0065,0.28519,0.56387,0.84256,1.1212,1.3999,1.6786,2.236;；(指定计算点)计算结果：用追赶法求得的弯矩量为要计算的16个节点出的值为M0=-4.5000S(-2.223)=0.83995M1=-1.9282S(-1.9443)=1.2163M2=-0.5336S(-1.6656)=1.4617M3=-0.8000S(-1.3869)=1.6485M4=-0.4188S(-1.1083)=1.7836M5=-0.4663S(-0.82956)=1.8818M6=-0.4264S(-0.55088)=1.9487M7=-0.3343S(-0.27219)=1.9876M8=-0.3677S(0.0065)=2M9=-0.3113S(0.28519)=1.9865M10=-0.3705S(0.56387)=1.9462M11=-0.4185S(0.84256)=1.8767M12=-0.6533S(1.1212)=1.7767M13=-0.7253S(1.3999)=1.6411M14=-1.0132S(1.6786)=1.4567M15=-1.2968S(2.236)=0.81662M16=-4.8967要计算的16个节点出的一阶导数值为：要计算的16个节点出的二阶导数值为：S1(-2.223)=1.8267S2(-2.223)=-4.5S1(-1.9443)=1.058S2(-1.9443)=-1.5044S1(-1.6656)=0.76267S2(-1.6656)=-0.62053S1(-1.3869)=0.57108S2(-1.3869)=-0.75438S1(-1.1083)=0.41615S2(-1.1083)=-0.46529S1(-0.82956)=0.29348S2(-0.82956)=-0.41006S1(-0.55088)=0.18877S2(-0.55088)=-0.34356S1(-0.27219)=0.089715S2(-0.27219)=-0.36728S1(0.0065)=-0.0028801S2(0.0065)=-0.31653S1(0.28519)=-0.095035S2(0.28519)=-0.34481S1(0.56387)=-0.19691S2(0.56387)=-0.37221S1(0.84256)=-0.30324S2(0.84256)=-0.39084S1(1.1212)=-0.41474S2(1.1212)=-0.40946S1(1.3999)=-0.56694S2(1.3999)=-0.60221S1(1.6786)=-0.76939S2(1.6786)=-0.90224S1(2.236)=-1.8916S2(2.236)=-4.8967画出的图形为:四、实验分析通过实验结果我们，知道三弯矩法求出满足初始条件的三次样条函数，与其 他插值函数的构造相比，三次样条插值法的计算量要小得多。其运行速度比较快 结果精确，在解决大量数据时可以采用，在工程，计算领域比较广泛。实验二基于正交多项式的最小二乘法拟合一、实验目的编制以离散点花=0,1,2,m）的正交多项式P G）为基的最小二乘拟合 程序,并用于对下列数据做次多项式最小二乘拟合.取权3 = 0 ,求出拟合曲线y = S G ）=Y a * P （x）=咒a xk，输出0mk kkk=0k=0a *, P （x）（系数），a（k = 0,1，2,3n）, S（t ）, i = 0,1m 的值及平方误差 |6 2，画 k k- ki112出y = S（x）的图形，并描出插值点（气,七）及（,S（”分别用o和*标记.二、实验原理根据给定节点,气,x及权函数W(X)> 0，造出带权3 G)正交的多项式p (x), (n < m)有递推公式n以LG=)GGpoplp”< <G1-p. p)qb.p o气 = )=1,2,n 1)£ 3(Xi )xP2 I.)£3(X )P2 (x )=1,n -1)-4=0£(x )P2 (x )ii =0仁)£(x )f (x )p (x )(k = 1,n).(f, P )i i k i并计算系数a *k= -= -4=0，时£ 成)P 2 (x )ik ii=0最后可得拟合曲线的函数y= S (x) = a * P (x) + a * P (x) + . + a * P (x) 0 01 1n n三、实验结果所用数据：x=-2.223,-1.987,-1.8465,-1.292,-1.2266,-1.1056,-0.8662,-0.6594,-0.2671,-0.0452,0.5385,1.2564,1.4398,1.5415,1.7646,1.9678,2.236;y=0.83995,1.1696,1.3141,1.6992,1.7312,1.7847,1.8708,1.9262,1.9881,1.9997,1.9511,1.7169,1.618,1.5543,1.3871,1.191,0.81662;d2s1= -4.5000;d2sn= -4.8967; %第二种边界条件t=-2.223,-1.9443,-1.6656,-1.3869,-1.1083,-0.82956,-0.55088,-0.27219,0.0065,0.28519,0.56387,0.84256,1.1212,1.3999,1.6786,2.236;计算结果:拟合系数a *为：平方误差|I；的值：a*(1)=0a*(2)=0.12109a*(3)=-6.3684e-017a*(4)=-0.41667a*(5)=6.3396e-016a*(6)=162.434959e-010拟合函数S3)在t上各点的函数值为：多项式P(i)的升幕系数为：S(0)=-1S(0.11719)=0.76699S(0.23438)=0.93856S(0.35156)=0.37655S(0.46875)=-0.33463S(0.58594)=-0.84526S(0.70313)=-0.99797S(0.82031)=-0.78541S(0.9375)=-0.30763S(1.0547)=0.27023S(1.1719)=0.76031S(1.2891)=0.99456S(1.4063)=0.86722S(1.5234)=0.37779S(1.6406)=-0.32858S(1.7578)=-0.91604S(1.875)=-0.81689P二100000P(2) =-110000P(3) =0.60159 -2 1 0 0 0P(4) =-0.32813 2.3281 -3100P(5) =0.14444 -2.0456 5.0228 -4 1 0P(6) =-0.0634771.6156 -6.3281 8.776 -5 1所画出的图形为:四、实验分析通过结果我们可以看出，当求得的平方误差比较小，通过最小二乘法可以比 较精确地拟合所给的数据；当平方误差比较大时，用最小二乘法拟合数据时就会 出现比较大的误差。实验三龙贝格与高斯求积公式一、实验目的（1）运用龙贝格求积公式计算上述积分I的值,要求到T（0）-T（0） wuch时结束,kk -1输出T表及I的近似值.（2）用5点高斯求积公式及复化3点高斯求积公式计算上述积分,并输出I的近似 值.二、实验原理龙贝格算法：(1) 取 k = 0, h = b - a,求 T(o)= -f (a)+ f (b ) 02令1 k (k记区间la,b的二分次数（2）求梯形值t fy01 2 k ）(k = 1,2,)（3）求加速值T (k) =4 皿 T (k+1)一 1 T (k)m 4 m 1 m-1 4 m 1) m-1（4）若 Tk（o）-Tk1（0）V8，就终止计算。下面是递推公式:T = 1T + 攵罗1 f (x)2 n2 n 2k+1/2k=05点高斯算法bf （x认=工A f （x ）根据课本提供的数据直接写出公式 ak kk = 0复化3点高斯算法与5点高斯算法原理一样三、实验结果所用数据：I = j 1cos(2x)ex3dx ,wuch=1e-9. 0计算结果为:龙贝格算龙贝格算法求得T为：法求的积分值I:0.0656000000-0.2733 -0.386300000-0.3865 -0.4242 -0.42670000-0.4178 -0.4282 -0.4285 -0.4285000-0.4258 -0.4285 -0.4285 -0.4285 -0.428500-0.4279 -0.4285 -0.4285 -0.4285 -0.4285 -0.42850-0.42855-0.4284 -0.4285 -0.4285 -0.4285 -0.4285 -0.4285 -0.4285复化三点五点高斯求的积分值I:高斯求的积分值I:0.428550:-0.430195四、实验分析通过实验，对于一个积分计算，虽然其最后算法是收敛的，但其速度是不同 的；并且其积分值是有点区别的，这是由于每次都会有误差，只是误差的大小不 同，最后的结果也会不同。实验四三次样条与数值积分求导法一、实验目的比较三次样条求导法与数值积分求导法：给定函数f (x)分别利用三次样条 求导法与数值积分求导法计算f K )x = a +生*i,i = 0,1 n二、实验原理1.三次样条求导先用追赶法求解第一类边界条件的三弯矩方程：一 21M0f x0, x0, xH12力1M1-6f x , x , x 0 1 2.H n-12.力n-1M xn-1f x :, x , x n-2 n-1 nL12 _Mnf x , x , x n - 1n nhj，h - x - xh + h jj+1j第一类边界条件：冗1,d f Ex ,x 广)日1,d (f， f lx ,x D00 h 0 10 n n h nn-1 n三次样条一阶导公式f'(x ) M S'(x )- M (xj.! - M '-x? + 七 Lj - Mj+! - Mj h k k j2hj+1 2hh6 j将x值代入公式即可2.数值积分法求导(辛普森求积公式做)4 11 41-1mm 2.1七24 mn-1mn即为所求的一阶导数值三、实验结果所用数据:3 Ef (x)- f (x )- f '(x ) h 2003 Ef (x)- f()h 31.3 f (x )- f (x )hn-1n-33 f (x )- f (x ) f f(x )h nn-2nI = 一1 cos(2x)ex3dx , wuch - 1e - 9 0x. - a + a i, i - 0,1ni n函数在x()处求出的一阶导数真值广为：用三次样条求导法在 x (i )处求出的一阶导数值f 1'为：用数值积分求导法在 xG)处求出的一阶导数值f 2'为：计算结果为：f(0.5)=-29.5562f(0.625)=-12.6798f(0.75)=-6.7443f(0.875)=-4.0956f(1)=-2.7183f(1.125)=-1.9219f(1.25)=-1.4243f(1.375)=-1.0946f(1.5)=-0.86566f(1.625)=-0.70073f(1.75)=-0.57822f(1.875)=-0.48487f(2)=-0.41218f1(0.5)=-29.5562f1(0.625)=-12.4933f1(0.75)=-6.7599f1(0.875)=-4.0826f1(1)=-2.7189f1(1.125)=-1.9207f1(1.25)=-1.4242f1(1.375)=-1.0943f1(1.5)=-0.86581f1(1.625)=-0.69989f1(1.75)=-0.58115f1(1.875)=-0.47915f1(2)=-0.44271f2(0.5)=-29.5562f2(0.625)=-12.4933f2(0.75)=-6.7599f2(0.875)=-4.0826f2(1)=-2.7189f2(1.125)=-1.9207f2(1.25)=-1.4242f2(1.375)=-1.0944f2(1.5)=-0.8656f2(1.625)=-0.70069f2(1.75)=-0.5782f2(1.875)=-0.48485f2(1.75)=-0.41218四、实验分析试验中用三次样条与数值积分求一阶导时，误差是不同的，而在试验中给定 的数据显示：在一定的误差范围下，用数值积分求导比用三次样条求得的导数值 更加接近真值。实验五LU和列主元高斯消去法解方程组一、实验目的用LU分解和列主元高斯消去法求解上述两个方程组，输出Ax=b中矩阵A 及向量b，A = LU分解的L，U及det（A），解向量x（1）用LU分解和列主元高斯消去法求解上述方程组（中的数据取小数点 后三位）输出Ax = b中矩阵A及向量b，A = LU分解的L, U , det（A）及解向量x.（2）将方程组中系数3.01改为3.00，0.987改为0.990 .用列主元高斯消去 法求解，输出列主元行交换次序、解向量x及det（A）,并与（1）中结果比较.二、实验原理LU分解：A = LU，L是单位下三角阵，U是上三角阵。求解Ax = b的问题就等价于求解两个三角阵方程组Ly = 8,求y , Ux - y U,求x。列主兀高斯消去法：1. det 12. 对于 k = 1,2 ,n -1(1) .按列选主元 a ,k = max. a(2) .如果a，=。，则停止计算det(A)= 0lk(3) .若果lk = k，则转换行：a f a (j = k, k +1,n) b f b det <detkllk,jklk(4) .消元计算对于l = k +1,n(a) a m =气匕kk(b) 对于 j = k +1,n(c) b b 一 m * b(5) . det a * det3.如果a = 0，则停止计算det (A)= 04.回代求解 bn "ann(2) 对于 l = n -1,-,2,1b - U alljb l j=l+1l5. det a * detnn三、实验结果所用数据：3.011.270.9871.99 -尤1了-1.23尤_129.34X316.034.16-4.81计算结果为:LU分解法对应的L为：LU分解法对应的U为：1.0000000.42191.000000.3279-4.20061.00003.01006.03001.990001.6158-2.069600-0.0063LU分解法对应的向量解x1.0e+003 *1.5926-0.6319-0.4936行列式A的值det-0.030547列主元高斯消去法消去后的矩阵A:把第2行和第3行互换列主元高斯消去法对应的向量解x:3.01006.03001.99000-6.78738.687500.0000-0.00151.0e+003 *1.5921-0.6318-0.4936行列式A的值det-0.030547将系数矩阵中3.01改为3.00，0.987 改为 0.990列主元高斯消去法消去后的矩阵A:把第2行和第3行互换列主元高斯消去法对应的向量解x3.00006.03001.99000-6.79998.683300-0.0200118.9960 -47.0441 -36.8403行列式A的值det-0.407014四、实验分析通过实验我们可以知道，用LU分解和列主元高斯消去法解方程组时，结果是比较精确地。而系数的微小改变会导致最后的结果有很大的变化。实验六迭代法求解方程一、实验目的给定的迭代误差要求下，通过对一个方程组使用雅克比迭代、高斯-塞德尔 迭代法以及SOR迭代法，观察得到的序列是否收敛。若收敛，记录收敛次数，并 研究迭代法的收敛速度和什么因素有关。二、实验原理J X (0)X (k +1) = Bx (k) + f雅克比迭代：B = D-1(L + D), f = D-ibX1xi = 1,2,n(k = 0,1,表示迭代次数)高斯-塞德尔迭代：B = (D - L)-1U , f = (D L)-1 bX（0） = （0）,，X（。）力f b-正 aX k+!）"'尸1i = 1,2,n(k = 0,1,表示迭代次数)SOR 迭代法：B = (D-L)-1(1-)D + U)，f =rn(D L)-1 bx(0)= X (o),x 1nj+ij=1(k+!)x (k+1)= x (k)iiAx.j=i+1i = 1,2,，n(k = 0,1,)g为松弛因子三、实验结果所用数据：A=diag(3*1:1:19,0)+diag(-1:1:18/3,-1)+diag(-ones(1,18)/3,1);b=ones(19,1);wuch=1e-6;计算结果为：Jocabi迭代法求解的解向量：Jocabi迭代次数为60.3548 0.1935 0.1290 0.0962 0.0766 0.0637 0.0544 0.04750.0422 0.0379 0.0345 0.0316 0.0291 0.0270 0.0252 0.02360.0222 0.0210 0.0198Gauss_Seidel 迭代法 求解的解向量：Gauss_Seidel 迭代次数为50.3548 0.1935 0.1290 0.0962 0.0766 0.0637 0.0544 0.04750.0422 0.0379 0.0345 0.0316 0.0291 0.0270 0.0252 0.02360.0222 0.0210 0.0198SOR迭代法当松弛因子为1.00时 SOR迭代次数为4 当松弛因子为1.05时 SOR迭代次数为4 当松弛因子为1.10时 SOR迭代次数为6 当松弛因子为1.15时 SOR迭代次数为7 当松弛因子为1.20时 SOR迭代次数为8 当松弛因子为1.25时 SOR迭代次数为10 当松弛因子为1.30时 SOR迭代次数为12 当松弛因子为1.35时 SOR迭代次数为13 当松弛因子为1.40时 SOR迭代次数为16 当松弛因子为1.45时 SOR迭代次数为19 当松弛因子为1.50时 SOR迭代次数为23 当松弛因子为1.55时 SOR迭代次数为28当松弛因子为1.60时 当松弛因子为1.65时 当松弛因子为1.70时 当松弛因子为1.75时 当松弛因子为1.80时 当松弛因子为1.85时 当松弛因子为1.90时 当松弛因子为1.95时 当松弛因子为2.00时SOR迭代次数为35 SOR迭代次数为45 SOR迭代次数为62 SOR迭代次数为89 SOR迭代次数为126 SOR迭代次数为167 SOR迭代次数为256SOR迭代法不收敛 SOR迭代法不收敛四、实验分析通过试验，给定的迭代误差要求下，使用雅克比迭代、高昕塞德尔迭代法、 SOR迭代法，其收敛次数是不同的；而如果所构建的迭代阵地谱半径越小，则其 收敛速度越快。实验七牛顿法解非线性方程组、实验目的用牛顿法求解非线性方程组的根，输出迭代次数及解向量尤,y,z的近似值、实验原理牛顿迭代法：If £ ,七）=0FG）= 0 （雅可比矩阵）%nn3, k=0,1 F'(x (k)三、实验结果尤2 - y2 - 4 = 0,< x2 + 5y2 + 4Z2 - 20 = 0,x, y,却。）=-1,1,1 x2 + y2 + 4z2 一 10 = 0.计算结果:Newton法第1次迭代的解向量x y z为：x=-1 y=1 z=1Newton法第2次迭代的解向量x y z为x=-3.75 y=1.75 z=0.625Newton法第3次迭代的解向量x y z为x=-2.7417 y=1.5893 z=0.5125Newton法第4次迭代的解向量x y z为x=-2.5562 y=1.5812 z=0.50015Newton法第5次迭代的解向量x y z为x=-2.5495 y=1.5811 z=0.5Newton迭代次为6由Newton法解得方程组近似解x y z为x=-2.5495 y=1.5811 z=0.5四、实验分析利用牛顿法解非线性方程组时，其迭代次数与其初始向量有关，其结果的 精确程度也是一样。实验八QR法求特征值一、实验目的(1) 根据QR算法原理编制求矩阵A全部特征值的程序并输出计算结果(要求误差 10 -6).(2) 直接用现有数学软件求矩阵A的全部特征值，并与(1)的结果比较二、实验原理带原点位移的QR算法设 A = A & Rnxn对A - 5 I进行QR分解A - s I = QR形成矩阵 A = RQ + s I = Qt(A - sI)2 + s I = QtAQ21111111111求的Ak后将Ak-skl进行QR分解A - sj = QR,k = 3,4形成矩阵 Ak+1 = RQ + sj = Q：A QkA = (QQ -Q )tA (QQ -Q ) = QtAQk+11 2 k 1 1 2 k k 1 kAk的QR分解式为Ak = Q Rk。三、实验结果所用数据：52212-41121311112计算结果为:此矩阵用数学软件的特征值为：用QR法求A的特征值为：77.64709.7001 - 9.5372i9.7001 + 9.5372i0.1956 - 0.0000i9.7001 - 9.5372i4.7571 + 0.0000i0.19569.7001 + 9.5372i4.757177.6470四、实验分析带原点位移的QR算法收敛速度快，算法比较稳定，而且其求得的矩阵特征 值也是相当精确的。实验九 改进欧拉及四阶R-K方法求初值问题一、实验目的分别用改进欧拉法及经典四阶R-K.方法求初值问题的数值解，输出两种方法在七的数值解七，与精确解比较并分析结果，输出精确解和两种近似解的图形二、实验原理改进欧拉公式:四阶R-K法公式:y = y + -(K + 2 K + 2 K + K ) n+1n 6 1234Ki=f(xn，y)L _/ h h )< K f x + ，y + K2 J " n 2 n 2 1)K = f x + h, y + 如K '3 " n 2 七 2 2)K = f (x + h, y + hK )三、实验结果所用数据:y '=1-晟精确解：y（x）=y (0) = 0.x( x 2 + 3) 3(1+ x 2)x = 0.05海=0,1, ,20）计算结果为：h = 0.05x,对应的精确解：x,对应改进欧拉法求初值问 题数值解为：x.对应经典四阶R-K方 法求初值问题数值解为f(0)=0y(0)=0y(0)=0f(0.05)=0.049917y(0.05)=0.049875y(0.05)=0.041604f(0.1)=0.09934y(0.1)=0.099258y(0.1)=0.091089f(0.15)=0.1478y(0.15)=0.14768y(0.15)=0.13965f(0.2)=0.19487y(0.2)=0.19472y(0.2)=0.18686f(0.25)=0.2402y(0.25)=0.24001y(0.25)=0.23235f(0.3)=0.28349y(0.3)=0.28327y(0.3)=0.27584f(0.35)=0.32454y(0.35)=0.3243y(0.35)=0.31711f(0.4)=0.36322y(0.4)=0.36296y(0.4)=0.35603f(0.45)=0.39948y(0.45)=0.39921y(0.45)=0.39255f(0.5)=0.43333y(0.5)=0.43306y(0.5)=0.42667f(0.55)=0.46484y(0.55)=0.46456y(0.55)=0.45845f(0.6)=0.49412y(0.6)=0.49384y(0.6)=0.48799f(0.65)=0.52129y(0.65)=0.52102y(0.65)=0.51544f(0.7)=0.54653y(0.7)=0.54626y(0.7)=0.54094f(0.75)=0.57y(0.75)=0.56973y(0.75)=0.56467f(0.8)=0.59187y(0.8)=0.59161y(0.8)=0.58679f(0.85)=0.61231y(0.85)=0.61206y(0.85)=0.60747f(0.9)=0.63149y(0.9)=0.63125y(0.9)=0.62689f(0.95)=0.64956y(0.95)=0.64933y(0.95)=0.64518f(1)=0.66667y(1)=0.66645y(1)=0.6625画出的图形如下所示:四、实验分析试验中我们选择在步长相同的情况下，四阶R-K法对于常微分方程所求的 解的精确程度不如用改进的欧拉方法；而结果的精确程度是与步长的选取有关 的。代码：function shuzhishiyanclcfprintf（' 数 值 分 析 实 验n'）;fprintf（'nn'）;fprintf（'实验一n'）;fprintf（'实验二n'）;fprintf（'实验三n'）;fprintf（'实验四n'）;fprintf（'实验五n'）;fprintf（'实验六n'）;fprintf（'实验七n'）;fprintf（'实验八n'）;fprintf（'实验九n'）;fprintf（'按任意键继续n'）;pause;sanciyangtiao;fprintf（'按任意键继续n'）;pause;closeErCheng;fprintf（'按任意键继续n'）;pause;closelongbeige;fprintf（'按任意键继续n'）;pause;sanciqiudao;fprintf（'按任意键继续n'）;pause;LU_Gauss;fprintf（'按任意键继续n'）;pausediedaifa;fprintf（'按任意键继续n'）;pause;Newton;fprintf（'按任意键继续n'）;pause;QR;fprintf（'按任意键继续n'）;pause;chuzhishuzhijie;fprintf(，数值分析实验结束n');endfunction sanciyangtiaoclearclcformat shortfprintf(，实验一：三次样条插值的三弯矩法，);x=-2.223,-1.987,-1.8465,-1.292,-1.2266,-1.1056,-0.8662,-0.6594,-0.2671,-0.0452,0.5385,1.2564,1.4398,1.5415,1.7646,1.9678,2.236;y=0.83995,1.1696,1.3141,1.6992,1.7312,1.7847,1.8708,1.9262,1.9881,1.9997,1.9511,1.7169,1.618,1.5543,1.3871,1.191,0.81662;d2s1=-4.5000;d2sn= -4.8967; %第二种边界条件t=-2.223,-1.9443,-1.6656,-1.3869,-1.1083,-0.82956,-0.55088,-0.27219,0.0065,0.28519,0.56387,0.84256,1.1212,1.3999,1.6786,2.236;n=length(x);h=zeros(n,1);for i=1:n-1h(i)=x(i+1)-x(i);enda,v,u,d=sanwanju(x,y,d2s1,d2sn);M=zhuigan(u,v,d);k=length(t);flag=0;%求插值点所对应的函数值for i=1:kfor j=1:n-1if (t(i)>=x