MATLAB大作业 给定一个时间序列,使用三次样条插值方法进行均匀内插.doc

MATLAB作业给定一个时间序列，使用三次样条插值方法进行均匀内插（题目的相关说明：按题目要求编写一个MATLAB程序函数，并把自己编制程序所得的结果与MATLAB库函数分析结果进行对比。）理论基础：时间序列的概念：时间序列是一种定量预测方法，又称简单外延法，时间序列分析是根据系统观测得到的时间序列数据，通过曲线拟合和参数估计来建立数学模型的理论与方法，时间序列分析可分为以下三种情况（1）把一个时间序列的值变动为N 个组成部分，通常可以分为四种 a、倾向变动，又称长期趋势变动 b、循环变动，又称周期变动 c、季节变动，即每年有规则的反复进行变动 d、不规则变动，即随机变动。然后把这四个综合到一起得出预测的结果。虽然分成这四部分，但这四部分之间的相互关系是怎么样的呢，目前一般采用相乘的关系，其实各个部分都是在其他部分作用的基础上进行作用的，所以采用相乘是有一定依据的，此种方法适合于短期预测和库存预测（2）把预测对象、预测目标和对预测的影响因素都看成为具有时序的，为时间的函数，而时间序列法就是研究预测对象自身变化过程及发展趋势，如果未来预测是线性的，其数学模型为YT+L=aT+bTL,YT+L为未来预测值，aT为截距，bT为斜率，L为由T到需要预测的单位时间数（如5年、10年等）（3）根据预测对象与影响因素之间的关系及影响程度来推算未来，与目标的相关因素很多，只能选择那些因果关系较强的为预测影响的因素，此时间序列法用于短期预测比较有效，若要用于长期预测，还需要结合其他方法才行。三次样条插值的实际应用：在制造船体和汽车外形等工艺中传统的设计方法是，首先由设计人员按外形要求，给出外形曲线的一组离散点值，施工人员准备好有弹性的样条（一般用竹条或有弹性的钢条）和压铁，将压铁放在点的位置上，调整竹条的形状，使其自然光滑，这时竹条表示一条插值曲线，我们称为样条函数。从数学上看，这一条近似于分段的三次多项式，在节点处具有一阶和二阶连续微商。样条函数的主要优点是它的光滑程度较高，保证了插值函数二阶导数的连续性，对于三阶导数的间断，人类的眼睛已难以辨认了。样条函数是一种隐式格式，最后需要解一个方程组，它的工作量大于多项式拉格朗日型式或牛顿型式等显式插值方法。三次样条插值函数的定义 设在区间上，被插函数y=f(x)在个节点 和这些点上的函数值分别为。若满足；在每个小区间上至多是一个三次多项式；在上有连续的二阶导数，则称为关于剖分的三次样条插值函数，称为样条节点。要在每个子区间上构造三次多项式，共需要个条件，由插值条件，提供了个条件；用每个内点的关系建立条件又得到个条件；再附加两个边界条件，即可惟一确定样条函数了。用待定系数法确定了构造样条函数的存在性和惟一性。在具体构造样条函数时一般都不使用计算量大的待定系数法。边界条件有三类：（1)m边界条件：给定曲线在两端点处的切线斜率，s(X0)=m0,s(Xn)=Mn.即两个边界节点的一阶导数值为给定值：m0,Mn. (2) M边界条件：令被插函数在两端点的二阶导数为零，s(X0)=m0, s(Xn)=Mn.即两个边界节点的二阶导数值为给定值：m0,Mn.(3) 周期性边界条件：令被插值函数为周期函数，周期为b-a=xn-x0,故有y0=yn.相应的也要求样条函数S(x)也是周期函数，在端点上有s(X0)=s(Xn); s(X0)=s(Xn).算法原理：例如：x=0 1 2 3 456 7 8 9 10; y=0 2 0 -4 0 40 -2 0 3 1;待求的三次多项式函数s(x)在0 10上有连续的一阶，二阶导数.我们通过简单的讨论来认识问题。在第一区间0 1、第二区间1 2上考虑两个三次多项式 s(x)=s1*x3+s2*x2+s3*x+s4 r(x)=r1*x3+r2*x2+r3*x+r4示意图:可以得到 s(0)=s1*03+s2*02+s3*0+s4=0 (1) s(1)=s1*13+s2*12+s3*1+s4=2 (2) r(1)=r1*23+r2*22+r3*2+r4=2 (3) r(2)=r1*23+r2*22+r3*2+r4=0 (4)一阶导函数 s(x)=3*s1*x2+2*s2*x+s3 r(x)=3*r1*x2+2*r2*x+r3由一阶导数的连续性且在1点处相等，有 3*s1*12+2*s2*1+s3=3*r1*12+2*r2*1+r3 (5)二阶导函数 s(x)=6*s1*x+2*s2 r(x)=6*r1*x+2*r2 由二阶导数的连续性且在1点处相等，有 6*s1*1+2*s2=6*r1*1+2*r2 (6) 由m边界条件 s(0)=1.6,r(2)=0.3有 3*s1*02+2*s2*0+s3=1.6 (7) 3*r1*22+2*r2*2+r3=0.3 (8) M边界条件 s(0)=-1,r(2)=1有 6*s1*0+2*s2=-1 (7) 6*r1*2+2*r2=1 (8) 由周期性边界条件 s(0)=r(2) s(0)=r(2)有 3*s1*02+2*s2*0+s3=-1 (7) 3*r1*22+2*r2*2+r3=1 (8)这样，对于两个多项式的8个未知量，我们给出了8个方程。三次样条曲线的难点在于，我们不能分段去求解方程，完成绘图。（本题中为了简化运算以及和MATLAB程序结合起来，使用下面两式3*s1*02+2*s2*0+s3=03*r1*22+2*r2*2+r3=0（把它当作自然条件来实现）2. 算法步骤： 流程图： 开始 输入一个时间序列输入一个初始i值不满足N满足Yi11函数语句函数曲线曲线 i=i+1结束步骤： 1)输入一个时间序列 2）输入一个初始值 3）进行循环 4)通过程序得出函数曲线 5）依次得出的曲线连接成一段完整的曲线 6)超出时间序列即刻终止，程序运行结束3. MATLAB程序 : x=0,100,200,300,400,500,600,700,800,900,1000;y=0,2,0,-4,0,4,0,-2,0,3,1;tici=1;while i10; A=x(1,i).3,x(1,i).2,x(1,i),1,0,0,0,0; x(1,i+1).3,x(1,i+1).2,x(1,i+1),1,0,0,0,0; 0,0,0,0, x(1,i+1).3,x(1,i+1).2,x(1,i+1),1; 0,0,0,0, x(1,i+2).3,x(1,i+2).2,x(1,i+2),1; 3*x(1,i+1).2,2*x(1,i+1),1,0,-3*x(1,i+1).2,-2*x(1,i+1),-1,0; 6*x(1,i+1),2,0,0,-6*x(1,i+1),-2,0,0; 3*x(1,i).2,2*x(1,i),1,0,0,0,0,0; 0,0,0,0,3*x(1,i+2).2,2*x(1,i+2),1,0;B=y(1,i),y(1,i+1),y(1,i+1),y(1,i+2),0,0,0,0;c=ABs1=c(1,1),s2=c(2,1),s3=c(3,1),s4=c(4,1); x1=100*(i-1):100*i; figure(1);y1=s1*x1.3+s2*x1.2+s3*x1+s4; plot(x1,y1);hold on;i=i+1;endtoc运行结果：c = -0.0000 0.0006 0 0 0.0000 -0.0018 0.2400 -8.0000s1 = -4.0000e-006s2 = 6.0000e-004s3 = 0c = -0.0000 0.0000 0.0150 1.0000 0.0000 -0.0024 0.4950 -31.0000s1 = -5.0000e-007s2 = 1.7319e-018s3 = 0.0150c = 0.0000 -0.0060 1.4400 -112.0000 -0.0000 0.0084 -2.8800 320.0000s1 = 8.0000e-006s2 = -0.0060s3 = 1.4400c = -0.0000 0.0024 -0.9000 104.0000 -0.0000 0.0024 -0.9000 104.0000s1 = -2.0000e-006s2 = 0.002s3 = -0.9000c = 1.0e+003 * -0.0000 0.0000 -0.0048 0.7040 0.0000 -0.0000 0.0072 -1.2960s1 = -8.0000e-006s2 = 0.0108s3 = -4.8000c = 0.0000 -0.0060 3.3750 -621.0000 -0.0000 0.0012 -0.9450 243.0000s1 = 3.5000e-006s2 = -0.0060s3 = 3.3750c = 1.0e+003 * 0.0000 -0.0000 0.0050 -1.0800 -0.0000 0.0000 -0.0067 1.6640s1 = 4.0000e-006s2 = -0.0078s3 = 5.0400c = 1.0e+003 * -0.0000 0.0000 -0.0007 0.1940 -0.0000 0.0000 -0.0045 1.2180s1 = -2.5000e-007s2 = 7.5000e-004s3 = -0.6825c = 1.0e+003 * -0.0000 0.0000 -0.0114 3.2160 0.0000 -0.0000 0.0129 -4.0740s1 = -5.2500e-006s2 = 0.0134s3 = -11.4000MATLAB库函数c程序： x=0,200,400,600,800,1000,1200,1400,1600,1800,2000;x1=0:1:2000;y=0,-20,0,16,0,-12,0,8,0,-4,0;ticy1=interp1(x,y,x1,spline);plot(x,y,ok,x1,y1,-r);toc库函数运行结果： 4.结果分析：如上两幅图进行对比看出编制的程序函数基本能达到题目要求，但与MATLAB库函数得到