三次样条插值的MATLAB实现

合集下载

matlab实现三次样条插值法

matlab实现三次样条插值法

题目背景:对y=1/(1+x^2)在[-1,1]区间以Xn=-1+0.1*(n-1),n=1 (21)为插值点做三次样条插值求解思路简析:以插值为四段三次函数为例进行说明(题干为插值20段三次函数),可看出方程组为q*x=d,其中q为方程组系数矩阵,x为所求三次函数的系数矩阵,其中方程组系数矩阵和d均呈规律性变化(边界点除外,首位两个点特殊堪虑)function qiujieyangtiao %%定义求解函数q=zeros(80); %%方程组的系数矩阵,赋初值为0n=-1:0.1:1; %%插值点的横坐标nd=zeros(80,1); %%插值点q*x=d中的dy=zeros(21,1); %%插值点的纵坐标向量a=1;for i=-1:0.1:1y(a)=1/(1+i^2);a=a+1; %%给插值点的纵坐标y通过原函数赋值endq(1,3)=2;q(1,4)=6*n(1);q(2,1)=1;q(2,2)=n(1);q(2,3)=n(1)^2;q(2,4)=n(1)^3;d(2)=y(1); %%给左端边界点的两个方程组系数赋值j=2;for i=3:4:75q(i,i-1)=1;q(i,i)=2*n(j);q(i,i+1)=3*n(j)^2;q(i,i+3)=-1;q(i,i+4)=-2*n(j);q(i,i+5)=-3*n(j)^2;d(i)=0;q(i+1,i)=2;q(i+1,i+1)=6*n(j);q(i+1,i+4)=-2;q(i+1,i+5)=-6*n(j);d(i+1)=0;q(i+2,i-2)=1;q(i+2,i-1)=n(j);q(i+2,i)=n(j)^2;q(i+2,i+1)=n(j)^3;d(i+2)=y(j);q(i+3,i+2)=1;q(i+3,i+3)=n(j);q(i+3,i+4)=n(j)^2;q(i+3,i+5)=n(j)^3;d(i+3)=y(j);j=j+1;end %%给系数矩阵赋值q(79,79)=2;q(79,80)=6*n(21);d(79)=0;q(80,77)=1;q(80,78)=n(21);q(80,79)=n(21)^2;q(80,80)=n(21)^3;d(80)=y(21); %%给右端边界点的两个方程组系数赋值result=q\d; %%求解系数矩阵function A=fun(x)if x>=-1&&x<-0.9A=result(1)+result(2)*x+result(3)*x*x+result(4)*x*x*x;elseif x>=-0.9&x<-0.8A=result(5)+result(6)*x+result(7)*x*x+result(8)*x*x*x;elseif x>=-0.8&x<-0.7A=result(9)+result(10)*x+result(11)*x*x+result(12)*x*x*x; elseif x>=-0.7&x<-0.6A=result(13)+result(14)*x+result(15)*x*x+result(16)*x*x*x; elseif x>=-0.6&x<-0.5A=result(17)+result(18)*x+result(19)*x*x+result(20)*x*x*x; elseif x>=-0.5&x<-0.4A=result(21)+result(22)*x+result(23)*x*x+result(24)*x*x*x; elseif x>=-0.4&x<-0.3A=result(25)+result(26)*x+result(27)*x*x+result(28)*x*x*x; elseif x>=-0.3&x<-0.2A=result(29)+result(30)*x+result(31)*x*x+result(32)*x*x*x; elseif x>=-0.2&x<-0.1A=result(33)+result(34)*x+result(35)*x*x+result(36)*x*x*x; elseif x>=-0.1&x<0A=result(37)+result(38)*x+result(39)*x*x+result(40)*x*x*x; elseif x>=0&x<0.1A=result(41)+result(42)*x+result(43)*x*x+result(44)*x*x*x; elseif x>=0.1&x<0.2A=result(45)+result(46)*x+result(47)*x*x+result(48)*x*x*x; elseif x>=0.2&x<0.3A=result(49)+result(50)*x+result(51)*x*x+result(52)*x*x*x; elseif x>=0.3&x<0.4A=result(53)+result(54)*x+result(55)*x*x+result(56)*x*x*x; elseif x>=0.4&x<0.5A=result(57)+result(58)*x+result(59)*x*x+result(60)*x*x*x; elseif x>=0.5&x<0.6A=result(61)+result(62)*x+result(63)*x*x+result(64)*x*x*x; elseif x>=0.6&x<0.7A=result(65)+result(66)*x+result(67)*x*x+result(68)*x*x*x; elseif x>=0.7&x<0.8A=result(69)+result(70)*x+result(71)*x*x+result(72)*x*x*x; elseif x>=0.8&x<0.9A=result(73)+result(74)*x+result(75)*x*x+result(76)*x*x*x; elseA=result(77)+result(78)*x+result(79)*x*x+result(80)*x*x*x; endend %%插值函数用子函数表达,方便调用x=linspace(-1,1);for i=1:length(x)A(i)=fun(x(i));endY=1./(1+x.^2);plot(x,Y,'--',x,A,':')legend('primitive','fitting') %%将原函数与该插值函数画在同一图上进行比较grid ontitle('三次样条插值')for m=1:20fprintf("S%d=%.3f+%.3f*x+%.3f*x.^2+%.3f*x.^3\n",m,result(4*m-3,1),result(4*m-2,1),result(4*m-1,1),result(4*m,1)) %%输出结果endend输出结果:S1=2.049+3.619*x+3.104*x.^2+1.035*x.^3S2=1.010+0.156*x+-0.743*x.^2+-0.390*x.^3S3=1.137+0.632*x+-0.149*x.^2+-0.143*x.^3S4=1.054+0.273*x+-0.660*x.^2+-0.386*x.^3S5=1.023+0.120*x+-0.916*x.^2+-0.528*x.^3S6=1.003+-0.002*x+-1.160*x.^2+-0.691*x.^3S7=0.997+-0.044*x+-1.265*x.^2+-0.779*x.^3S8=0.998+-0.034*x+-1.233*x.^2+-0.743*x.^3S9=1.000+-0.010*x+-1.113*x.^2+-0.543*x.^3S10=1.000+-0.000*x+-1.010*x.^2+-0.200*x.^3S11=1.000+-0.000*x+-1.010*x.^2+0.200*x.^3S12=1.000+0.010*x+-1.113*x.^2+0.543*x.^3S13=0.998+0.034*x+-1.233*x.^2+0.743*x.^3S14=0.997+0.044*x+-1.265*x.^2+0.779*x.^3S15=1.003+0.002*x+-1.160*x.^2+0.691*x.^3S16=1.023+-0.120*x+-0.916*x.^2+0.528*x.^3S17=1.054+-0.273*x+-0.660*x.^2+0.386*x.^3S18=1.137+-0.632*x+-0.149*x.^2+0.143*x.^3S19=1.010+-0.156*x+-0.743*x.^2+0.390*x.^3S20=2.049+-3.619*x+3.104*x.^2+-1.035*x.^3对比图。

