《MATLAB数值计算》PPT课件.ppt

第6章 MATLAB数值计算,Page 2,6.1.1 数据统计与分析 1. 求矩阵最大元素和最小元素 （1）求向量的最大值和最小值 y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。 y,I=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。 求向量X的最小值的函数是min(X)，用法和max(X)完全相同。,6.1 数据处理与多项式计算,Page 3,例 求向量x的最大值。 命令如下： x=-43,72,9,16,23,47; y=max(x) %求向量x中的最大值 y,l=max(x) %求向量x中的最大值及其该元素的位置,多输入函数 例如： x=max(z , y) 多输出函数多个输出值用 括起来，且输出值之间用逗号隔开。,6.1 数据处理与多项式计算,Page 4,（2）求矩阵的最大值和最小值 max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值。 Y,U=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。,max(A,dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。,6.1 数据处理与多项式计算,Page 5,6.1 数据处理与多项式计算, A=13,-56,78;25,63,-235;78,25,563;1,0,-1; a1 = max(A,2); a2 = min(A,2); a3 = max(A); a4 = min(A); a5 = max(max(A); a6 = min(A(:);,Page 6,（3）两个向量或矩阵对应元素的比较 函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为： U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。 U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。 例 求两个23矩阵x, y所有同一位置上的较大元素构成的新矩阵p。,6.1 数据处理与多项式计算, x=4,5,6;1,4,8 x = 4 5 6 1 4 8 y=1,7,5;4,5,7 y = 1 7 5 4 5 7 p=max(x,y) p = 4 7 6 4 5 8 p2=max(y,5) p2 = 5 7 5 5 5 7,Page 7,2. 求矩阵的平均值和中值 3. 矩阵元素求和与求积 7. 排序,6.1 数据处理与多项式计算,Page 8,6.1.2 数据插值 一维数据插值,插值：是在认定所给“基准数据”完全正确的情况下，研 究如何“平滑”的估算出“基准数据”之间其它点的函数值。因此，插值所得曲线一定穿过“基准数据”。,6.1 数据处理与多项式计算,Page 9,函数根据x,y的值，计算函数在xs处的值。x,y是两个等长的已知向量，分别描述采样点和样本值，xs是一个向量或标量，描述欲插值的点，ys是一个与xs等长的插值结果。 其中:(1) x,y是量测数据对； (2) xs是需要内插的点所构成的向量。 (3) method是指所使用的内插方法。,ys=interp1(x,y,xs,method);,说明：interp1仅是插值指令的一种，还有interp2 、interp3等。,插值算法：nearest，linear，spline，cubic,注意：xs的取值范围不能超出x的给定范围，否则，会给出“NaN”错误。,6.1 数据处理与多项式计算,Page 10,%给定数据对 x0=0:0.1:1; y0=-.44,1.97,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22; %采用三次多项式进行插值 xi=0:0.02:1; yi=interp1(x0,y0,xi,cubic); %绘图 plot(xi,yi,-b,x0,y0,.r,MarkerSize,20),xlabel(x),6.1 数据处理与多项式计算,Page 11,例:某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度()，用3次样条插值分别求得该日室内外6:00至18:00时之间每隔0.5小时各点的近似温度()。 设时间变量h为一行向量，温度变量t为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如下： h =6:2:18; t=18,20,22,25,30,28,24;15,19,24,28,34,32,30; XI =6:0.5:18; YI=interp1(h,t,XI, spline); %用3次样条插值计算 plot(XI,YI, -,h,t, .,MarkerSize,20),xlabel(时间),6.1 数据处理与多项式计算,Page 12,1）最邻近插值方法（nearest）,插值点的值与其最邻近的点的函数值相等。,*,*,*,*,*,6.1 数据处理与多项式计算,Page 13,2）线性插值方法（ linear ）,插值点的值在前，后两个数据点所构成的直线上。,*,*,*,*,*,6.1 数据处理与多项式计算,Page 14,3）三次样条插值方法(spline),利用一系列样条函数获得内插数据点，从而确定已有数据点之间的函数。,4）三次曲线内插方法（cubic）,构造三次曲线函数来拟合已知数据x、y，从而确定内插点的值。,6.1 数据处理与多项式计算,Page 15,说明: 1、四种插值方法中，x中的数据是单调但不一定距均匀的。 2、若已知x为均匀的，则在method前加*，可使执行速度加快。 3、按nearest linear cubic spline的顺序，对内存要求 从小到大，执行速度由快到慢，平滑度由差到好。,6.1 数据处理与多项式计算,Page 16,6.1.3 曲线拟合,p=polyfit(x,y,n);,多项式拟合,6.1 数据处理与多项式计算,Page 17,6.1 数据处理与多项式计算,Page 18,在MATLAB中，用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。,p=polyfit(x,y,n);,函数根据采样点x和采样点函数值y，产生一个n次多项式p。其中x,y是两个等长的向量，p是一个长度为n+1的向量，p的元素为多项式系数。 polyval函数的功能是按多项式的系数计算x点多项式的值.,6.1 数据处理与多项式计算,Page 19,6.1 数据处理与多项式计算,%给定数据对 x0=0:0.1:1; y0=-.44,1.97,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22; %求拟合多项式系数 n=3; P=polyfit(x0,y0,n) %图示拟合情况 xx=0:0.01:1; yy=polyval(P,xx); plot(xx,yy,-b,x0,y0,.r,MarkerSize,20),xlabel(x),拟合：寻找一条“平滑”曲线来最好的表现带噪声的“测量数据”。 但并不要求拟合曲线穿过这些“测量数据”点。,Page 20,例6.11 用一个3次多项式在区间0，2内逼近函数。 命令如下： X=linspace(0,2*pi,50); Y=sin(X); P=polyfit(X,Y,3) %得到3次多项式的系数和误差 Y1 = polyval(P,X); plot(X,Y,:o,X,Y1,-*); legend(原始曲线,拟合曲线);,6.1 数据处理与多项式计算,Page 21,6.1.4 多项式计算,1、多项式的表示： MATALB中，用一个向量来表示多项式。这个 向量中按照降幂的顺序排列多项式的各项系数。,在MATLAB中为：p=1,0,-2,5;,说明：如果多项式中缺某幂次项，则认为该项的系数为零。,6.1 数据处理与多项式计算,Page 22,1. 多项式的四则运算 （1）多项式的加减运算 （2）多项式乘法运算 因为多项式是用其系数构成的离散序列表示的，因此多项式的乘法和除法就对应于卷积和去卷积的操作。,6.1 数据处理与多项式计算, a=1,-2,5,3; b=0,0,6,-1; c=a+b c = 1 -2 11 2,Page 23,1）多项式乘法卷积：,c=conv(a,b);,例如: (x-1)(x-2)=x2-3x+2,求多项式a和b的乘积/卷积。a、b是两个多项式系数向量。, a=1,-1; b=1,-2; c=conv(a,b) c = 1 -3 2,6.1 数据处理与多项式计算,Page 24,例如：,q,r=deconv(c,b);, c=1,-3,5; b=1,-2; q,r=deconv(c,b) q = 1 -1 r = 0 0 3 conv(q,b)+r ans = 1 -3 5,其中q返回多项式c除以b的商式，r返回余式。q和r仍是多项式系数向量。,deconv是conv的逆函数， 即有c=conv(b,q)+r。,（3）多项式除法解卷积,6.1 数据处理与多项式计算,Page 25,2. 多项式的导函数,q=polyder(p);,具体用法: (1)q=polyder(p)可以得到多项式p的导数。,(2) q=polyder(a,b)可以求出多项式a,b之积的导数。,(3)q,d=polyder(b,a)可以求出多项式之比b(s)/a(s)的导数。,6.1 数据处理与多项式计算, c=1,-3,5; polyder(c) ans = 2 -3,Page 26,3. 多项式求值,y=polyval(p,x);,y=polyvalm(p,x);,功能：按数组运算规则计算多项式的值。 其中x可以是标量和数组。,功能：按矩阵运算规则计算多项式的值。 其中x必须为方阵。,6.1 数据处理与多项式计算,Page 27,P= 2 3 1 ; X= 1 2 3 ; 3 2 1 ; 2 1 3 ;,y=polyval ( P , X ),y=polyvalm ( P , X ),y=2*X.*X+3*X+1,y=2*X*X+3*X+eye(3),6.1 数据处理与多项式计算,Page 28,4. 多项式求根,R=roots(P), p=1,-3,2; R = roots(p) R = 2 1 p1=poly(R),功能：计算多 项式P的根。,6.1 数据处理与多项式计算,若已知多项式的全部根，则可以用poly函数建立起该多项式：P=poly(x) 当x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋给向量P。,Page 29,6.2 数值微积分,6.2.1 数值微分 1. 数值差分与差商 2. 数值微分的实现 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为： DX=diff(X)：计算向量X的向前差分 DX=diff(X,n)：计算X的n阶向前差分。,Page 30,clear;clc; p = 1 2 1; x=-1:0.01:1; dp = polyder(p); dpx = polyval(dp,x); y = x.2 + 2*x +1; dx = diff(y)/0.01; x1 = x(1:end-1); gx = 2*x+2; plot(x,dpx,r,x1,dx,g,x,gx,b); figure(2) subplot(1,3,1); plot(x,dpx,r); subplot(1,3,2); plot(x1,dx,g); subplot(1,3,3); plot(x,gx,b);,例：f(x) = x2+2x+1 用不同方法求其倒数，并画出图象,Page 31,(1)被积函数是一个解析式,例：求以下的定积分。其精确值为0.7468204,%（1）使用字符串表示被处理函数 fun=exp(-x.*x);,%（2）调用积分指令求积分值 E1 = quad (fun,0,1) E2 = quadl (fun,0,1),区别：求数值积分是采取不同的方法,6.2 数值微积分,6.2.2 数值积分2. 数值积分的实现,Page 32,6.4.1 直接解法 1利用左除运算符的直接解法 对于线性方程组Ax=b，可以利用左除运算符“”求解： x=Ab,6.4 线性方程组求解,Page 33,6.4.3 求线性方程组的通解 线性方程组的求解分为两类：一类是求方程组的惟一解即特解，另一类是求方程组的无穷解即通解。这里对线性方程组 Ax=b的求解理论作一个归纳：p172,6.4 线性方程组求解,Page 34,6.5.1 非线性方程数值求解 1. 单变量非线性方程求解 z=fzero(fname,x0,tol,trace) 求单变量非线性方程的根 其中fname是待求根的函数文件名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。,6.5 非线性方程与最优化问题求解,Page 35,例6.33 求f(x)=x-1/x+5在x0=-5和x0=1作为迭代初值时的零点。 先建立函数文件fz.m： function f=fz(x) f=x-1/x+5; 然后调用fzero函数求根。： fzero(fz,-5) %以-5作为迭代初值 fzero(fz,1) %以1作为迭代初值,x=-6:0.1:2; f=x-1./x+5; plot(x,f,b-,x,f,r.) hold on plot(x,zeros(1,length(x),k-),6.5 非线性方程与最优化问题求解,Page 36,常微分方程：,6.6 常微分方程的数值求解,6.6.2 龙格库塔法的实现 基于龙格库塔法的求常微分方程数值解的函数： t,y