MATLAB数值计算3.ppt

第三章 MATLAB数值计算,3.1矩阵运算 3.2稀疏矩阵 3.3 数据处理与多项式运算,3.1矩阵运算,3.1.1矩阵变换 3.1.2矩阵的逆与秩 3.1.3矩阵的特征值与特征向量,3.1.1矩阵变换,MATLAB提供了几个矩阵变换函数及调用方式: B=flipud(A);A对于水平轴沿上下方向进行列维翻转 B=fliplr(A); A对于垂直轴沿左右方向进行列维翻转 B=flipdim(A,dim); A沿着特定维进行翻转. B=rot90(A);矩阵A逆时针方向旋转90度 B=rot90(A,k);矩阵A逆时针方向旋转k*90度(k为整数),例,A=1,2;3,4;5,6; B=fliplr(A); B = 2 1 4 3 6 5,B=flipud(A); B = 5 6 3 4 1 2,B=rot90(A); B = 2 4 6 1 3 5,B=rot90(A,2); B = 6 5 4 3 2 1,3.1.2矩阵的逆与秩,3x1+ x2 - x3 = 3.6 x1+2x2+4x3 = 2.1 -x1+4x2+5x3 = -1.4 A=3 1 -1;1 2 4;-1 4 5;b=3.6;2.1;-1.4; Ax=b 解为:x=A-1b x=inv(A)*b x = 1.4818 -0.4606 0.3848,矩阵的秩,函数:rank() c=rank(A); c= 3,3.1.3矩阵的特征值与特征向量,矩阵的特征值 函数:eig 格式: E= eig(A) 求A的全部特征值构成向量E V,D= eig(A) 先对A作相似变换,然后求矩阵A的全部特征值,构成对角阵D,并求它的特征向量以构成V的列向量 V,D= eig(A,nobalance) 直接求A的特征向量和特征值.,例:,A=1 2 3;4 5 6;7 8 9 E=eig(A) E = 16.1168 -1.1168 -0.0000,V,D= eig(A,nobalance) V = -0.2833 -1.0000 0.5000 -0.6417 -0.1104 -1.0000 -1.0000 0.7792 0.5000 D = 16.1168 0 0 0 -1.1168 0 0 0 -0.0000,求根函数,例： a=4,9,0,-5,3,-20; x=roots(a); x = -1.8403 + 0.6250i -1.8403 - 0.6250i 1.1422 0.1442 + 1.0668i 0.1442 - 1.0668i,3.2稀疏矩阵,定义：矩阵中含有大量的零元素 存储方式：全元素存储和稀疏存储 全元素存储：每个元素需要相同的存储空间 稀疏存储：MATLAB仅存储那些非零元素及下标,对于有较多零元素的大矩阵这种存储方式可以减少存储空间.,3.2.2创建稀疏矩阵,全元素存储与稀疏存储的转换 例: A=0 0 1;0 0 0;0 0 2; s=sparse(A); f=full(s); s = f= 0 0 1 (1,3) 1 0 0 0 (3,3) 2 0 0 2,直接创建稀疏矩阵,格式: s=sparse(i,j,s,m,n); 说明: i和j分别是矩阵非零元素的行和列的下标向量,s是对应行列处的非零元素值,m和n分别是矩阵的行数和列数 例: s=sparse(1,2,3,2,3,2,5,3,6,3,3); s = (1,2) 5 (3,2) 6 (2,3) 3,3.2.3稀疏矩阵的查看与运算,函数nnz: 返回稀疏矩阵非零元素个数 例: nnz(s) ans = 3 函数nonzeros: 返回稀疏矩阵所有非零元素值 nonzeros(s) ans = 5 6 3,find函数,作用: 返回非零元素的下标和数值 格式:i,j,k=find(s) i = j= k= 1 2 5 3 2 6 2 3 3,spy函数,作用: 将矩阵中的每个非零元素用一个小点表示从点的分布可看出矩阵非零元素结构 spy(s),3.3数据处理与多项式,3.3.1数据统计与分析 求向量的最大最小元素 格式: y=max(X); y,i=max(X); i最大元素的序号 如果包含复数元素则按模取最大值. 求最小值函数为min格式与上相同,求矩阵的最大最小值,格式: y=max(A);返回一个向量,第i个元素为第 i列的最大元素 y,u=max(A);返回两个向量,y记录每列的最大值,u为每列最大元素的行号 max(A,dim);dim取1或2,取1与max相 同,取2返回一个列向量,第i个元素是A 矩阵的第i行的最大元素 两个向量或矩阵元素的比较 u=max(A,B) ;A,B为两个同型矩阵或向量 u=max(A,n);n为标量,例,A=3,2, 11,1;7, 15,4,3; 5,10, 9,18; 求每行每列的最大元素及矩阵的最大元素 y1=max(A,2); y2=max(A); y3=max(y1);,运行结果,y1 = 11 15 18 y2 = 7 15 11 18 y3 = 18,求矩阵的平均值和中值,mean(X);返回向量X的算术平均值 median(X);返回向量X的中值 mean(A); median(A); mean(A,dim); median(A,dim);,A=3 5 6;10 15 20;1 4 6; b=mean(A) c=mean(A,1) d=mean(A,2) b = 4.6667 8.0000 10.6667 c = 4.6667 8.0000 10.6667 d = 4.6667 15.0000 3.6667,矩阵的求和求积,sum(X); prod(X); sum(A); prod(A); sum(A,dim); prod(A,dim);,A=3 5 6;10 15 20;1 4 6; b=sum(A) c=sum(A,1) d=sum(A,2) b = 14 24 32 c = 14 24 32 d = 14 45 11,矩阵元素的累加和累积,设X=(x1,x2xn),U,V是与X等长的两个向量,并且 matlab用cumsum,cumprod函数表示,格式同sum,prod. 例 X=cumprod(2:2:10); X = 2 8 48 384 3840,A=3 5 6;10 15 20;1 4 6; b=cumprod(A) c=cumprod(A,1) d=cumprod(A,2),b = 3 5 6 30 75 120 30 300 720,c = 3 5 6 30 75 120 30 300 720,d = 3 15 90 10 150 3000 1 4 24,标准方差,对于有n个元素的序列x1,x2xn公式为 MATLAB提供标准方差函数为std 格式: std(X); std(A); std(A,FLAG,dim);FLAG=0时按S1计 算,FLAG=1时按S2计算.默认 FLAG=0,dim=1,例,A=3,7,14;4,7,11; y1=std(A);,y1 = 0.7071 0 2.1213,y2=std(A,0,1);,y2 = 0.7071 0 2.1213,y3=std(A,0,2);,y3 = 5.5678 3.5119,y4=std(A,1,1);,y4 = 0.5000 0 1.5000,矩阵元素排序,函数:sort 格式Y,I= sort(A,dim, mode); 若dim=1，则按列排；若dim=2，则按行排。mode值为 ascend或descend。Y是排序后的矩阵，而I记录Y中的元 素在A中位置。,例: A=3,15,7;9,6,8; Y,I=sort(A,2); Y = 3 7 15 6 8 9 I = 1 3 2 2 3 1,A=3,15,7;9,6,8;4 7 10; Y,I=sort(A,1),Y = 3 6 7 4 7 8 9 15 10,I = 1 2 1 3 3 2 2 1 3,Y,I=sort(A,2),Y = 3 7 15 6 8 9 4 7 10,I = 1 3 2 2 3 1 1 2 3,A=3,15,7;9,6,8;4 7 10; Y,I=sort(A,2,descend),Y = 15 7 3 9 8 6 10 7 4,I = 2 3 1 1 3 2 3 2 1,3.3.2曲线拟合,曲线拟合的目的是根据给定区间或区域上的有限个采样点的函数值,构造一个较简单的函数g(x)去逼近一个复杂的或未知的函数f(x)。一般采用最小二乘法,其基本思想是要求逼近函数在采样点与被逼近函数的值之差的平方和最小。 MATLAB提供最小二乘法曲线拟合函数,其格式为: p,s=polyfit(x,y,m); p为多项式的系数向量,长度 m+1 S为误差向量 函数: y=polyval(p,x); 作用:求拟合函数在采样点的函数值,例:用一个8次多项式在0,3逼近sin(x),x=0:0.1:3; y=sin(x); p=polyfit(x,y,8); y1=polyval(p,x); plot(x,y1,k*,x,y,r-);,x1=0:0.2:15; y1=sin(x1); p=polyfit(x1,y1,8); y2=polyval(p,x1); plot(x1,y2,k*,x1,y1,r-);,3.3.3数值插值,假设已知在n个自变量值x1,x2,.,xn处的函数值为y1,y2,.,yn，或称有n个样本点： (x1, y1),(x2 ,y2),.,(xn, yn)，由这些样本点的信息获得在其他点上的函数值称为插值，如果在给定点的范围内进行插值，称为内插；否 则称为外插。,3.3.3数值插值,一维数值插值: 格式:Y1=interp1(X,Y,X1,method); 函数根据X,Y的值，计算函数在X1处的值。X,Y是两个等长 的已知向量，分别描述采样点和样本值，X1是一个向量或 标量，描述欲插值的点，Y1是一个与X1等长的插值结果。 method是插值方法，允许的取值有linear（线性插值）、 nearest（最近插值),spline(三次样条插值),cubic(三次多 项式插值)，默认插值是linear(线性插值）。,例 用不同的插值方法计算sinx在/2点的值。 x=0:0.2:pi; y=sin(x); Y1=interp1(x,y,pi/2)%默认线性插值 Y2=interp1(x,y,pi/2,nearest)%最近邻插值 Y3=interp1(x,y,pi/2, spline)%样条插值 Y4=interp1(x,y,pi/2,cubic)%立方插值,例 某观测站测得某日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; XI =6.5:2:17.5 YI=interp1(h,t,XI,spline) %用3次样条插值计算,二维网格数据插值: 格式为： Z1=interp2(X,Y,Z,X1,Y1,method) 其中X,Y是两个向量，分别描述两个参数的采样点，Z是与 参数采样点对应的函数值，X1,Y1是两个向量或标量或网格 点形式的矩阵，描述欲插值的点。Z1是根据相应的插值方法 得到的插值结果。,例 某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度()。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。 命令如下： x=0:2.5:10; h=0:30:60; T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41; xi=0:10; hi=0:20:60; TI=interp2(x,h,T,xi,hi),二维一般分布数据的插值问题,格式:Z1=griddata(X,Y,Z,X1,Y1,method) 功能: X,Y,Z是已知样本点数据, X,Y,Z由向量给出时，要求长度相同；X1,Y1为需要插值处的自变量值，可以为标量、数组矩阵，z1为对应的插值结果；方法可以选：linear(线性插值，默认),cubic(三次hermite插值),v4 ,nearest(最近点等值方式)等,两者的区别: interp2的插值数据必须是矩形域，一般使用meshgid生成; griddata函数的插值数据X和Y没有那么多数据，特别是对试验中随机没有规律采取的数据进行插值具有很好的效果 .,例：clear;x,y=meshgrid(-3:0.6:3,-2:0.4:2); z=(x.2-2*x).*exp(-x.2-y.2); x1,y1=meshgrid(-3:0.2:3,-2:0.2:2); z1=interp2(x,y,z,x1,y1,spline),例：clear; x= -3+6*rand(200,1); y= -2+4*rand(200,1); z=(x.2-2*x).*exp(-x.2-y.2); x1,y1=meshgrid(-3:0.2:3,-2:0.2:2); z1=griddata(x,y,z,x1,y1,splinecubic),3.3.4多项式及其运算,多项式的建立 按降幂用向量表示.如: p(x)=2x4-3x-10 p=2 0 0 -3 -10; 多项式的求根 命令:x=roots(p); x = 1.6538 -0.1677 + 1.5050i -0.1677 - 1.5050i -1.3185,多项式求值 函数格式: polyval(p,x);x为标量或矩阵 例: x=1,2;3,4; polyval(p,x); ans = -11 16 143 490,多项式四则运算,多项式的加减法 次数相同 次数不同 例: p1(x)=5x2-2x+20, p2(x)=x3-3x2+x+13 则p1=0 5 -2 20; p2=1 -3 1 13; p=p1+p2 p= 1 2 -1 33,低阶多项式用首零填补,使其与高阶多项式有相同的阶数,多项式四则运算,多项式的乘法 函数：conv (P1,P2) P1、P2是两个多项式的系数向量,a=3 -5 2 -7 5 6;b=3 5 -3 c=conv(a,b),例：已知a(x)=3x5-5x4+2x3-7x2+5x+6和 b(x)=3x2+5x-3 求a(x)/b(x),多项式四则运算,多项式的除法 格式 :Q,r=deconv(P1,P2) 表示P1除以P2,给出商多项式Q和余数r，Q和r仍为多项式系 数向量。 deconv是conv的逆函数，即有P1=conv(P2,Q)+r。,例：已知a(x)=3x5-5x4+2x3-7x2+5x+6和 b(x)=3x2+5x-3 求a(x)/b(x),a=3 -5 2 -7 5 6;b=3 5 -3 Q,r=deconv(a,b),多项式的导数,格式: p=polyder(a);计算多项式a的导数 p= polyder(a,b); 求a*b的导数 p,q=polyder(a,b);求a/b的导数，导函数的分子存入p，分母存入q。,例：求有理分式的导数：,命令：a=3 5 0 -8 1 -5; b=10 5 0 0 6 0 0 7 -1 0 -100; p,q=polyder(a,b),多项式求值,函数：y=polyval(p,x).其中p是多项式的系数向量,x是自变量的取值(可以是向量或矩阵).,代数多项式求值,矩阵多项式求值,函数：y=polyvalm(p,x). X为方阵,设A为方阵，p=x3-5x2+8,那么polyval(p,A)的含义A.*A.*A-5*A.*A+8*ones(size(A). polyval(p,A)的含义:A*A*A-5*A*A+8*eye(size(A).,例. 已知多项式x3-3x+5,求多项式在x=5和在矩阵A=9 1;6 8;2 7的值;,p=1 0 -3 5； y=polyval(p,5) %自变量为数 y= 115,A=9 1;6 8;2 7; c=polyval(p,A) %自变量为矩阵 c= 707 3 203 493 7 327,p=1 0 -3 5; A=9 1;6 8; c=polyvalm(p,A) %自变量为方阵 c= 863 220 1320 643,c=polyval (p,A) %自变量为方阵 c= 707 3 203 493,多项式求根,函数：x=roots(p).其中p是多项式的系数向量,求得根赋给向量x。,例：求多项式p(x)=2x4-3x-10的根 命令：p=2 0 0 -3 -10; x=roots(p) x = 1.6538 -0.1677 + 1.5050i -0.1677 - 1.5050i -1.3185,多项式求根,已知多项式的根，可用函数：p=poly(x)建立多项式,例：求多项式p(x)=2x4-3x-10的根 命令：p=2 0 0 -3 -10; x=roots(p) x = 1.6538 -0.1677 + 1.5050i -0.1677 - 1.5050i -1.3185,P1=poly(x) p1 = 1.0000 -0.0000 -0.0000 -1.5000 -5.0000,3.4数值微分,在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为： dx=diff(X)：计算向量X的向前差分，dxi=X(i+1)- X(i)，i=1,2,n-1。 dx=diff(X,n)：计算X的n阶向前差分。例如， diff(X,2)=diff(diff(X)。 dx=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。,例：设x由0,2pi间均匀分布的10个点组成，求sinx的1-3阶差分。,X=linspace(0,2*pi,10);,Y=sin(X); DY=diff(Y); D2Y=diff(Y,2); D3Y=diff(Y,3);,例: 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f(x)的图像。 程序如下： f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2); g=inline(3*x.2+4*x-1)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5); x=-3:0.01:3; p=polyfit(x,f(x),5); %用5次多项式p拟合f(x) dp=polyder(p); %对拟合多项式p求导数dp dpx=polyval(dp,x); %求dp在假设点的函数值 dx=diff(f(x,3.01)/0.01; %直接对f(x)求数值导数 gx=g(x); %求函数f的导函数g在假设点的导数 plot(x,dpx,x,dx,.,x,gx,-); %作图,一、 数值积分基本原理 求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。,3.5数值积分,一元函数数值积分： (1)被积函数是一个解析式 q,n=quad(fun,a,b,tol，trace) -采用辛普森计算积分 q,n=quad8(fun,a,b,tol,trace) -采用newton cotes方法计算积分 q,n=quadl(fun,a,b,tol,trace) -采用lobatto方法计算 fun可以是字符串、内联函数或M函数名， a,b是定积分的下限和上限，tol表示绝对误差限，默认10-6， trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值 ,返回参数q即定积分值，n为被积函数的调用次数。,(1) 建立被积函数文件ex.m。 function ex=ex(x) ex=exp(-x.2); 然后在命令窗口中，输入 format long I=quad(ex,0,1) I = 0.74682418072642,例 求定积分,I=quadl(ex,0,1),I = 0.74682413398845,(2)使用内联函数,g=inline(exp(-x.2); I=quadl(g,0,1);,format short,(2)被积函数由一个表格定义,函数：trapz(X,Y) 。其中向量X,Y定义函数关系Y=f(X)。,例：X=1:0.01:2.5; Y=exp(-X); %生成函数关系数据向量 trapz(X,Y) ans = 0.28579682416393,（3）二重定积分的数值求解 函数： I=dblquad(f,a,b,c,d,tol,trace) 该函数求f(x,y)在a,bc,d区域上的二重定积分。参数 tol，trace的用法与函数quad完全相同。,注意本函数不允许返回被积函数的调用次数，如果需要，可以在被积函数中设置一个计数变量，从而统计被积函数的调用次数,例 计算二重定积分,(1) 建立一个函数文件fxy.m： function f=fxy(x,y) global ki; ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.2/2).*sin(x.2+y); (2) 调用dblquad函数求解。 global ki;ki=0; I=dblquad(fxy,-2,2,-1,1) I = 1.57449318974494 ki = 1050,3.6离散傅立叶变换 一、 离散傅立叶变换的实现 一维离散傅立叶变换函数，其调用格式与功能为： fft(X) 返回向量X的离散傅立叶变换。设X的长度(即元素个数)为 N，若N为2的幂次，则为以2为基数的快速傅立叶变换，否 则为运算速度很慢的非2幂次的算法。 对于矩阵X，fft(X)应用于矩阵的每一列。,(2) fft(X,N) 计算N点离散傅立叶变换。它限定向量的长度为N，若X的长 度小于N，则不足部