MATLAB数值实验一(数据的插值运算及其应用完整版)

MATLAB数值实验一(数据的插值运算及其应用完整版)

佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 插值法与数据拟合 专业班级 机械工程 姓 名 余红杰 学 号 10 指导教师 陈剑 成 绩 日 期 月 日一、实验目的1、学会Lagrange 插值、牛顿插值和三次样条插值等基本插值方法;2、讨论插值的Runge 现象3、学会Matlab 提供的插值函数的使用方法,会用这些函数解决实际问题。

二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、三次样条插值 三、实验步骤1、用MATLAB 编写独立的拉格朗日插值多项式函数2、用MATLAB 编写独立的牛顿插值多项式函数3、用MATLAB 编写独立的三次样条函数(边界条件为第一、二种情形)4、已知函数在下列各点的值为:根据步骤1,2,3编好的程序,试分别用4次拉格朗日多项式4()L x 、牛顿插值多项式4()P x 以及三次样条函数()S x (自然边界条件)对数据进行插值,并用图给出 {(,),0.20.08,0,1,2,,10i i i x y x i i =+=},4()L x 、4()P x 和()S x 。

5、在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数21(),(11)125f x x x=-≤≤+作多项式插值,对不同n 值,分别画出插值函数及()f x 的图形。

6、下列数据点的插值可以得到平方根函数的近似,在区间[0,64]上作图。

(1)用这9个点作8次多项式插值8()L x 。

(2)用三次样条(第一边界条件)程序求()S x 。

7、对于给函数21()125f x x =+在区间[-1,1]上取10.2(0,1,,10)i x i i =-+=,试求3次曲线拟合,试画出拟合曲线并打印出方程,与第5题的结果比较。

