CHAPT4 数值计算

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,Chapt4 数值计算,4.1 插值与数据拟合,4.2 多项式,4.3 微分与积分,4.4 线性方程与非线性方程,4.5 常微分方程,1,4.1 插值与数据拟合,在分析试验数据中，常常要面临将试验数据作解析描述的任务，这个问题有插值和曲线拟合两种方法,。,当数据量不够，需要补充，且认定已有数据可信时, 通常利用函数插,值，插值是一个非常有价值的工具。,在曲线拟合中，假定已知曲线的规律，作曲线的最佳逼近，但不需要经过所有的数据点,；,4.1.1,插值函数,多项式插值原理,已知函数,f(x),在区间a,b上有定义且已知，则可以求出点,x,i,（a,x,i,b）,上的函数值,y,i,=,f,(,x,i,), (,i,=0,1,2,n)。,如果,x,i,出界，则对应值为,NaN,2,1、曲线插值函数,interp1,一维数值插值,t=interp1(x,y,x,0,method),x、y：,原始数据点,，,是两个等长的已知向量。,x,0,为进行插值的数组,，,计算函数在,x,0,处的值。,method,为插值算法,:,该缺省的使用假定为线性插值,3,当x,0,为等距时可以用快速插值法，使用快速插值法的格式为*nearest、*linear、*spline、*cubic,以上五种方法得出的数据值一个比一个精确,而所需内存及计算时间也一个比一个要大要长.,nearest,方法,:,将插值点的数据设为与之最接近的邻点的数值,.,执行速度最快，输出结果为直角转折,linear,方法,:,利用每两个相邻插值节点作线性插值,.,默认值，在样本点上斜率变化很大,样条插值方法：分段多项式的在连接处光滑连接。最花时间，但输出结果也最平滑,Hermite,方法：对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值。,Cubic,：三次,插值，要求,x,0,必须为等距。最占内存，输出结果与,spline,差不多,4,例1：已知1900年到2010年每隔十年的数据如下：,75.995 91.972 105.711 123.203 131.669 150.697,179.323 203.212 226.505 249.633 256.344 267.893,用插值法求1995年的数据。,命令如下：,T =1900:10:2010;,Y=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893;,Y,1,=interp1(T,Y,1995),设时间变量T为一行向量，数据用Y为一个行,向量,x,0,=1995.,Y,1,=,252.9885,plot(T,Y,1995,Y,1,Or),t=interp1(x,y,x,0,method),5,例2：已知1900年到2010年每隔十年的数据如下：,75.995 91.972 105.711 123.203 131.669 150.697,179.323 203.212 226.505 249.633 256.344 267.893,用三种不同的插值方法进行比较。,T =1900:10:2010;,Y=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893;,T1=1900:5:2010,Y1=interp1(T,Y,T1,nearest);,Y2=interp1(T,Y,T1,linear);,Y3=interp1(T,Y,T1,spline);,plot(T,Y,-b,T1,Y1,-r,T1,Y2,xk,T1,Y3,-y);,legend(原曲线,最近点插值,线性插值,样条插值)；,6,注意：,1.首先x,0,的取值范围不能超出x的给定范围，否则，会给出“NaN”错误。,2.其次，独立变量必须是单调的。即独立变量在值上必须总是增加的或总是减小的。,3.最后，不同的插值方法得到的结果是不相同的，所以，在确定方法之前，要考虑该插值的应用领域以及插值的曲线形状。,7,例3 某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度()，用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度()。,设时间变量h为一行向量，温度变量t为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。,h =6:2:18;,t=18,20,22,25,30,28,24;15,19,24,28,34,32,30;,X1=6.5:2:17.5,Y1=interp1(h,t,XI,spline),plot(h,t,-b,X1,Y1,-+r);,legend(原曲线,样条插值)；,YI =,18.5020 15.6553,20.4986 20.3355,22.5193 24.9089,26.3775 29.6383,30.2051 34.2568,26.8178 30.9594,8,2,、曲面插值函数,interp2,二维数值插值,z,i,=interp2(x,y,z,x,i,y,i,method),其中,X,，,Y,是两个向量，分别描述两个参数的采样点，,Z,是与参数采样点对应的函数值。,x,i,，,y,i,是两个向量或标量，描述欲插值的点，,Z,i,是根据相应的插值方法得到的插值结果。,method,的取值与一维插值函数相同。,X,Y,Z,也可以是矩阵形式。,9,clc;clear;,x,y=meshgrid(-3:3); %创建x和y矩阵,z=peaks(x,y);,surf(x,y,z);,xi,yi=meshgrid(-3:0.5:3);对数据插值加密,z1=interp2(x,y,z,xi,yi,nearest);,z2=interp2(x,y,z,xi,yi,linear);,z3=interp2(x,y,z,xi,yi,bicubic);,figure,subplot(121); surf(xi,yi,z1);,subplot(122);contour(xi,yi,z1);等值线,title(nearest);,figure,subplot(121); surf(xi,yi,z2);,subplot(122);contour(xi,yi,z2);,title(linear);,figure,subplot(121); surf(xi,yi,z3);,subplot(122);contour(xi,yi,z3);,title(bicubic);,);,10,11,4.1.2,数据拟合,所谓数据拟合是求一个简单的函数,例如是一个低次多项式,不要求通过已知的这些点,而是要求在整体上“尽量好”的逼近原函数。这时,在每个已知点上就会有误差,数据拟合就是从整体上使误差,尽量的小一些。,已知数据表,x,x,1,x,2,x,m,f,(,x,),y,1,y,2,y,m,最小,12,p=polyfit(x,y,n),p,s=polyfit(x,y,n),p,s,mu=polyfit(x,y,n),多项式拟合命令,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;,p=polyfit(x, y, 2),其中,X,，,Y,是两个向量，是已知的数据点。,n,是阶次，如果我们选择,n=1,，,得到最简单的线性近似。通常称为线性回归。相反，如果我们选择,n=2,，,得到一个,2,阶多项式。,p =,-9.8108 20.1293 -0.0317,polyfit 的输出是一个多项式系数的行向量。,其解是y = 9.8108x,2,20.1293x0.0317,多项式阶次的选择是有点任意的,，但是,选择的过,小，拟合效果不好；阶次选择的过高，虽然数据点上看到效果好，数据点之间会,出现有数据振荡的问题，阶数不宜过高，小于5阶。,13,xi=linspace(0, 1, 100);,z=polyval(p, xi);函数值,plot(x,y,o,x,y,xi,z,:),为了将曲线拟合解与数据点比较，让我们把二者都绘成图。,14,4.2 多项式,1,、多项式表示,多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。 （,表示中需要包含零系数的项,）,f(x)=a,n,x,n,+a,n-1,x,n-1,+a,0,可用行向量 p=a,n,a,n-1, a,0,表示,例如:x,4,12x,3,25x11,p=1 -12 0 25 11,15,规定,：,多项式用行向量，根用列向量。,给出多项式的根，使用,poly,函数也可以构造出相应的多项式,2、roots 求多项式的根,p=1 -12 0 25 11,r=roots(p),r =,11.8142,1.7471,-1.0573,-0.5041,16,1给出一个多项式的根（向量），也可以构造相应的多项式。,3、poly 产生特征多项式系数向量,p=1 -12 0 25 11,r=roots(p),p2=poly(r),p2=,1.0000 -12.0000 -0.0000 25.0000 11.0000,注意：,特征多项式一定是n+1维的。,特征多项式第一个元素一定是1。,17,a=1 2 3;4 5 6;7 8 0;,p=poly(a),2给出一个矩阵，返回的是特征多项式的系数,p =,1.00 -6.00 -72.00 -27.00p是矩阵a的特征多项式的系数,r=roots(p),r = 12.12,-5.73 r是矩阵a的特征值,-0.39,p2=poly2sym(p,x),表达成多项式形式,p2=,x3-6*x2-72*x-27,p1=poly2str(p,x),将多项式系数向量表示成符号多项式,p1 =,x3 - 6 x2 - 72 x - 27,18,4,、多项式运算,对多项式加法，MATLAB不提供一个直接的函数。如果两个多项式向量大小相同，标准的数组加法有效。,多项式运算,命令,乘法,（,卷积）,conv,，,convs,除法,（,反卷积）,deconv,微分,polyder(p),求,值,polyval,d=a+b,19,1conv多项式乘法运算,a=1 2 3;b=4 5 6;,c=conv(a,b）,p=poly2str(c,y),例:a(x)=x2+2x+3; b(x)=4x2+5x+6;,c = (x2+2x+3)(4x2+5x+6),c =,4 13 28 27 18,p =,4 y4 + 13 y3 + 28 y2 + 27 y + 18,20,2deconv多项式除法运算,a=1 2 3;,c = 4 10 28 26 18,b,d=deconv(c,a),Q,r=deconv(P1,P2),返回多项式P1除以P2的商式,r返回P1除以P2的余式,c=conv(a,b)+d,b =,4 2 12,d =,0 0 0 -4 -18,c =,4 10 28 26 18,21,3多项式微分,命令格式：,polyder(p): 求p的微分,polyder(a,b): 求多项式a,b乘积的微分,p,q=polyder(a,b): 求多项式a,b商的微分,a=1 2 3 4 5;,b=polyder(a),b = 4 6 6 4,22,4多项式求值,可以求出多项式在某一点的值。,MATLAB,提供了两种求多项式值的函数：,polyval,与,polyvalm,，,它们的输入参数均为多项式系数向量,P,和自变量,x,。,两者区别：前者是代数多项式求值，而后者是矩阵多项式求值。,23,已知多项式x,4,+8x,3,-10，分别取x=1.2和一个23矩阵为自变量计算该多项式的值。,a=1 8 0 0 -10;,y1=polyval(a,1.2),x=2 3;4 5;,y2=polyval(a,x),y1 =,5.8976,y2 =,70 287,758 1615,1代数多项式求值,polyval函数用来求代数多项式的值，其调用格式为：,若,x,为一数值，则求多项式在该点的值；,若,x,为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。,Y=polyval(P,x),24,2矩阵多项式求值,polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。,polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。,Y=polyval(P,x),a=1 2 3;4 5 6;7 8 0;,p=poly(a)；,x=pascal(3);创建一个正定矩阵,polyvalm(p,x)；,ans =,-98.0000 -63.0000 -51.0000,-63.0000 -146.0000 -168.0000,-51.0000 -168.0000 -374.0000,25,第三次上机作业,1.复习本次的课件，熟悉本次课程所提到的全部内容。,2.根据所学的内容，与本专业相联系，动手操作书上的例题。,3.给出下表数据，求二次和三次拟合多项式,。,x,-1,-0.75,-0.5,-0.25,0,25,0.5,0.75,1,y,-0.2209,0.3295,0.8826,1.4392,2.0003,2.5645,3.1334,3.7061,4.2836,4、以函数 , 为例，用 “nearest,linear,spline”三插值方法进行插值，并对插值结果绘图进行比较。,26,5.给出概论积分 的数据表,x,0.46,0.47,0.48,0.49,y,0.484656,0.493745,0.50275,0.511668,求（1）x=0.472时，该积分值等于多少？,（2）当x为何值时积分值等于0.5？,6.（1）求方程,x,3, 4,x +,5 = 0 的解,（2）求,x,3,8,x,2,+6,x,30=0的解,（3）求 的“商”及“余”多项式,。,27,