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

matlab三次样条插值函数
MATLAB是一个高级技术计算软件,它拥有强大的数值计算、绘图和数据分析功能。

在MATLAB中进行插值时,三次样条插值函数是一种优秀的插值方法。

三次样条插值函数是一种基于多项式的插值方法。

它通过使用离散数据集的低次多项式来拟合数据,然后利用这些多项式按照一定的规律来生成整个插值函数。

三次样条插值方法采用的多项式是二次函数,它会在两个数据点之间生成一条平滑的曲线。

在MATLAB中,三次样条插值函数为spline(x,y,xq),其中x和y分别代表已知离散数据集的x坐标和y坐标,xq为要进行插值的新数据点。

在调用spline函数时,需要将x和y从小到大进行排序。

在进行三次样条插值方法时,需要注意两个点。

第一,在插值函数两端需要加上边界条件,即插值函数在区间的两个端点处一阶导数值要相等,这使得插值函数在两端呈现出更平滑的性质。

第二,如果数据点x中有相同的值,则需要对数据进行处理以避免出现插值函数不确定的情况。

在使用MATLAB进行三次样条插值时,也可以借助MATLAB提供的其他函数,如ppval函数可以求解xq处插值函数的值,而fnplt函数可以绘制插值函数的图像。

总之,三次样条插值函数是一种十分有效的插值方法,可以通过MATLAB轻松地进行实现。

在进行插值时,需要注意边界条件和数据处理的细节问题。

用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)为已知的点的坐标。

计算方法上机作业——求三次样条插值函数的matlab程序

计算方法上机作业——求三次样条插值函数的matlab程序
25
附录 3 求三次样条插值函数的 matlab 程序 for f = 2:n-1; ly = 0; for g = 1:f-1 ly = ly+l(f,g)*yy(g); end yy(f) = D(f)-ly; end M1(n-1) = yy(n-1)/u(n-1,n-1); for rr=1:n-2 r = n-1-rr; uM1 = 0; for s=r+1:n-1 uM1 = uM1+u(r,s)*M1(s); end M1(r) = (yy(r)-uM1)/u(r,r); end M = [M1(n-1,1);M1]; end ss = 0; for t=1:n-1 S(t,1) = (M(t+1)-M(t))/(6*h(t)); S(t,2) = (M(t)*x(t+1)-M(t+1)*x(t))/(2*h(t)); S(t,3) = (M(t+1)*x(t)^2-M(t)*x(t+1)^2)/(2*h(t))+(y(t+1)-y(t))/h(t)+h(t)*(M(t)-M(t+1))/6; S(t,4) = (M(t)*x(t+1)^3-M(1)*x(t)^3)/(6*h(t))+(y(t)*x(t+1)-y(t+1)*x(t))/h(t)+h(t)*(M(t+1)* x(t)-M(t)*x(t+1))/6; for x1 = x(t):(x(t+1)-x(t))/100:x(t+1) ss = ss+1; xx(ss) = x1; SS(ss) = S(t,1)*x1^3+S(t,2)*x1^2+S(t,3)*x1+S(t,4); end end plot(xx,SS,'-k','linewidth',2); hold on plot(x,y,'*k','markersize',10); hold on xlabel('x'); ylabel('S(x)'); grid; fprintf('\n 所求的三次样条插值函数为:\n'); for uu=1:n-1 fprintf('S(x) = %10.5f*x^3+%10.5f*x^2+%10.5f*x+%10.5f, %8.4f<= x <=%8.4f\n',S(uu,1),S(uu,2),S(uu,3),S(uu,4),x(uu),x(uu+1)); end

三次样条插值函数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 ----------⨯⎛⎫⎪ ⎪ ⎪ ⎪⎪⎪ ⎪-- ⎪-+- ⎪--⎝⎭为了使样条插值问题有惟一解,我们在原有方程基础上增加两个边界条件。

三次样条函数插值

三次样条函数插值

计算方法第二次大作业01087.00='Y ,1007='Y ,画出插值曲线的图像:(2) 时针旋转坐标轴45度保持(1)中的节点和边界条件的几何关系不变,再次作三次样条函数插值,画出插值曲线的图像;(3) 比较(1)和(2)的结果,能得到什么结论?解:Matlab 程序如下:x=[8.125 8.4 9.0 9.485 9.6 9.959 10.166 10.2]; y=[0.0774 0.099 0.280 0.60 0.708 1.200 1.800 2.177]; pp=csape(x,y,'complete',[0.01087,100]); xx=8.125:0.001:10.2; yy=ppval(pp,xx); plot(x,y,'*',xx,yy); pp.coefsans =1.0e+004 *0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0002 -0.0000 0.0000 0.0000 -0.0069 0.0003 0.0002 0.0001 0.0074 -0.0021 -0.0001 0.0001 -0.0512 0.0058 0.0013 0.0001 4.2279 -0.0260 -0.0029 0.0002得到插值函数的图像如下:图1 三次样条函数插值图像(2)逆时针旋转坐标轴45°保持(1)中节点和边界条件的几何关系不变,再次作三次样条函数插值,画出插值曲线的图像。

解:设坐标逆时针旋转45°后,点的坐标为x ′,y ′,则有''cossin44sin cos 44x x y y ππππ⎡⎤⎢⎥⎛⎫⎛⎫=⎢⎥ ⎪ ⎪ ⎪⎝⎭⎢⎥⎝⎭-⎢⎥⎣⎦在新坐标系下:i 0 1 2 3 4 5 67X i 5.799973 6.009701 6.561951 7.131172 7.288857 7.890605 8.46124 8.751861 Y i-5.69051-5.86969-6.16597-6.28264-6.28759-6.19355-5.91566 -5.67312两个边界条件为:'0Y =-0.9785,'7Y = 0.9802。

matlab三次样条插值的方法

matlab三次样条插值的方法

matlab三次样条插值的方法Cubic spline interpolation is a common method used in MATLAB to approximate values between specified data points. This technique involves fitting a piecewise cubic polynomial to the data points, ensuring that the function is smooth and continuous at the knots. Through this process, the spline curve can accurately represent the overall trend of the data, making it particularly useful in various scientific and engineering applications.三次样条插值是MATLAB中常用的一种方法,用于在指定数据点之间近似数值。

这种技术涉及将分段三次多项式拟合到数据点,确保在节点处函数平滑连续。

通过这个过程,样条曲线可以准确地表示数据的总体趋势,使其特别适用于各种科学和工程应用。

One advantage of cubic spline interpolation is its ability to capture the local behavior of the data while maintaining global smoothness. This is achieved by constructing individual cubic polynomials between adjacent data points, ensuring that the interpolated curve passes through each data point without introducing significant oscillations or deviations. As a result, cubic splines provide a reliableand visually appealing way to interpolate data that may exhibit complex patterns or fluctuations.三次样条插值的一个优点是能够捕捉数据的局部行为,同时保持全局的平滑性。

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