四、实验过程与结果:1、Lagrange 插值多项式源代码:function ya=lag(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点,自变量 %ya 所求点插值估计量 ya=0; mu=1; %初始化%循环方式求L 系数,并求和: for i = 1:length(y) for j = 1:length(x) if i ~= jmu = mu * (xa - x(j) ) / ( x(i) - x(j) ); else continue end endya = ya + y(i) * mu ; mu = 1; end2、Newton 源代码:function ya = newton(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点,自变量 %ya 所求点插值估计量 %建立系数零矩阵D 及初始化:D = zeros(length(x)-1);ya = y(1);xi = 1;%求出矩阵D,该矩阵第一行为牛顿插值多项式系数:for i=1:(length(x)-1)D(i,1) = (y(i+1) -y(i))/(x(i+1) -x(i));endfor j=2:(length(x)-1)for i=1:(length(x)-j)D(i,j) = (D(i+1,j-1) - D(i,j-1)) / (x(i+j) - x(i)); endend%xi为单个多项式(x-x(1))(x-x(2))...的值for i=1:(length(x)-1)for j=1:ixi = xi*(xa - x(j));endya = ya + D(1,i)*xi;xi = 1;end3、三次样条插值多项式(1)(第一边界条件)源代码: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]4、——————————————PS:另建了一个f方程文件,后面有一题也有用到。

Matlab实验报告六(三次样条与分段线性插值)范文

Matlab实验报告六(三次样条与分段线性插值)范文
1.分析问题
本题是给出粗略等分点让你插入更多点用双线性插值法来作出更清晰的山区地貌图。
2.问题求解
x=0:400:2800;
y=0:400:2400;
z=[1430 1450 1470 1320 1280 1200 1080 940;
1450 1480 1500 1550 1510 1430 1300 1200;
2.分段线性插值与计算量与n无关;n越大,误差越小.
3.三次样条插值比分段线性插值更光滑。
4.‘linear’:分段线性插值;‘spline’:三次样条值。
【实验环境】
MatlabR2010b
二、实验内容
问题1对函数 ,x[-5,5],分别用分段线性插值和三次样条插值作插值(其中插值节点不少于20),并分别作出每种插值方法的误差曲线.
本次实验因为是我们课本没有的内容,心理上给了我很大的压力,幸好我们还能根据老师的课件以及例题去掌握这次实验所需要的各种插值法,但结果还好,两道题都做出来了。
plot(x,y,'*',x1,yl,'r',x1,y2,'b')
y0=1./(1+x1.^2);
y3=yl-y0;
y4=ys-y0;
holdon
plot(x1,y3,'y',x1,y4,'g')
3.结果
4误。
问题2山区地貌图在某山区(平面区域(0,2800)(0,2400)内,单位:米)测得一些地点的高程(单位:米)如表1,试作出该山区的地貌图.
1.分析问题
本题先取出少量的插值节点并作出图形,再用分段线性插值法和三次样条插值法做出更精确的图形,最后在作出误差曲线。

matlab三次样条插值例题解析

matlab三次样条插值例题解析

文章标题:深度解析Matlab三次样条插值1. 前言在数学和工程领域中,插值是一种常见的数值分析技术,它可以用来估计不连续数据点之间的值。

而三次样条插值作为一种常用的插值方法,在Matlab中有着广泛的应用。

本文将从简单到复杂,由浅入深地解析Matlab中的三次样条插值方法,以便读者更深入地理解这一技术。

2. 三次样条插值概述三次样条插值是一种利用分段三次多项式对数据点进行插值的方法。

在Matlab中,可以使用spline函数来进行三次样条插值。

该函数需要输入数据点的x和y坐标,然后可以根据需要进行插值操作。

3. 三次样条插值的基本原理在进行三次样条插值时,首先需要对数据点进行分段处理,然后在每个分段上构造出一个三次多项式函数。

这些多项式函数需要满足一定的插值条件,如在数据点处函数值相等、一阶导数相等等。

通过这些条件,可以得到一个关于数据点的插值函数。

4. Matlab中的三次样条插值实现在Matlab中,可以使用spline函数来进行三次样条插值。

通过传入数据点的x和y坐标,可以得到一个关于x的插值函数。

spline函数也支持在已知插值函数上进行插值点的求值,这为用户提供了极大的灵活性。

5. 三次样条插值的适用范围和局限性虽然三次样条插值在许多情况下都能够得到较好的插值效果,但也存在一些局限性。

在数据点分布不均匀或有较大噪音的情况下,三次样条插值可能会出现较大的误差。

在实际应用中,需要根据具体情况选择合适的插值方法。

6. 个人观点和总结通过对Matlab中三次样条插值的深度解析,我深刻地理解了这一插值方法的原理和实现方式。

在实际工程应用中,我会根据数据点的情况选择合适的插值方法,以确保得到准确且可靠的结果。

我也意识到插值方法的局限性,这为我在实际工作中的决策提供了重要的参考。

通过以上深度解析,相信读者已经对Matlab中的三次样条插值有了更加全面、深刻和灵活的理解。

在实际应用中,希望读者能够根据具体情况选择合适的插值方法,以提高工作效率和准确性。

三次样条插值函数MATLAB编程实现

三次样条插值函数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 x x x x x x x x x x x x x x x x x x x x x x x x xx ααββ+++++++++++++++''=---+''=-+--+''=---+''=---则有()()()()()()()()()()()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 jj j j j j j jj j j j jj j j j j jjj jj jj jj x x y x x y x x x x m m x x xx xx xx x x y x x y x x x x m m x x xx xx xx ++++++++++--------------+++--------=+++----(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 ----------⨯⎛⎫⎪ ⎪ ⎪ ⎪⎪⎪ ⎪-- ⎪-+- ⎪--⎝⎭为了使样条插值问题有惟一解,我们在原有方程基础上增加两个边界条件。

MATLAB三次样条插值之三弯矩法

MATLAB三次样条插值之三弯矩法

MATLAB三次样条插值之三弯矩法首先说这个程序并不完善,为了实现通用(1,2,…,n)格式解题,以及为调用追赶法程序,没有针对节点数在三个以下的情况进行分类讨论。

希望能有朋友给出更好的方法。

首先,通过函数sanwanj得到方程的系数矩阵,即追赶法方程的四个向量参数,接下来调用追赶法(在intersanwj函数中),得到三次样条分段函数系数因子,然后进行多项式合并得到分段函数的解析式,程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算改点的值。

附:追赶法程序chase%%%%%%%%%%%%%%function [newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b)ﻫ%三弯矩样条插值ﻫ%将插值点分两次输入,x0y0单独输入ﻫ% 边值条件a的二阶导数 y1a 和b的二阶导数y1b,这里建议将y1a和y1b换成y2a和y2b,以便于和三转角代码相区别ﻫn=length(x);m=length(y);if m~=nﻫerror('x or y 输入有误,再来');endﻫv=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1);ﻫw=2*ones(n+1);ﻫh0=x(1)-x0;ﻫh=zeros(n-1,1);for k=1:n-1ﻫh(k)=x(k+1)-x(k);ﻫendv(1)=h0/(h0+h(1));u(1)=1-v(1);d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1));ﻫ%for k=2:n-1ﻫv(k)=h(k-1)/(h(k-1)+h(k));ﻫu(k)=1-v(k);ﻫd(k)=6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k));endnewv=[v;1];ﻫnewu=[1;u];d0=6*((y(1)-y0)/h0-y1a)/h0;d(n)=6*(y1b-(y(n)-y(n-1))/h(n-1))/h(n-1);newd=[d0;d];%%%%%%%%%%%%function intersanwj(x,y,x0,y0,y1a,y1b)%三弯矩样条插值ﻫ%第一部分ﻫn=length(x);m=length(y);if m~=nﻫerror('xory 输入有误,再来');endﻫ%重新定义hﻫh=zeros(n,1);h(1)=x(1)-x0;for k=2:nh(k)=x(k)-x(k-1);ﻫend%sptep1调用三弯矩函数ﻫ[a,b,c,d]=sanwj(x,y,x0,y0,y1a,y1b);% 三对角方程ﻫM=chase(a,b,c,d);% 求插值函数ﻫfprintf('三次样条(三弯矩)插值的函数表达式\n');syms X ;ﻫfprintf('S0--1:\n');S(1)=collect(((1/6)*M(2)*(X-x0).^3-(1/6)*M(1)*(X-x(1)).^3+(y(1)-(M(2)*h(1).^2)/6)*(X-x0)-(y0-(M(1)*h(1).^2)/6)*(X-x(1)))/h(1));ﻫfor k=2:nfprintf('S%d--%d:\n',k-1,k);S(k)=collect(((1/6)*M(k+1)*(X-x(k-1)).^3-(1/6)*M(k)*(X-x(k)).^3+(y(k)-(M(k+1)*h(k).^2)/6)*(X-x(k-1))-(y(k-1)-(M(k)*h(k).^2)/6)*(X-x(k)))/h(k));endﻫS=S.';ﻫdisp(S);ﻫfprintf('以上为样条函数(三弯矩)解析式,显示为手写如下:\n');ﻫpretty(S);%第二部分%是否继续运行程序ﻫmyloop=input('继续运行程序输入“1”,否则输入“0”\n');ﻫifmyloopwhile myloopxi=input('输入需要计算的点的值,并按回车键\n');if xi>x0|xi<x(n)ﻫfprintf('现在开始计算输入点的插值函数值……\n');elseﻫfprintf('输入数值不在插值范围内,请重新输入\n');ﻫxi=input('输入需要计算的点的值,并按回车键……\n');end%确定输入的数值应该使用哪个解析式newx=[x0;x];[r,suoy]=min(abs(newx-xi));ﻫfprintf('输入点的插值函数值为:\n\n');fpr intf('\t');if xi<=newx(suoy)ﻫf=subs(S(suoy-1),X,xi);ﻫelsef=subs(S(suoy),X,xi);enddisp(f);ﻫmyloop=input('继续计算输入“1”,终止计算输入“0”\n');endelsereturn;end%%%%%%%%%%%%function[x]=chase(a,b,c,d)%追赶法解性方程组 a是下三角b是对角线c是上三角 d是常数项%输入的a bc d 均为列向量ﻫn=length(b);ﻫu=zeros(n,1);ﻫv=zeros(n,1);ﻫx=zeros(n,1);%追v(1)=c(1)/b(1);u(1)=d(1)/b(1);for i=2:n-1v(i)=c(i)/(b(i)-v(i-1)*a(i-1));u(i)=(d(i)-u(i-1)*a(i-1))/(b(i)-v(i-1)*a(i-1));endﻫu(n)=(d(n)-u(n-1)*a(n-1))/(b(n)-v(n-1)*a(n-1));ﻫ%赶ﻫx (n)=u(n);ﻫfori=n-1:-1:1x(i)=u(i)-v(i)*x(i+1);ﻫend。

