matlab在数值分析中的应用.pptx

《matlab在数值分析中的应用.pptx》由会员分享，可在线阅读，更多相关《matlab在数值分析中的应用.pptx（41页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、17.1 多项式插值原理7.2 分段线性插值7.3 离散数据的多项式拟合7.3.1 线性最小二乘法拟合7.3.2 非线性最小二乘法拟合第1页/共41页2 插值与拟合是来源于实际、又广泛应用于实际的两种重要方法。随着计算机的不断发展及计算水平的不断提高，它们已在国民生产和科学研究等方面扮演着越来越重要的角色。第2页/共41页37.1 多项式插值原理多项式插值原理已知 f(x)在点 xi 上的函数值 yi=f(xi),(i=0,1,2,n)求多项式 P(x)=a0+a1x+anxn满足:P(xk)=yk (k=0,1,n)P10(t)f(t)f(x)第3页/共41页47.2 分段线性插值分段线性插

2、值 所谓分段线性插值就是通过插值点用折线段连接起来逼近原曲线，这也是计算机绘制图形的基本原理。实现分段线性插值不需编制函数程序，MATLAB自身提供了内部函数interp1其主要用法如下：interp1(x,y,xi)一维插值第4页/共41页5 yi=interp1(x,y,xi)对一组点(x,y)进行插值，计算插值点xi的函数值。x为节点向量值，y为对应的节点函数值。如果y 为矩阵，则插值对y 的每一列进行，若y 的维数超出x 或 xi 的维数，则返回NaN。yi=interp1(y,xi)此格式默认x=1:n，n为向量y的元素个数值，或等于矩阵y的size(y,1)。第5页/共41页6 y

3、i=interp1(x,y,xi,method)method用来指定插值的算法。默认为线性算法。其值常用的可以是如下的字符串。nearest 线性最近项插值。linear 线性插值。spline 三次样条插值。cubic 三次插值。pchip 分段三次Hermite插值。所有的插值方法要求x是单调的。x 也可能并非连续等距的。第6页/共41页7正弦曲线的插值示例：x=0:10;y=x.*sin(x);xx=0:.25:10;yy=interp1(x,y,xx);plot(x,y,kd,xx,yy)第7页/共41页8例 产品产量拟合（离散数据的一维插值）year=1900:10:2010;pro

4、duct=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893;p1995=interp1(year,product,1995)x=1900:1:2010;y=interp1(year,product,x,pchip);plot(year,product,o,x,y)第8页/共41页9插值结果为：p1995=252.9885第9页/共41页10 Matlab也能够完成二维插值的运算，相应的函数为interp2，使用方法与interpl基本相同，只是输入和输出的参

5、数为矩阵，对应于二维平面上的数据点。二维插值zi=interp2(x,y,z,xi,yi,method)第10页/共41页11 返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)Xi(i,j),yi(i,j)。用户可以输入行向量和列向量Xi与Yi，此时，输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回NAN（Not a Number）。第11页/共41

6、页12二维插值图例：X,Y=meshgrid(-3:.25:3);Z=peaks(X,Y);XI,YI=meshgrid(-3:.125:3);ZZ=interp2(X,Y,Z,XI,YI);surfl(X,Y,Z);hold on;surfl(XI,YI,ZZ+15)axis(-3 3-3 3-5 20);shading flathold off第12页/共41页13例years=1950:10:1990;service=10:10:30;wage=150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767

7、226.505 153.706 426.730 249.633 120.281 598.243;w=interp2(service,years,wage,15,1975)插值结果为：w=190.6288第13页/共41页14三维插值vi=interp3(x,y,z,v,xi,yi,zi,method)找出由参量X,Y,Z决定的三元函数V=V(X,Y,Z)在点（XI,YI,ZI）的值。参量XI,YI,ZI是同型阵列或向量。若向量参量XI,YI,ZI是不同长度，不同方向（行或列）的向量，这时输出参量VI与Y1,Y2,Y3为同型矩阵。其中Y1,Y2,Y3为用命令meshgrid(XI,YI,ZI)生

8、成的同型阵列。若插值点(XI,YI,ZI)中有位于点(X,Y,Z)之外的点，则相应地返回特殊变量值NaN。第14页/共41页15例x,y,z,v=flow(20);xx,yy,zz=meshgrid(.1:.25:10,-3:.25:3,-3:.25:3);vv=interp3(x,y,z,v,xx,yy,zz);slice(xx,yy,zz,vv,6 9.5,1 2,-2.2);shading interp;colormap cool 第15页/共41页167.3 离散数据的多项式拟合 在科学实验的统计方法研究中,往往要从一组实验数据(xi,yi)中寻找出自变量x 和因变量y之间的函数关系y

