第31讲 matlab插值计算.docx

1、第31讲 matlab插值计算插值方法晚上做一个曲线拟合，结果才开始用最小二乘法拟合时，拟合出来的东西太难看了！于是尝试用其他方法。经过一番按图索骥，终于发现做曲线拟合的话，采用插值法是比较理想的方法。尤其是样条插值，插完后线条十分光滑。方法付后，最关键的问题是求解时要积分，放这里想要的时候就可以直接过来拿，不用死去搜索啦。呵呵插值方法的Matlab实现一维数据插值MATLAB中用函数interp1来拟合一维数据，语法是YI = INTERP1(X,Y,XI,方法)其中（X， Y） 是已给的数据点， XI 是插值点，其中方法主要有 linear -线性插值，默认 pchip -逐段三次Herm

2、ite插值 spline -逐段三次样条函数插值其中最后一种插值的曲线比较平滑例：x=0:.12:1; x1=0:.02:1;%（其中x=0:.12:1表示显示的插值点，x1=0:.02:1表示插值的步长）y=(x.2-3*x+5).*exp(-5*x).*sin(x);plot(x,y,o); hold on;y1=interp1(x,y,x1,spline);plot(x1,y1,:)如果要根据样本点求函数的定积分，而函数又是比较光滑的，则可以用样条函数进行插值后再积分，在MATLAB中可以编写如下程序：function y=quadspln(x0,y0,a,b)f=inline(inte

3、rp1(x0,y0,x,spline),x,x0,y0);y=quadl(f,a,b,1e-8,x0,y0);现求sin(x)在区间0，pi上的定积分，只取5点x0=0,0.4,1,2,pi;y0=sin(x0);I=quadspln(x0,y0,0,pi)结果得到的值为 2.01905， 精确值为2求一段matlab插值程序悬赏分：20 - 解决时间：2009-12-26 19:57 已知5个数据点：x=0.25 0.5 0.75 1 y=0 0.3104 0.6177 0.7886 1 ，求一段matlab插值程序，求过这5个数据点的插值多项式，并在x-y坐标中画出y=f（x）图形，并且求

4、出f（x）与x轴围成图形的面积（积分），不胜感激！ 使用Lagrange 插值多项式的方法：首先把下面的代码复制到M文件中，保存成lagranfunction C,L=lagran(X,Y)% input - X is a vector that contains a list of abscissas% - Y is a vector that contains a list of ordinates% output - C is a matrix that contains the coefficients of the lagrange interpolatory polynomial%