三次样条插值端点约束条件的构造与Matlab实现

三次样条插值端点约束条件的构造与Matlab实现

三次样条插值端点约束条件的构造与Matlab实现邢丽【摘要】Spline interpolation techniques are increasingly important in engineering calculations. The boundary conditions of the cubic spline interpolation are given according to the actual problem in the state of the endpoint. Through researching cubic spline function interpolation constraints for different endpoints, using Matlab computational analysis, each interval segment cubic spline function body expression is showed. The point of interpolation is calculated and each interval segment graph is displayed which is applied to practical problems. Endpoint constraints as well as mixed boundary conditions is focused on.% 在工程计算中,样条插值技术的研究越来越重要。

三次样条插值的边界条件是根据实际问题在端点的状态给出。

通过研究三次样条函数插值,针对不同的端点约束,用 Matlab 计算分析,显示各区间段三次样条函数体表达式,计算出已给点插值并显示各区间分段曲线图,并应用到实际问题中。

重点讨论端点约束条件以及混合边界条件。

【期刊名称】《上海第二工业大学学报》【年(卷),期】2012(000)004【总页数】5页(P319-323)【关键词】计算数学;三次样条插值;端点约束;Matlab【作者】邢丽【作者单位】上海第二工业大学理学院,上海201209【正文语种】中文【中图分类】P315.31在工程计算中,插值技术的研究越来越重要。

三次样条插值函数的构造与Matlab实现

三次样条插值函数的构造与Matlab实现

自动测量与控制 Automatic Measurement and Control
O. I. Automation 2006, Vol. 25, No. 11
三次样条插值函数的构造与 Matlab 实现
许小勇 1 ,钟太勇 1,2 ( 1. 云南民族大学 数学与计算机科学学院, 云南 昆明 650031 ; 2. 郧阳师范高等专科学校 数学系, 湖北 丹江口 442700 ) 摘要: 三次样条插值函数边界条件由实际问题对三次样条插值在端点的状态要求给出。以第 1 边界条件为例, 用节点处二阶导数表示三次样条插值函数,用追赶法求解相关方程组。通过 Matlab 编制三次样条函数的通用程序, 可直接显示各区间段三次样条函数体表达式,计算出已给点插值并显示各区间分段曲线图。 关键词: 三次样条;插值函数; Matlab 程序 中图分类号: O242.1 文献标识码: A
注意到 S(x) 在 [x j, x j+1 ]( j=1,2,… ,n- 1 )上是三 次多项式,于是 S"(x)在 [x j, x j+1 ] 上是一次多项式, 如果 S"(x) 在 [x j,x j+1 ]( j=1,2,… ,n -1)两端点上的值 已知,设 S"(x j)=M j,S"(x j+1 )=M j+1 ,则 S"(x) 的表达 x j+1 − x x −xj Mj + M j+1 , 其 中 h j = 式 为 : S'' ( x ) = hj hj x j+1 -x j,对 S"(x) 进行两次积分,则得到 1 个具有 2
S j (x) = a j x 3 + b j x 2 + c j x + d j , (j = 1,2, … ,n - 1) (1)

MATLAB 三次样条

MATLAB 三次样条

12.1
基本特征
在三次样条中,要寻找三次多项式,以逼近每对数据点间的曲线。在样条术语中,这 些数据点称之为断点。因为,两点只能决定一条直线,而在两点间的曲线可用无限多的三 次多项式近似。因此,为使结果具有唯一性。在三次样条中,增加了三次多项式的约束条 件。通过限定每个三次多项式的一阶和二阶导数,使其在断点处相等,就可以较好地确定 所有内部三次多项式。此外,近似多项式通过这些断点的斜率和曲率是连续的。然而,第 一个和最后一个三次多项式在第一个和最后一个断点以外,没有伴随多项式。因此必须通 过其它方法确定其余的约束。最常用的方法,也是函数 spline 所采用的方法,就是采用非 扭结(not-a-knot)条件。这个条件强迫第一个和第二个三次多项式的三阶导数相等。对最后 一个和倒数第二个三次多项式也做同样地处理。 基于上述描述,人们可能猜想到,寻找三次样条多项式需要求解大量的线性方程。实 际上,给定 N 个断点,就要寻找 N-1 个三次多项式,每个多项式有 4 个未知系数。这样, 所求解的方程组包含有 4*(N-1)个未知数。把每个三次多项式列成特殊形式,并且运用各种 约束,通过求解 N 个具有 N 个未知系数的方程组,就能确定三次多项式。这样,如果有 50 个断点,就有 50 个具有 50 个未知系数的方程组。幸好,用稀疏矩阵,这些方程式能够简 明地列出并求解,这就是函数 spline 所使用的计算未知系数的方法。
0 7.0000 0.0007 -0.0083 0.0042 0.3542 0.1635 4.9136 0.9391
1.0000 8.0000 0.0007 0.1068 0.0072 -0.2406 0.1925 0 1.2088
2.0000 9.0000 0.0010 -0.1982 0.0109 4.2439 0.2344 0.1263 1.5757

