三次样条插值多项式matlab
- 格式:doc
- 大小:60.00 KB
- 文档页数:6
三次样条插值函数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 ----------⨯⎛⎫⎪ ⎪ ⎪ ⎪⎪⎪ ⎪-- ⎪-+- ⎪--⎝⎭为了使样条插值问题有惟一解,我们在原有方程基础上增加两个边界条件。
样条插值函数 matlab样条插值函数是一种广泛应用于数值计算中的插值技术。
它的主要功能是使用已经给定的数据点来构建一个连续且平滑的函数,该函数可以通过插值来预测在给定数据点之间的值。
在 MATLAB 中,我们可以使用 spline() 函数来实现样条插值函数的操作。
样条插值函数的主要思想是在给定的数据点上构建一组多项式,以实现对在数据点之间位置的函数的插值。
这些多项式通常在相邻数据点之间的区间上定义,并满足一系列的平滑条件。
在实践中,可以使用不同的插值多项式来构建样条插值函数,例如线性样条,二次样条或三次样条。
为了使用 MATLAB 中的 spline() 函数进行样条插值的操作,我们需要先确定要插值的数据点以及在这些数据点上的值。
通常我们会使用一维数组来存储这些数据点的值,例如:x = [0 0.5 1 1.5 2];接下来,我们需要为这些数据点上的值生成一个一维数组,例如:y = [0 0.4794 0.8415 0.9975 0.9093];为了使用 spline() 函数进行样条插值,我们需要将这些数据点作为输入参数传递给该函数,例如:yy = spline(x,y);这将返回一个样条插值函数,可以使用该函数来计算在数据点之间位置的函数值。
例如,要计算在 0.25 处的函数值,我们可以使用以下代码:val = ppval(yy,0.25);上述代码将返回插值函数在 0.25 处的函数值。
我们还可以使用 plot() 函数来可视化样条插值函数的结果,例如:plot(x,y,'o',0:0.01:2,ppval(yy,0:0.01:2));该代码将绘制原始数据点以及样条插值函数在给定范围内的函数值。
样条插值函数是一种强大的数值计算技术,可以用于许多科学和工程应用。
在 MATLAB 中,我们可以使用 spline() 函数轻松实现样条插值函数的操作并可视化结果。
例5.6.1 给定以下数据, 求出三次样条函数,并计算函数分别在-0.15,-0.05,0.05,0.18,0.25处的近似值,并作图。
解:编程如下:clearx=[0.1,0.2,0.15,0,-0.2,0.3];y=[0.95,0.84,0.86,1.06,1.50,0.72];pp=spline(x,y);pp.coefsxx=[-0.15,-0.05,0.05,0.18,0.25];yy=ppval(pp,xx) %or:yy=spline(x,y,xx)fnplt(pp,'k')hold onplot(x,y,'o')hold onplot(xx,yy,'r*')运行结果:ans =-36.3850 21.8592 -5.1164 1.5000-36.3850 0.0282 -0.7390 1.0600227.6995 -10.8873 -1.8249 0.9500-143.0047 23.2676 -1.2059 0.8600-143.0047 1.8169 0.0484 0.8400yy =1.2943 1.1016 1.0186 0.8409 0.8291写出三次样条函数为:⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤+-+-+--≤≤+---+--≤≤+-----≤≤+-+-≤≤-++-+++-=3.02.0,84.0)2.0(0484.0)2.0(8169.1)2.0(0047.1432.015.0,86.0)15.0(2059.1)15.0(2676.23)15.0(0047.14315.01.0,95.0)1.0(8249.1)1.0(8873.10)1.0(6995.2271.00,06.1739.00282.0385.3602.0,5.1)2.0(1164.5)2.0(859.21)2.0(385.36)(2323232323x x x x x x x x x x x x x x x x x x x x x s练习: 给定函数11)(2+=x x f ,55≤≤-x ,节点)10,,1,0(5 =+-=k k x k ,用三次样条插值求)(x S 及在某些点处的近似值.MATLAB 实现:yi=interp1(x,y,xi) %根据数据(x,y)给出分段线性插值函数的xi 的值yi.yi=interp1(x,y,xi,’spline’) %使用三次样条插值.yi=interp1(x,y,xi,’cubic’) %使用分段三次插值.Yi=spline(x,y,xi) %同yi=interp1(x,y,xi,’spline’)pp=spline(x,y) %返回三次样条插值的分段多项式(pp)形式结构(使用“非扭结”端点条件).MATLAB 中的spline 函数可以对数据点进行三次插值,默认的边界条件强制使得插值第一段的三次项系数与第二段的三次项系数相同,最后一段与倒数第二段的多项式系数相同,这就所谓的非扭结(not-a-knot )条件. 1一维插值函数Interp1()命令格式: yi=interp1(x,y,xi,’method’)x 为插值节点构成的向量,y 为插值节点函数值构成的向量,yi 是被插值点xi 的插值结果,‘method ‘是采用的插值方法,缺省时表示分线段性插值,’nearest ‘为最邻近插值;’linear ‘为分线段性插值;’spline’为三次样条插值;’pchip’为分段Hermite 插值;’cubic’为分段Hermite 插值例子:画出y=sin(x)在区间[0 10]的曲线,并在曲线上插值节点xk=k,k=0,1…10及函数值,画出分段线性插值折线图x=0:10;y=sin(x);xi=0:0.25:10;yi1=interp1(x,y,xi,'nearest');yi2=interp1(x,y,xi,'linear');yi3=interp1(x,y,xi,'spline');yi4=interp1(x,y,xi,'pchip');yi5=interp1(x,y,xi,'cubic');subplot(1,5,1)plot(x,y,'o',xi,yi1,'k--',xi,sin(xi),'k:');title('\bfNearest');subplot(1,5,2)plot(x,y,'o',xi,yi2,'k--',xi,sin(xi),'k:');title('\bfLinear');subplot(1,5,3)plot(x,y,'o',xi,yi3,'k--',xi,sin(xi),'k:');title('\bfSpline');subplot(1,5,4)plot(x,y,'o',xi,yi4,'k--',xi,sin(xi),'k:');title('\bfPchip');subplot(1,5,1)plot(x,y,'o',xi,yi5,'k--',xi,sin(xi),'k:');title('\bfCubic');spline()为三次样条函数命令格式1:yi=spline(x,y,xi),意义等同于yi=interp1(x,y,xi,'spline')命令格式2:pp=spline(x,y) ,输出三次样条函数分段表示的结构pchip()命令格式与spline()完全相同csape()为可输入边界条件的三次样条函数命令格式:pp=csape(x,y,conds,valconds),x为插值节点构成的向量,y为插值节点函数值构成的向量;conds为边界类型,缺省为非扭结边界条件;valconds表示边界值。
MATLAB 中的曲线拟合和插值在大量的使用领域中,人们经常面临用一个分析函数描述数据(通常是测量值)的任务。
对这个问题有两种方法。
在插值法里,数据假定是正确的,要求以某种方法描述数据点之间所发生的情况。
这种方法在下一节讨论。
这里讨论的方法是曲线拟合或回归。
人们设法找出某条光滑曲线,它最佳地拟合数据,但不必要经过任何数据点。
图11.1说明了这两种方法。
标有'o'的是数据点;连接数据点的实线描绘了线性内插,虚线是数据的最佳拟合。
11.1 曲线拟合曲线拟合涉及回答两个基本问题:最佳拟合意味着什么?应该用什么样的曲线?可用许多不同的方法定义最佳拟合,并存在无穷数目的曲线。
所以,从这里开始,我们走向何方?正如它证实的那样,当最佳拟合被解释为在数据点的最小误差平方和,且所用的曲线限定为多项式时,那么曲线拟合是相当简捷的。
数学上,称为多项式的最小二乘曲线拟合。
如果这种描述使你混淆,再研究图11.1。
虚线和标志的数据点之间的垂直距离是在该点的误差。
对各数据点距离求平方,并把平方距离全加起来,就是误差平方和。
这条虚线是使误差平方和尽可能小的曲线,即是最佳拟合。
最小二乘这个术语仅仅是使误差平方和最小00.20.40.60.81-2024681012xy =f (x )Second O rder C urv e Fitting图11.1 2阶曲线拟合在MATLAB 中,函数polyfit 求解最小二乘曲线拟合问题。
为了阐述这个函数的用法,让我们以上面图11.1中的数据开始。
» x=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1]; » y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];为了用polyfit ,我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度。
如果我们选择n=1作为阶次,得到最简单的线性近似。
三次样条众所周知,使用高阶多项式的插值常常产生病态的结果。
目前,有多种消除病态的方法。
在这些方法中,三次样条是最常用的一种。
在MATLAB中,实现基本的三次样条插值的函数有spline,ppval,mkpp和unmkpp。
在这些函数中,仅spline 在《MATLAB参考指南》中有说明。
下面几节,将展示在M文件函数中实现三次样条的基本特征基本特征在三次样条中,要寻找三次多项式,以逼近每对数据点间的曲线。
在样条术语中,这些数据点称之为断点。
因为,两点只能决定一条直线,而在两点间的曲线可用无限多的三次多项式近似。
因此,为使结果具有唯一性。
在三次样条中,增加了三次多项式的约束条件。
通过限定每个三次多项式的一阶和二阶导数,使其在断点处相等,就可以较好地确定所有内部三次多项式。
此外,近似多项式通过这些断点的斜率和曲率是连续的。
然而,第一个和最后一个三次多项式在第一个和最后一个断点以外,没有伴随多项式。
因此必须通过其它方法确定其余的约束。
最常用的方法,也是函数spline所采用的方法,就是采用非扭结(not-a-knot)条件。
这个条件强迫第一个和第二个三次多项式的三阶导数相等。
对最后一个和倒数第二个三次多项式也做同样地处理。
基于上述描述,人们可能猜想到,寻找三次样条多项式需要求解大量的线性方程。
实际上,给定N个断点,就要寻找N-1个三次多项式,每个多项式有4个未知系数。
这样,所求解的方程组包含有4*(N-1)个未知数。
把每个三次多项式列成特殊形式,并且运用各种约束,通过求解N个具有N个未知系数的方程组,就能确定三次多项式。
这样,如果有50个断点,就有50个具有50个未知系数的方程组。
幸好,用稀疏矩阵,这些方程式能够简明地列出并求解,这就是函数spline 所使用的计算未知系数的方法。
《数值分析》课外课堂大作业论文题目:基于多项式插值与三次样条插值曲线拟合的比较姓名:学号:学院:专业方向:联系方式:(QQ号)(手机号)导师姓名:完成人(亲笔)签字基于多项式插值与三次样条插值曲线拟合的比较摘要:在数值计算中经常要计算函数,当函数只在有限点集上给定函数值要包含改点集的区间上用公式给出函数的简单表达式,这就涉及在已知区间上用简单函数逼近已知复杂函数问题。
本文为了解决这类问题就采用多项式插值与三次样条插值两种插值法并利用MATLAB数值分析软件进行编程,实现相应数据的曲线拟合以获得最佳曲线模型与相应数据的曲线拟合,选出最优的插值法以解决所给数据的曲线拟合问题。
关键词:函数;多项式插值;三次样条插值;曲线拟合;MATLABAbstract:In numerical analysis ,the function value is often calculated .when the function is only given a function point set ,the simple expression of the function is given by the interval .which involves the use of a simple function to approximate the known complex function .in order to solve this problem ,we use polynomial interpolation and cubic spline interpolation tow kind of interpolation method and use MATLAB numerical analysis software to program ,to achieve the curve fitting of the corresponding date to obtain the best cure fitting ,and to choose the best interpolation method to solve the problem of curve fitting to the date.Keyword: Function ; Polynomial interpolation ; Cubic spline interpolation ; Fitting of a curve ; MATLAB前言现代科学研究中,物理量之间的相互关系通量是用函数来描述的,许多实际问题都用函数y=f(x)来表示某种内在规律的数量关系其中相当一部分函数是通过试验或观测得到的也有少量函数关系是由经典物理分析推导得到的,但许多实际问题很难用经典理论分析得出,因为虽然f(x)在某个区间[a,b]上是存在的,有的还是连续的,但往往这个f(x)并不包含我们所得函数表的所有值因此我们希望根据给定的函数表做一个即能反应函数f(x)的特行,又便于计算的简单函数p(x),用p(x)近似f(x),这样确定的p(x)就是我们希望得得到的插值函数。
Matlab中的插值和平滑方法1. 引言在数值分析和数据处理中,插值和平滑是常用的技术手段,可以用于填补数据的空缺以及降低数据中的噪声。
Matlab作为一种强大的数值计算和数据处理软件,提供了丰富的插值和平滑方法,本文将介绍其中的一些常用方法及其应用。
2. 插值方法2.1 线性插值线性插值是最简单的一种插值方法,它假设待插值函数在相邻数据点之间是线性变化的。
Matlab中提供了interp1函数实现线性插值,可以通过设定插值点的横坐标向量和已知数据点的横坐标向量,以及对应的纵坐标向量,得到插值结果。
2.2 分段插值分段插值是一种更精确的插值方法,它假设待插值函数在相邻数据点之间是分段线性变化的。
Matlab中的interp1函数也可以实现分段插值,通过指定'linear'插值方法和 'pchip'插值方法,可以得到不同的插值结果,前者得到的结果比较平滑,而后者更接近原始数据的形状。
2.3 样条插值样条插值是一种更高阶的插值方法,它假设待插值函数在相邻数据点之间是多项式变化的。
Matlab中的spline函数可以实现三次样条插值,它通过计算每个数据点处的二阶导数,得到一个以每个数据点为节点的三次多项式函数。
样条插值可以更加精确地还原数据,但也容易受到离群点的干扰。
3. 平滑方法3.1 移动平均移动平均是一种常用的平滑方法,它通过计算数据点周围一定范围内的平均值,得到平滑后的结果。
Matlab中的smoothdata函数提供了不同的平滑方法,包括简单移动平均、指数移动平均和加权移动平均等,可以根据具体需求选择适当的方法。
3.2 Savitzky-Golay滤波Savitzky-Golay滤波是一种基于最小二乘法的平滑方法,它通过拟合多项式曲线来实现数据的平滑。
Matlab中的sgolay函数可以实现Savitzky-Golay滤波,通过指定不同的拟合阶数和窗口大小,可以得到不同程度的平滑结果。
三次样条插值多项式
——计算物理实验作业四
陈万物理学2013级
主程序:
clear,clc;
format rat
x = [1,4,9,16,25,36,49,64];
y = [1,2,3,4,5,6,7,8];
f1 = ;
fn = 1/16;
[a,b,c,d,M,S] = spline(x,y,f1,fn);
子程序1:
function [a,b,c,d,M,S]=spline(x,y,f1,fn)
% 三次样条插值函数
% x是插值节点的横坐标
% y是插值节点的纵坐标
% u是插值点的横坐标
% f1是左端点的一阶导数
% fn是右端点的一阶导数
% a是三对角矩阵对角线下边一行
% b是三对角矩阵对角线
% c是三对角矩阵对角线上边一行
% S是插值点的纵坐标
n = length(x);
h = zeros(1,n-1);
deltay = zeros(1,n);
miu = zeros(1,n-1);
lamda = zeros(1,n-1);
d = zeros(1,n-1);
for j = 1:n-1
h(j) = x(j+1)-x(j);
deltay(j) = y(j+1)-y(j);
end % 得到h矩阵
for j = 2:n-1
sumh = h(j-1) + h(j);
miu(j) = h(j-1) / sumh;
lamda(j) = h(j) / sumh;
d(j) = 6*( deltay(j)/h(j)-(deltay(j-1)/h(j-1)))/sumh; end
% 根据第一类边界条件,作如下规定
lamda(1) = 1;
d(1) = 6*(deltay(1)/h(1)-f1)/h(1);
miu(1) = 1;
d(n) = 6*(fn-deltay(n-1)/h(n-1))/h(n-1);
% 输出三对角矩阵的a,b,c
a = miu;
b = 2*ones(1,n);
c = lamda;
M = chase(a,b,c,d); %调用chase函数得到M
sym u;
for j = 1:n-1
u = x(j)::x(j+1);
v = ones(size(u));
S = (M(j)*(x(j+1)*v-u).^3/(6*h(j))+M(j+1)*(u-x(j)*v).^3/(6*h(j))... +(y(j)-M(j)*h(j)^2/6)*(x(j+1)*v-u)/h(j)+(y(j+1)...
-M(j+1)*h(j)^2/6)*(u-x(j)*v)/h(j));
plot(u,S,'-k');
hold on
end
plot(x,y,'-.*r');
xlabel('x'),ylabel('y'),title('cubic spline interp');
end
子程序2:
function M = chase(a,b,c,f)
% 追赶法求解三对角矩阵方程,Ax=f
% a是对角线下边一行的元素
% b是对角线元素
% c是对角线上边一行的元素
n = length(b);
beta = ones(1,n-1);
y = ones(1,n);
M = ones(1,n);
for i = (n-1):(-1):1
a(i+1) = a(i);
end
% 将a矩阵和n对应
beta(1) = c(1)/b(1);
for i = 2:(n-1)
beta(i) = c(i)/( b(i)-a(i)*beta(i-1) );
end
y(1) = f(1)/b(1);
for i = 2:n
y(i) = (f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1)); end
M(n) = y(n);
for i = (n-1):(-1):1
M(i) = y(i)-beta(i)*M(i+1);
end
end
三次样条插值结果:
与拉格朗日插值作对比:
分析图一知,三次样条插值结果与预期结果吻合得很好,曲线平滑连续性好,在左右短点处的
小区间也吻合得很好,可以延伸到区间[a,b]外的一小段,用最邻近的小区间插值函数可以近似求得[a,b]区间外一小段范围内的函数值。
分析图一和图二,三次样条插值在x(i)和x(i+1)间隔很小的地方插值效果稍差,间隔稍大时插值效果非常好。
拉格朗日插值则不同,x(i)和x(i+1)间隔越大插值效果越差。