9、=f(x)。第16页/共41页17第17页/共41页18第18页/共41页19已知数据表 x x1 x2 xmf(x)y1 y2 ym求拟合函数:(x)=a0+a1x+anxn使得达到最小通俗的讲：第19页/共41页207.3.1 线性最小二乘法拟合第20页/共41页21P=polyfit(X,Y,n)求出(最小二乘意义下)n次拟合多项式P(x)=a0 xn+a1xn-1+an-1x+an计算结果为系数:P=a0,a1,an-1,an 在MATLAB中实现最小二乘法拟合通常采用polyfit函数进行。第21页/共41页22 函数polyfit是指用一个多项式函数来对已知数据进行拟合：x=0:0

10、.1:1;y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2 P=polyfit(x,y,2)P=-9.8108 20.1293 -0.0317函数返回的是一个多项式系数的行向量，写成多项式形式为：第22页/共41页23 为了使用polyfit，首先必须指定我们希望以多少阶多项式对以上数据进行拟合，如果我们指定一阶多项式，结果为线性近似，通常称为线性回归。我们选择二阶多项式进行拟合。第23页/共41页24多项式求值命令：y1=polyval(P,x)其中，P是n次多项式的系数，x是自变量的值，y1是多项式在x处的值。第24页/

11、共41页25 为了比较拟合结果，我们绘制两者的图形：xi=linspace(0,1,100);%绘图的X-轴数据。Z=polyval(p,xi);%得到多项式在数据点处的值。当然，我们也可以选择更高幂次的多项式进行拟合，如10阶：p=polyfit(x,y,10);xi=linspace(0,1,100);z=ployval(p,xi);第25页/共41页26 可以上机绘图进行比较，曲线在数据点附近更加接近数据点的测量值了，但从整体上来说，曲线波动比较大，并不一定适合实际使用的需要，所以在进行高阶曲线拟合时，“越高越好”的观点不一定对的。第26页/共41页27例 血药浓度问题 某人用快速静脉注

12、射方式一次注入药物300mg，在一定时刻采集血样，测得血药浓度如下：t=0.25,0.5,1,1.5,2,3,4,6,8c=19.21,18.15,15.36,14.10,12.89,9.32,7.45,5.24,3.01拟合血药浓度随时间的变化规律。第27页/共41页28第28页/共41页29用如下Matlab命令（xueyao.m）t=0.25,0.5,1,1.5,2,3,4,6,8;c=19.21,18.15,15.36,14.10,12.89,9.32,7.45,5.24,3.01;y=log(c);aa=polyfit(t,y,1);a=aa(1)b=aa(2)k=-ad=300;v

13、=d/exp(b)cc=exp(b)*exp(a*t);plot(t,c,k+,t,cc,r)第29页/共41页30a=-0.2347b=2.9943k=0.2347v=15.0219第30页/共41页317.3.2 非线性最小二乘法拟合第31页/共41页32Matlab命令：（1）Lsqnonlin(f,a0)其中，f.m是描述函数R(a)的函数文件名，a0是初值。第32页/共41页33例 血药浓度问题function f=ct(x)t=0.25,0.5,1,1.5,2,3,4,6,8;c=19.21,18.15,15.36,14.10,12.89,9.32,7.45,5.24,3.01;f

14、=c-x(1)*exp(x(2)*t);nongdu.mx0=10,0.5lsqnonlin(ct,x0)ans=20.2413 -0.2420第33页/共41页34（2）lsqcurvefit(fun,a0,x,y,LB,UB)其中，f.m是描述函数f(x,a)的函数文件名，a0是初值，x=x1,xn，y=y1,yn为数据第34页/共41页35例 血药浓度问题function f=xue(x,t)f=x(1)*exp(x(2)*t);xuenon.mt=0.25,0.5,1,1.5,2,3,4,6,8;c=19.21,18.15,15.36,14.10,12.89,9.32,7.45,5.2

15、4,3.01;x0=10,0.5x=lsqcurvefit(xue,x0,t,c);yy=xue(x,t);plot(t,c,k*,t,yy,r)第35页/共41页36第36页/共41页37例 药物的吸收与排除（口服）假设口服剂量为d，中心室血药浓度数据为t=0.083,0.167,0.25,0.50,0.75,1.0,1.5,2.25,3.0,4.0,6.0,8.0,10.0,12.0c=10.9,21.1,27.3,36.4,35.5,38.4,34.8,24.2,23.6,15.7,8.2,8.3,2.2,1.8 由机理分析（吸收室-中心室模型）得，中心室的血药浓度其中排除速率和吸收速率

16、k1由上述数据拟合。第37页/共41页38function f=xipai(x,t)f=x(1)*x(3)*(exp(-x(2)*t)-exp(-x(1)*t)/(x(1)-x(2);xipainon.mt=0.083,0.167,0.25,0.50,0.75,1.0,1.5,2.25,3.0,4.0,6.0,8.0,10.0,12.0;c=10.9,21.1,27.3,36.4,35.5,38.4,34.8,24.2,23.6,15.7,8.2,8.3,2.2,1.8;x0=1,0,40;x=lsqcurvefit(xipai,x0,t,c);yy=xipai(x,t);plot(t,c,k*,t,yy,r)第38页/共41页39第39页/共41页40第40页/共41页41谢谢您的观看！第41页/共41页