Matlab Spline 三次样条插值多项式表达式问题

Matlab Spline 三次样条插值多项式表达式问题

如何运用MATLAB 三次样条插值的问题,今天做作业,突然想用Matlab搞搞。

题目如下:清华大学出版社的《数值分析(第5版)》P49,20题。

x=[0.25 0.3 0.39 0.45 0.53];y=[ 0.5 0.5477 0.6245 0.6708 0.7280 ]pp=csape(x,y,'second',[0,0.0]);disp(pp.coefs);其中COEFS的含义是在Xi-Xi+1区间上的多项式是,例如COEFS数组第一行的意思是在X=0.25到X=0.3的区间上时表达式是-6.2652*(X-0.25)^3+0.9697*(X-0.25)^1+0.5;-6.2652 0.0000 0.9697 0.50001.8813 -0.9398 0.9227 0.5477-0.4600 -0.4318 0.7992 0.62452.1442 -0.5146 0.7424 0.6708关于csape的用法引用自:/ck436/blog/item/6fe40c46400d3c046b63e52b.htmlcsape,是计算在各种边界条件下的三次样条插值。

pp = csape(x,y,conds)其中conds主要有以下的选项variational(自然边界条件,首末点二阶导数均为0),second (指定首末点的二阶导数),periodic(周期性边界条件,首末点的0~2阶导数相等),complete (给定导数情况,默认)function pp = csape(x,y,conds,valconds)%pp=csape(x,y,'变界类型','边界值'),生成各种边界条件的三次样条插值. 其中,(x,y)为数据向量%边界类型可为:'complete',给定边界一阶导数.% 'not-a-knot',非扭结条件,不用给边界值.% 'periodic',周期性边界条件,不用给边界值.% 'second',给定边界二阶导数.% 'variational',自然样条(边界二阶导数为0)% .%例考虑数据% x | 1 2 4 5% ---|-------------% y | 1 3 4 2%边界条件S''(1)=2.5,S''(5)=-3,% x=[1 2 4 5];y=[1 3 4 2];% pp=csape(x,y,'second',[2.5,-3]);pp.coefs % xi=1:0.1:5;yi=ppval(pp,xi);% plot(x,y,'o',xi,yi);。

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

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

MATLAB作业给定一个时间序列,使用三次样条插值方法进行均匀内插(题目的相关说明:按题目要求编写一个MATLAB程序函数,并把自己编制程序所得的结果与MATLAB库函数分析结果进行对比。

)理论基础:时间序列的概念:时间序列是一种定量预测方法,又称简单外延法,时间序列分析是根据系统观测得到的时间序列数据,通过曲线拟合和参数估计来建立数学模型的理论与方法,时间序列分析可分为以下三种情况(1)把一个时间序列的值变动为N 个组成部分,通常可以分为四种 a、倾向变动,又称长期趋势变动 b、循环变动,又称周期变动 c、季节变动,即每年有规则的反复进行变动 d、不规则变动,即随机变动。

然后把这四个综合到一起得出预测的结果。

虽然分成这四部分,但这四部分之间的相互关系是怎么样的呢,目前一般采用相乘的关系,其实各个部分都是在其他部分作用的基础上进行作用的,所以采用相乘是有一定依据的,此种方法适合于短期预测和库存预测(2)把预测对象、预测目标和对预测的影响因素都看成为具有时序的,为时间的函数,而时间序列法就是研究预测对象自身变化过程及发展趋势,如果未来预测是线性的,其数学模型为YT+L=aT+bTL,YT+L为未来预测值,aT为截距,bT为斜率,L为由T到需要预测的单位时间数(如5年、10年等)(3)根据预测对象与影响因素之间的关系及影响程度来推算未来,与目标的相关因素很多,只能选择那些因果关系较强的为预测影响的因素,此时间序列法用于短期预测比较有效,若要用于长期预测,还需要结合其他方法才行。

三次样条插值的实际应用:在制造船体和汽车外形等工艺中传统的设计方法是,首先由设计人员按外形要求,给出外形曲线的一组离散点值,施工人员准备好有弹性的样条(一般用竹条或有弹性的钢条)和压铁,将压铁放在点的位置上,调整竹条的形状,使其自然光滑,这时竹条表示一条插值曲线,我们称为样条函数。

从数学上看,这一条近似于分段的三次多项式,在节点处具有一阶和二阶连续微商。

用matlab实现3次样条曲线插值的算法程序

用matlab实现3次样条曲线插值的算法程序

注意:采用的是Lagrange插值的第二种情况实验二Lagrange插值【实验目的】1.了解插值法及Lagrange插值的基本概念.2.学习、掌握MATLAB软件有关的命令。

【实验原理】插值法定义:设函数y=f(x)在区间[a,b]上有定义,且已知f(x)在[a,b]上n+1个互异点a x0<x1<…<xn b处的值yi=f(xi),i=0,1,2,…,n.若存在一个简单函数P(x),使P(xi)=yi (i=0,1,2,…,n)(2.1)成立,则称P(x)为f(x)的插值函数,f(x)称为被插函数,点xi(i=0,1,2,…,n)称为插值节点,[a,b]称为插值区间,(2.1)式为插值条件。

求插值函数的方法为插值法。