5、- L is a matrix that contains the lagrange coefficients polynomialw=length(X);n=w-1;L=zeros(w,w);for k=1:n+1V=1;for j=1:n+1if k=jV=conv(V,poly(X(j)/(X(k)-X(j);endendL(k,:)=V;endC=Y*L;然后在命令窗口中输入以下内容： x=0 0.25 0.5 0.75 1;y=0 0.3104 0.6177 0.7886 1;lagran(x,y)ans = 3.3088 -6.3851 3.3164 0.7599 0得到的数据就是

6、多项式各项的系数，注意最后一个是常数项，即x0，所以表达式为:f=3.3088*x.4-6.3851*x.3+3.3164*x.2 +0.7599*x求面积就是积分求解 f=(x)3.3088*x.4-6.3851*x.3+3.3164*x.2 +0.7599*x; quad(f,0,1)ans = 0.5509 这些点肯定是通过这个多项式的！ MATLAB插值与拟合1曲线拟合实例：温度曲线问题气象部门观测到一天某些时刻的温度变化数据为：t012345678910T1315171416192624262729试描绘出温度变化曲线。曲线拟合就是计算出两组数据之间的一种函数关系，由此可描绘其变化曲

7、线及估计非采集数据对应的变量信息。曲线拟合有多种方式，下面是一元函数采用最小二乘法对给定数据进行多项式曲线拟合，最后给出拟合的多项式系数。1.线性拟合函数：regress()调用格式： b=regress(y,X) b,bint,r,rint,stats= regress(y,X) b,bint,r,rint,stats= regress(y,X,alpha)说明：b=regress(y,X)返回X与y的最小二乘拟合值，及线性模型的参数值、。该函数求解线性模型：y=X+是p1的参数向量；是服从标准正态分布的随机干扰的n1的向量；y为n1的向量；X为np矩阵。bint返回的95%的置信区间。r中

8、为形状残差，rint中返回每一个残差的95%置信区间。Stats向量包含R2统计量、回归的F值和p值。例1：设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。即y=10+x+ ；求线性拟合方程系数。程序： x=ones(10,1) (1:10);y=x*10;1+normrnd(0,0.1,10,1);b,bint=regress(y,x,0.05)结果： x = 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 1 10y = 10.9567 11.8334 13.0125 14.0288 14.8854 16.1191 17.1189 17.9962 1

9、9.0327 20.0175b = 9.9213 1.0143bint = 9.7889 10.0537 0.9930 1.0357即回归方程为：y=9.9213+1.0143x2.多项式曲线拟合函数：polyfit( )调用格式： p=polyfit(x,y,n) p,s= polyfit(x,y,n)说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。(见下一函数polyval)例2：由离散数据x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.52拟合出多项式。程序： x=0:.1:1;y=.3

10、 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2;n=3;p=polyfit(x,y,n)xi=linspace(0,1,100);%linspace用于创建向量，如：x=linspace(a1,a2,a3);a1为第一个元素，a2为最末一个元素，a3表示x共有a3个元素，每个元素间距相等。z=polyval(p,xi); %多项式求值plot(x,y,o,xi,z,k:,x,y,b)legend(原始数据,3阶曲线)结果：p = 16.7832 -25.7459 10.9802 -0.0035多项式为：16.7832x3-25.7459x2+10.9802x-0.0035曲线

11、拟合图形：如果是n=6，则如下图：也可由函数给出数据。例3：x=1:20,y=x+3*sin(x)程序： x=1:20; y=x+3*sin(x); p=polyfit(x,y,6) xi=linspace(1,20,100); z=polyval(p,xi); %多项式求值函数 plot(x,y,o,xi,z,k:,x,y,b) legend(原始数据,6阶曲线)结果：p =0.0000 -0.0021 0.0505 -0.5971 3.6472 -9.7295 11.3304再用10阶多项式拟合 程序：x=1:20;y=x+3*sin(x);p=polyfit(x,y,10)xi=lins

12、pace(1,20,100);z=polyval(p,xi);plot(x,y,o,xi,z,k:,x,y,b)legend(原始数据,10阶多项式)结果：p = Columns 1 through 7 0.0000 -0.0000 0.0004 -0.0114 0.1814 -1.8065 11.2360 Columns 8 through 11 -42.0861 88.5907 -92.8155 40.2671可用不同阶的多项式来拟合数据，但也不是阶数越高拟合的越好。3. 多项式曲线求值函数：polyval( )调用格式： y=polyval(p,x) y,DELTA=polyval(p,

13、x,s)说明：y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。y,DELTA=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则Y DELTA将至少包含50%的预测值。(未完)4. 多项式曲线拟合的评价和置信区间函数：polyconf( )调用格式： Y,DELTA=polyconf(p,x,s) Y,DELTA=polyconf(p,x,s,alpha)说明：Y,DELTA=polyconf(p,x,s)使用polyfit函数的选项输出s给出Y的95%置信区

14、间Y DELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。1-alpha为置信度。例4：给出上面例1的预测值及置信度为90%的置信区间。程序： x=0:.1:1;y=.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2n=3;p,s=polyfit(x,y,n)alpha=0.05;Y,DELTA=polyconf(p,x,s,alpha) 结果：p = 16.7832 -25.7459 10.9802 -0.0035s =R: 4x4 double df: 7normr: 1.1406Y = Columns 1 through 9 -0.0035

15、0.8538 1.2970 1.4266 1.3434 1.1480 0.9413 0.8238 0.8963 Columns 10 through 11 1.2594 2.01405. 稳健回归函数：robust( )稳健回归是指此回归方法相对于其他回归方法而言，受异常值的影响较小。调用格式： b=robustfit(x,y) b,stats=robustfit(x,y) b,stats=robustfit(x,y,wfun,tune,const)说明：b返回系数估计向量；stats返回各种参数估计；wfun指定一个加权函数；tune为调协常数；const的值为on(默认值)时添加一个常数项

16、；为off 时忽略常数项。例5：演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数y=10-2x加上一些随机干扰的项生成数据集，然后改变一个y的值形成异常值。调用不同的拟合函数，通过图形观查影响程度。程序：x=(1:10);y=10-2*x+randn(10,1);y(10)=0;bls=regress(y,ones(10,1) x) %线性拟合brob=robustfit(x,y) %稳健拟合scatter(x,y)hold onplot(x,bls(1)+bls(2)*x,:)plot(x,brob(1)+brob(2)*x,r)结果 ： bls = 8.4452 -1.47

17、84brob = 10.2934 -2.0006分析：稳健拟合(实线)对数据的拟合程度好些，忽略了异常值。最小二乘拟合(点线)则受到异常值的影响，向异常值偏移。6. 向自定义函数拟合对于给定的数据，根据经验拟合为带有待定常数的自定义函数。所用函数：nlinfit( )调用格式： beta,r,J=nlinfit(X,y,fun,beta0)说明：beta返回函数fun中的待定常数；r表示残差；J表示雅可比矩阵。X,y为数据；fun自定义函数；beta0待定常数初值。例6：在化工生产中获得的氯气的级分y随生产时间x下降，假定在x8时，y与x之间有如下形式的非线性模型：现收集了44组数据，利用该数

18、据通过拟合确定非线性模型中的待定常数。x y x y x y8 0.49 16 0.43 28 0.418 0.49 18 0.46 28 0.4010 0.48 18 0.45 30 0.4010 0.47 20 0.42 30 0.4010 0.48 20 0.42 30 0.3810 0.47 20 0.43 32 0.4112 0.46 20 0.41 32 0.4012 0.46 22 0.41 34 0.4012 0.45 22 0.40 36 0.4112 0.43 24 0.42 36 0.3614 0.45 24 0.40 38 0.4014 0.43 24 0.40 38

19、0.4014 0.43 26 0.41 40 0.3616 0.44 26 0.40 42 0.3916 0.43 26 0.41 首先定义非线性函数的m文件：fff6.mfunction yy=model(beta0,x) a=beta0(1); b=beta0(2); yy=a+(0.49-a)*exp(-b*(x-8); 程序：x=8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.00 14.00 14.00 14.00. 16.00 16.00 16.00 18.00 18.00 20.00 20.00 20.00 20.00 22.00

20、 22.00 24.00. 24.00 24.00 26.00 26.00 26.00 28.00 28.00 30.00 30.00 30.00 32.00 32.00. 34.00 36.00 36.00 38.00 38.00 40.00 42.00; y=0.49 0.49 0.48 0.47 0.48 0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43. 0.43 0.46 0.42 0.42 0.43 0.41 0.41 0.40 0.42 0.40 0.40 0.41 0.40 0.41 0.41. 0.40 0.40 0.40 0

21、.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39 0.39; beta0=0.30 0.02;betafit = nlinfit(x,y,model ,beta0)结果：betafit = 0.38960.1011 即：a=0.3896 ，b=0.1011 拟合函数为：yy=0.3896+0.1004*exp(-0.1011*(x-8)excell数据调入matlab中的方法：import datadatax=data(:,所调出数据的列数)如果想调出某列数据中间隔相等的数据则：x=data(1:2:10,colume2)则调用原始数据的第二列中前十个数据

22、中1、3、5、7、9行数据。2 插值问题在应用领域中，由有限个已知数据点，构造一个解析表达式，由此计算数据点之间的函数值，称之为插值。实例：海底探测问题某公司用声纳对海底进行测试，在55海里的坐标点上测得海底深度的值，希望通过这些有限的数据了解更多处的海底情况。并绘出较细致的海底曲面图。一、一元插值一元插值是对一元数据点(xi,yi)进行插值。1 线性插值：由已知数据点连成一条折线，认为相临两个数据点之间的函数值就在这两点之间的连线上。一般来说，数据点数越多，线性插值就越精确。调用格式：yi=interp1(x,y,xi,linear) %线性插值zi=interp1(x,y,xi,splin

23、e) %三次样条插值wi=interp1(x,y,xi,cubic) %三次多项式插值说明：yi、zi、wi为对应xi的不同类型的插值。x、y为已知数据点。例1：已知数据：x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.52求当xi=0.25时的yi的值。程序：x=0:.1:1;y=.3 .5 1 1.4 1.6 1 .6 .4 .8 1.5 2;yi0=interp1(x,y,0.025,linear)xi=0:.02:1;yi=interp1(x,y,xi,linear);zi=interp1(x,y,xi,spline);wi=interp1(x,

24、y,xi,cubic);plot(x,y,o,xi,yi,r+,xi,zi,g*,xi,wi,k.-)legend(原始点,线性点,三次样条,三次多项式)结果：yi0 = 0.3500要得到给定的几个点的对应函数值，可用：xi = 0.2500 0.3500 0.4500yi=interp1(x,y,xi,spline)结果：yi =1.2088 1.5802 1.3454(二) 二元插值二元插值与一元插值的基本思想一致，对原始数据点(x,y,z)构造见世面函数求出插值点数据(xi,yi,zi)。一、单调节点插值函数，即x,y向量是单调的。调用格式1：zi=interp2(x,y,z,xi,yi,linear)liner