利用n次插值基函数可以将满足条件Ln(xn)=yk(k=0,1,2,…,n)的插值多项式Ln(x)表示为yk (2.2)称Ln(x)为拉格朗日插值多项式【实验内容】在区间[-5,5]上取截点n=11,等距间隔h=1的节点为插值节点,对函数f(x)=1/(1+x2)进行拉格朗日插值,并绘图。

下面的程序基于公式(2.2),且在xi点,Ln(xi)= yk相应的Matlab代码为function yi=lagrange(x,y,xi)%Lagrange插值%x,y为已知节点及其函数值向量%xi为插值点(可以是多个),yi为插值n=length(x);m=length(xi);%length是x的长度,for i=1:mz=xi(i);s=0;for k=1:n%for循环语句p=1.0;for j=1:nif j~=kp=p*(z-x(j))/(x(k)-x(j));endends=p*y(k)+s;endyi(i)=s;end下面用MATLAB运行程序得到以下结果,如下图:用Matlab实现了3次样条曲线插值的算法。

边界条件取为自然边界条件,即:两个端点处的2阶导数等于0;共包含3各个函数文件,主函数所在文件(即使用的时候直接调用的函数)为spline3.m,另外两个函数文件是在splin3函数文件中被调用的自定义函数。

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

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

MATLAB作业给定一个时间序列,使用三次样条插值方法进行均匀内插(题目的相关说明:按题目要求编写一个MATLAB程序函数,并把自己编制程序所得的结果与MATLAB库函数分析结果进行对比。

)理论基础:时间序列的概念:时间序列是一种定量预测方法,又称简单外延法,时间序列分析是根据系统观测得到的时间序列数据,通过曲线拟合和参数估计来建立数学模型的理论与方法,时间序列分析可分为以下三种情况(1)把一个时间序列的值变动为N 个组成部分,通常可以分为四种 a、倾向变动,又称长期趋势变动 b、循环变动,又称周期变动 c、季节变动,即每年有规则的反复进行变动 d、不规则变动,即随机变动。

然后把这四个综合到一起得出预测的结果。

虽然分成这四部分,但这四部分之间的相互关系是怎么样的呢,目前一般采用相乘的关系,其实各个部分都是在其他部分作用的基础上进行作用的,所以采用相乘是有一定依据的,此种方法适合于短期预测和库存预测(2)把预测对象、预测目标和对预测的影响因素都看成为具有时序的,为时间的函数,而时间序列法就是研究预测对象自身变化过程及发展趋势,如果未来预测是线性的,其数学模型为YT+L=aT+bTL,YT+L为未来预测值,aT为截距,bT为斜率,L为由T到需要预测的单位时间数(如5年、10年等)(3)根据预测对象与影响因素之间的关系及影响程度来推算未来,与目标的相关因素很多,只能选择那些因果关系较强的为预测影响的因素,此时间序列法用于短期预测比较有效,若要用于长期预测,还需要结合其他方法才行。

三次样条插值的实际应用:在制造船体和汽车外形等工艺中传统的设计方法是,首先由设计人员按外形要求,给出外形曲线的一组离散点值,施工人员准备好有弹性的样条(一般用竹条或有弹性的钢条)和压铁,将压铁放在点的位置上,调整竹条的形状,使其自然光滑,这时竹条表示一条插值曲线,我们称为样条函数。

从数学上看,这一条近似于分段的三次多项式,在节点处具有一阶和二阶连续微商。

三次样条插值matlab代码实现

三次样条插值matlab代码实现

三次样条插值matlab代码实现三次样条插值是一种常用的插值方法,可以用于曲线拟合和数据逼近。

在Matlab中,可以使用内置函数`interp1`来实现三次样条插值。

下面是一个简单的示例代码,演示了如何在Matlab中实现三次样条插值:matlab.% 创建一些示例数据。

x = 1:5;y = [3 6 5 8 2];% 生成更密集的x值,用于插值。

xi = 1:0.1:5;% 使用interp1进行三次样条插值。

yi = interp1(x, y, xi, 'spline');% 绘制原始数据和插值结果。

plot(x, y, 'o', xi, yi, '-');legend('原始数据', '三次样条插值');在这个示例中,我们首先创建了一些示例数据`x`和`y`,然后生成了更密集的`xi`值,用于插值。

接下来,我们使用`interp1`函数进行三次样条插值,并将结果存储在`yi`中。

最后,我们使用`plot`函数将原始数据和插值结果可视化出来。

需要注意的是,`interp1`函数中的第四个参数'spline'表示我们使用三次样条插值方法。

除了'spline'外,还可以选择'linear'(线性插值)或'pchip'(分段立方插值)等方法,具体选择取决于实际情况和数据特点。

以上就是在Matlab中实现三次样条插值的简单示例代码。

当然,实际应用中可能涉及到更复杂的数据和情况,需要根据具体问题进行相应的调整和处理。

希望这个示例能够帮助到你理解如何在Matlab中实现三次样条插值。

数值分析实验报告-插值、三次样条

数值分析实验报告-插值、三次样条

实验报告:牛顿差值多项式&三次样条... . (1)问题:在区间[-1,1]上分别取n=10、20用两组等距节点对龙格函数f (x)---作多项式插25 x 2值及三次样条插值对每个n值,分别画出插值函数矽(x)的图形。

实验目的:通过编程实现牛顿插值方法和三次样条方法,加深对多项式插值的理解。

应用所编程序解决实际算例。

实验要求:1.认真分析问题,深刻理解相关理论知识并能熟练应用;2.编写相关程序并进行实验;3.调试程序,得到最终结果;4.分析解释实验结果;5.按照要求完成实验报告。

实验原理:详见《数值分析第5版》第二章相关容。

实验容:(1)牛顿插值多项式1.1 当 n=10 时:在Matlab下编写代码完成计算和画图。

结果如下:代码:clear allclcx1=-1:0.2:1;y1=1./(1+25.*x1.八2);n=length(x1);f=y1(:);for j=2:nfor i=n:-1:jf(i) = (f(i)-f(i-1))/(x1(i)-x1(i-j+1));endendsyms F x p;F(1)=1;p(1)=y1(1);for i=2:nF(i)=F(i-1)*(x-x1(i-1));p(i)=f(i)*F(i);endsyms PP=sum(p);P10=vpa(expand(P),5);x0=-1:0.001:1;y0=subs(P,x,x0);y2=subs(1/(1+25火x八2),x,x0);plot(x0,y0,x0,y2)grid onxlabel('x')ylabel('y')P10即我们所求的牛顿插值多项式,其结果为:P10(x )=-220.94*x A10+494.91*x A8-9.5065e-14*x A7-381.43*x A6-8.504e-14*x A5+123.36*x A4+2.0202e-14*x A3-16.855*x A2-6.6594e-16*x+1.0并且这里也能得到该牛顿插值多项式的在[-1,1]上的图形,并和原函数进行对比(见Fig.1)。

分段低次插值及三次样条函数插值MATLAB实验报告

分段低次插值及三次样条函数插值MATLAB实验报告

分段低次插值及三次样条函数插值MATLAB实验报告1、 问题描述:自己编写程序或者用interp1(),及样条函数来实现分段低次插值及三次样条函数插值。

利用插值函数的图像来分析插值效果。

2、 实验步骤(过程):clc,clear;syms t;x=0:0.1:3;y=1/(1+t^2);df=diff(y);n=length(x);for i=1:ny(i)=1/(1+x(i)^2);m(i)=subs(df,t,x(i));endfor h=1:n-1x0(h)=(x(h)+x(h+1))/2;a(h)=fix(x0(h));b(h)=ceil(x0(h));k(h)=find(x==a(h));w(h)=find(x==b(h));I(h)=(x0(h)-b(h))/(a(h)-b(h))*y(k(h))+(x0(h)-a(h))/(b(h)-a(h))*y(w(h));S(h)=(1+2*(x0(h)-a(h)))*(x0(h)-b(h))^2*y(k(h))+(1-2*(x0(h)-b(h)))*(x0(h)-a(h))^2*y(w(h))+(x0(h)-a(h))*(x0(h)-b(h))^2*m(k(h))+(x0(h)-b(h))*(x0(h)-a(h))^2*m(w(h));endI=eval(I);S=eval(S);y=eval(y);plot(x,y,'k')hold onplot(x0,I,'r',x0,S,'b')hold onZ=interp1(x,y,x0,'spline');plot(x0,Z,'g')3、 结论:由图像可以得出相应的结论:分段线性插值的插值效果没有3次Hermite插值和3次样条插值的效果好,且3次Hermite插值和3次样条插值与原函数的误差相对较小,插值效果好。

当x值取值越大时,这三种方法的插值效果并不会有太大区别。

自编的三次样条插值matlab程序(含多种边界条件)

自编的三次样条插值matlab程序(含多种边界条件)
-5.6 -5.6 -5.7 -5.7 -5.8 -5.8 -5.9 -5.9 -6 -6 -6.1 -6.1 -6.2 -6.2 -6.3 -6.3 -6.4 5.5 -6.4 5.5 6 6.5 6.5 7 7.5 7.5 8 8.5 8.5 9
d2 0.67221 0.43138 0.28102 0.22141 0.24664 0.2551 0.43028
0.9802
则旋转后的三次样条的系数及图像为:
xx2=[x2(1):0.001:x2(end)]'; [yy2 b2 c2 d2]=spline3(x2,y2,xx2,1,v1,vn); fprintf('\t\t\tb2\t\t\tc2\t\t\td2\n'); disp([b2 c2(1:end-1,1) d2]); plot(x2,y2,'*b',xx2,yy2,'-.k'); grid on; b2 c2 -0.97849 -0.74704 -0.35362 -0.067629 0.0061747 0.3081 0.6992
" 又 设 cn Sn 1 xn 2 , 记 i xi 1 xi , i yi 1 yi , i 1, 2, , n 1 , 则由 (1.3)可 得 :
ci 1 ci , i 1, 2, , n 1. 3 i 从(1.2)解得: bi i ci i di i2 i i 2ci ci 1 , i 1, 2, , n 1. i i 3 将(1.4)与(1.5)代入(1.3)得: di


ai yi , i 1, 2, n 1. 2 3 y2 y1 b1 x2 x1 c1 x2 x1 d1 x2 x1 , 2 3 y y b x x n 1 n 1 n n 1 cn 1 xn xn 1 d n 1 xn xn 1 n 由节点处的一阶与二阶光滑性可知: Si'1 xi Si' xi , Si"1 xi Si" xi , i 1, 2, , n.

三次样条插值的MATLAB实现

三次样条插值的MATLAB实现

MATLAB 程序设计期中考查在许多问题中,通常根据实验、观测或经验得到的函数表或离散点上的信息,去研究分析函数的有关特性。

其中插值法是一种最基本的方法,以下给出最基本的插值问题——三次样条插值的基本提法:对插值区间[]b a ,进行划分:b x x x a n ≤<⋯⋯<<≤10,函数()x f y =在节点i x 上的值()()n i x f y i i ⋯⋯==,2,1,0,并且如果函数()x S 在每个小区间[]1,+i i x x 上是三次多项式,于[]b a ,上有二阶连续导数,则称()x S 是[]b a ,上的三次样条函数,如果()x S 在节点i x 上还满足条件()()n i y x S i i ⋯⋯==,1,0则称()x S 为三次样条插值函数。

三次样条插值问题提法:对[]b a ,上给定的数表如下.求一个分段三次多项式函数()x S 满足插值条件()()n i y x S i i ⋯⋯==,1,0 式,并在插值区间[]b a ,上有二阶连续导数。

这就需要推导三次样条插值公式:记()x f '在节点i x 处的值为()i i m x f ='(n i ⋯⋯=,1,0)(这不是给定插值问题数表中的已知值)。

在每个小区间[]1,+i i x x 利用三次Hermite 插值公式,得三次插值公式:()()()()1111+++++++=i i i i i i i i i m m x y x y x x S ββαα,[]1,+∈i i x x x 。

为了得到这个公式需要n 4个条件:(1).非端点处的界点有n 2个;(2).一阶导数连续有1-n 个条件;(3).二阶导数连续有1-n 个条件,其中边界条件:○1()()n n m x S m x S ='=' 00 ○2()()αα=''=''n x S x S 00 ○3()()()()16500403βααβαα=''+'=''+'n n x S x S x S x S○4n y y =0 ()()()()000000-''=+''-'=+'n nx S x S x S x S 其中:()⎩⎨⎧=≠=j i j i x j i,1,0α ()0='j i x α ()0=j i x β 且(1,0,=j i )。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

MATLAB 程序设计期中考查
在许多问题中,通常根据实验、观测或经验得到的函数表或离散点上的信息,去研究分析函数的有关特性。

其中插值法是一种最基本的方法,以下给出最基本的插值问题——三次样条插值的基本提法:
对插值区间[]b a ,进行划分:b x x x a n ≤<⋯⋯<<≤10,函数()x f y =在节点
i x 上的值()()n i x f y i i ⋯⋯==,2,1,0,并且如果函数()x S 在每个小区间[]1,+i i x x 上
是三次多项式,于[]b a ,上有二阶连续导数,则称()x S 是[]b a ,上的三次样条函数,如果()x S 在节点i x 上还满足条件 则称()x S 为三次样条插值函数。

三次样条插值问题提法:对[]b a ,上给定的数表如下.
求一个分段三次多项式函数()x S 满足插值条件()()n i y x S i i ⋯⋯==,1,0 式,并在
插值区间[]b a ,上有二阶连续导数。

这就需要推导三次样条插值公式: 记()x f '在节点i x 处的值为()i i m x f ='(n i ⋯⋯=,1,0)(这不是给定插值问题数表中的已知值)。

在每个小区间[]1,+i i x x 利用三次Hermite 插值公式,得三次插值公式:
()()()()1111+++++++=i i i i i i i i i m m x y x y x x S ββαα,[]1,+∈i i x x x 。

为了得到这个公式需要n 4个条件:
(1).非端点处的界点有n 2个;(2).一阶导数连续有1-n 个条件;(3).二阶导数
连续有1-n 个条件,其中边界条件:○1()()n
n m x S m x S ='=' 00
其中:()⎩⎨⎧=≠=j i j i x j i
,1,0α ()0='j i x α ()0=j i x β 且(1,0,=j i )。

()⎩⎨
⎧=≠='j i j i x j i
,1,0β ,i m 为对应变量的一阶导数。

其推导过程如下:
为了确定i m 的值,把()x S 展开为: +
()()
()()
,12
122
21+++--+
--i i
i i i i
i i m h x x x x m h x x x x
这里i i i x x h -=+1,对()x S 连续求两次导,得:
()()
()i i i
i i i i
i i i i
i i y y h x x x m h x x x m h x x x x S --++
--+
--=
''+++++13
112
1
2
1
26246426。

于是
考虑()x S ''在节点i x 处的右极限值,得: ()()i i i
i i i i i y y h m h m h x S -+--
=+''++1216
240。

同理,在相邻小区间[]i i x x ,1-上可得()x S ''的表达式为: 及()x S ''在节点i x 处的左极限值为:
()()12
1
1116
420-------+=
-''i i i i i i i i y y h m h m h x S 。

利用()x S 二阶导数于节点i x 处的连续性条件()()00-''=+''i i x S x S ,这里1,2,1-⋯⋯=n i ,有下式成立:
⎪⎪⎭⎫ ⎝⎛-+-=+⎪⎪⎭⎫ ⎝⎛++--++---211211111311121i i i i i i i i i i i i i h y y h y y m h m h h m h ,用i i h h 111+-除等式两边,并注意[]11,,++=-=i i i
i
i i i x x f h y y f y
,上式可简记为: 且[][]()1111
1,,31+----+=+=-=+=
i i i i i i i i
i i i i i i i i x x f x x f g h h h h h h μλλμλ
最后求得n m m ⋯⋯1的线性方程组为:
⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡⋯=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⋯⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯----n n n n n n n n g g g g m m m m 1211211122112000200000020002 λμμλμλλμ (**) 通过以上复杂的求解和迭代,就可以求解出插值函数的近似表达式。

得出来的表达式就可以用MATLAB 软件来求解。

具体求解过程如下:
已知n 对数据点()()()(),,,,,,,,332211n n y x y x y x y x ⋯⋯,假设函数关系为
()x f y =,但解析式不确定,数据插值就是构造函数关系式()x g y =,使()n i x i ,,3,2,1⋯⋯=,满足关系()()i i x f x g =。

例题:求满足下面函数表所给出的插值条件的三次自然样条函数。

分析:表中所列出的是函数对点,首先要把对应的插值函数求出来,再用
MATLAB 软件来求区间[]5,1上间隔为0.5的各点的值。

求解过程如下:
因自然样条插值函数的边界条件为
这里3=n ,故确定3210,,,m m m m 的方程组形式形如上面的(**)式,其中系数i i μλ,
和i g 可按如下步骤进行:
将上述参数带入(**)式,得到以下方程组: 解得: 由公式 +()()
()()
,12
122
21+++--+
--i i
i i i i
i i m h x x x x m h x x x x
可知,
由所求出的表达式可知区间[]5,1可分为[][]5,44,1⋃,对两个区间分别用MATLAB 命令即可: 针对第一个区间:
147
828123-++-=x x x y ; 其图像如下
命令如下:
x=1:4;y=(-1/8)*x.^3+(2/8)*x.^2+(7/4)*x-1;xi=1:0.5:4; y1=interp1(x,y,xi,'spline') 其运行结果如下: y1 =
Columns 1 through 6
0.8750 1.7656 2.5000 2.9844 3.1250 2.8281 Column 7
2.0000 针对第二个区间:。

相关文档
最新文档