《Matlab应用51》由会员分享，可在线阅读，更多相关《Matlab应用51（31页珍藏版）》请在金锄头文库上搜索。
1、3、数值计算、数值计算v矩阵特征值和恰定方程求解矩阵特征值和恰定方程求解v稀疏矩阵的存贮与运算稀疏矩阵的存贮与运算v数值微分和积分数值微分和积分v常微分方程求解常微分方程求解数组数组数组运算矩阵运算s./B,B.s 标量s分别除B中对应的元素s*inv(B) B阵的逆乘A.n A的每个元素自乘n次An A必为方阵，自乘n次p.A 以p为底，分别以A的每个元素为指数求幂值pA A必为方阵，标量的矩阵乘方A.*B 对应元素之积A*B 内维相同矩阵的乘积A./B A的元素被B的对应元素除A/B A右除B (=A*inv(B)B.A （一定与 A./B 结果相同）BA A左除B（一般与 A/B 结果不
2、同）exp(A) 以自然对数e为底，分别以A的元素为指数，求幂expm(A) A的矩阵指数函数log(A) 对A的各元素求对数logm(A) A的矩阵对数函数sqrt(A) sqrtm(A)f(A)funm(A)矩阵特征值函数及其功能矩阵特征值函数及其功能inv矩阵求逆trace求矩阵的迹det求矩阵行列式cond条件数rank求矩阵的秩condest估计范-1条件数poly矩阵特征多项式condeig特征值有关条件数norm 或 normest矩阵和向量范数eig 或 eigs求特征值、特征向量3-1、矩阵特征参数运算、矩阵特征参数运算3-1、矩阵特征参数运算、矩阵特征参数运算1、逆运算：i
3、nv a=1 2 3;5 8 6;7 2 9, b=inv(a)2、特征值计算、特征值计算eig（全元素阵）d=eig(a) %仅计算a的特征值c,d=eig(a) %满足a*c=c*d; c,d=eig(a,nobalance) %a中有截断误差量级数eigs（大型稀疏矩阵） d=eigs(a),c,d=eigs(a) %与eig同 c0,d0=eigs(a,2,10) %计算a的 2个在10附近的特征向量和特征值 eigs(a,2,lm) % lm最大幅值; sm最小幅值;lr最大实部;sr最小实部;be谱的两端矩阵中有元素与截断误差相当时的特征值问题 A=3 -2 -0.9 2*eps;
4、-2 4 -1 -eps; -eps/4 eps/2 -1 0; -0.5 -0.5 0.1 1 %A中含有形式“零”元素 V1,D1=eig(A);ER1=A*V1-V1*D1 %ER1 最大值0.5216 V2,D2=eig(A,nobalance);ER2=A*V2-V2*D2 % ER2 最大值1.0e-014%eig指令执行过程中首先调用使原矩阵各元素大致相当的“平衡”程序，使原方阵中的小元素作用被放大。指令eig与eigs的比较 rand(state,1),A=rand(100,100)-0.5; % 产生100阶随机方阵 t0=clock;V,D=eig(A);T_full=et
5、ime(clock,t0) % 指令eig的运算时间v% T_full = 0.1570 options.tol=1e-8;options.disp=0; % 设置eigs运算精度,不显示中间迭代结果 t0=clock;V2,D2=eigs(A,1,lr,options);T_part=etime(clock,t0) % 计算最大实部特征值和特征向量时指令eigs的运算时间vT_part = 0.5160veigs比eig用时长，主要用于大型稀疏矩阵稀疏矩阵稀疏矩阵3-1、矩阵特征参数运算、矩阵特征参数运算3、行列式运算（方阵）： det(a)4、秩 rank(a)5、迹 trace(a) %
6、对比 trace(d)6、范数运算 e=norm(a) %(2-范数)； norm(a,1) %(1-范数) ；norm(a,inf) %(无穷范数)； norm(a,fro) %(Frobenius 范数) 7、条件数运算 d=cond(a) % 2-范数的条件数 d=cond(a,1) %1范数的条件数，或2、inf 相应范数的条件数) 线性方程组的求解线性方程组的求解v含有n个未知数的n个方程构成的方程组称为恰定方程组，当det(A)0时存在唯一解；当 det(A)=0时由增广矩阵作行的最简形式,从中找出方程组的基础解系及特解v如： a a1,11,1x x1 1 + a + a1,21
7、,2x x2 2 + + a + + a1,n1,nx xn n = = b b1 1 a a2,12,1x x1 1 + a + a2,12,1x x2 2 + + + + a a2,n2,nx xn n = b = b2 2 a an,1n,1x x1 1 + a + an,1n,1x x2 2 + + + + a an,nn,nx xn n = = b bn n 记为 AX=b ，其中： 恰定方程组的解恰定方程组的解v线性代数中最常见的该方程的解法有：1、利用公式求解： 式中Aij是元素aij的余子式2、利用逆阵求解法，即：x=A-1b3、利用Gaussian消元法4、利用LU法求解 对
8、于维数不高、条件数不大的矩阵，以上4种解法结果相近。对于高阶大型矩阵主要应用LU分解法进行求解。恰定方程组的解恰定方程组的解vLU分解（也称三角分解），满足：分解（也称三角分解），满足：LU=PA式中，式中，L为主对角元为为主对角元为1的下三角阵的下三角阵,U上三角阵上三角阵,P是由是由0或或1组成组成的交换矩阵。的交换矩阵。 a=4 2 -1;3 -1 2;11 3 1 ; b=2 10 8 l,u,p=lu(a) %满足l*u=p*a v行列式和逆：行列式和逆：detA=det(P-1LU)=detU= A-1=U-1L-1P det(a); inv(a); %x=inv(a)*bv在在m
9、atlab中求解中求解Ax=b可以直接输入指令可以直接输入指令： Ab %x=Ab其他矩阵分解计算其他矩阵分解计算v正交分解正交分解： q,r=qr(a) %r为上三角阵，q正交阵，满足a=q*r q,r,p=qr(a) %r上三角阵，q正交阵，p交换矩阵，满足a*p=q*rv特征值分解特征值分解：d,e=eig(a) %d以a的特征向量作为列向量组成矩阵，e以特征向量作为主对角元素的对角阵，满足a*d=d*ed,e=eig(a,b) %广义特征值分解，满足a*d=b*d*e vChollesky分解：d=chol(a) %a必须为对称正定矩阵，满足a=d*d v奇异值分解奇异值分解：u,s,
10、v=svd(a) %s为对角阵，u、v正交阵，满足a=u*s*vu,s,v=svd(a,0) %奇异值最佳分解 求解非齐次线性方程组求解非齐次线性方程组例例：解方程组：1、LU分解法求解A=4 2 -1;3 -1 2;11 3 1;b=2 10 8;L,U=lu(A)X=U(Lb)解:2、正交分解法求解Q,R=qr(A)X=R(Qb)3、直接求解 X=Ab练习练习v求解线形方程组：习题答案习题答案 x=-2:0.2:2 ;y=x; X,Y=meshgrid(x,y) %以0.2步长，-2 2区间划分计算网格点z=X.*exp(-X.2-Y.2) %计算网格点上z值的大小px,py=gradie
11、nt(z,.2,.2) %0.2为网格步长，求z的导数 cz=linspace(min(min(z),max(max(z),8) %通过8条等值线描绘z的大小 contourf(X,Y,z,cz,k-) %contour(z)无格式hold on %保持图形 AE=sqrt(px.2+py.2);px=px./AE;py=py./AE %箭头等长化，可略 quiver(X,Y,px,py,0.7) %方向导数可视化3-2、稀疏矩阵、稀疏矩阵vn个未知数的线性方程的系数组成n*n的矩阵（解方程需n2内存，n3时间）。v通常，一个未知数只与少量其他变量有关。多数系数为0，关系矩阵为稀疏矩阵。v不存
12、贮0元素，也不对其操作，从而节省内存和运算时间。稀疏存储特点稀疏存储特点v两种存贮矩阵方式：全元素存储和稀疏存储。v稀疏存储特点：（1）只存储矩阵的非0元素，并且存储按列进行；（2）每列用一个实（复）数组记述非0元素值，再用一个整数数组记述相应的非0元素的行下标v如果某个矩阵以稀疏方式存储，由它参与的运算结果仍以稀疏方式存储，除非运算本身使稀疏性消失。v两种存储方式转换指令： c=sparse(b) % 满矩阵转换成稀疏矩阵 d=full(c) %稀疏矩阵转换成满矩阵常用稀疏矩阵指令常用稀疏矩阵指令v稀疏矩阵生成： a=eye(8);b=speye(size(a) 或 b=speye(8) %
13、稀疏单元阵 b=sprand(a)或 b=sprand(8,6,0.1) %稀疏随机阵,0.1非0元素密度v返回矩阵非0元素位置、值： i,j=find(b) %返回非0元素位置-i行j列;i,j,v=find(d) % i行j列v元素值v稀疏矩阵操作： b=sprand(8,6,0.2),nz=nnz(b) %返回非0元素的个数 dz=nnz(b)/prod(size(b) %返回非0元素的密度 a=spones(b) %生成与b同样矩阵，非0元素均为1v稀疏矩阵判断： issparse(a) %稀疏矩阵=1，全元素阵=0常用稀疏矩阵指令常用稀疏矩阵指令nonzeros矩阵的非0元素值sym
14、rcm反向Cuthill-McKee排序find找非0元素下标sprank结构秩spones用1置换非0元素etree消去树spfun求各非0元素的函数值spy画稀疏结构图colmmd列最小度排序spalloc为非0元素配置内存colperm计算列排序置换向量nnz矩阵的非0元素总数dmperm矩阵Dulmage-Mendelsohn分解nzmax指定存放非0元素所需内存randperm随机置换向量treelayout展开树symmmd对称最小度排序treeplot画树图3-3、数值微分和积分、数值微分和积分1、数值微分：vdiff(a) %求数组a相邻行元素间的1阶（按列方向）差分vdiff
15、(a,n) %求数组a相邻行元素间的n阶（按列方向）差分vdiff(a,n,1) %求数组a相邻行元素间的n阶（按列方向）差分vdiff(a,n,2) %求数组a相邻列元素间的n阶（按行方向）差分 F=1,2,3;6,5,8;7,8,9 ;Dx=diff(F) %列方向1阶差分 Dx_2=diff(F,2) %列方向2阶差分 Dx_2=diff(F,2,1) %列方向2阶差分 Dx_2=diff(F,2,2) %行方向2阶差分3-3、数值微分和积分、数值微分和积分2、近似梯度：vax,ay=gradient(a) %x-行向梯度，y-列向梯度vax,ay=gradient(a,h1,h2) %
16、h1行向步长，h2列向步长 FX_2,FY_2=gradient(F,0.5) %行列向步长均取0.5 FX_2,FY_2=gradient(F,0.5,5) %行步长取0.5，列步长取53、方向导数的可视化:vquiver(X,Y,U,V,scale) %u,v为平面坐标位置x,y点的方向导数，scale为绘图比例尺vquiver3(X,Y,Z,U,V,W,scale) %u,v,w为平面坐标位置x,y,z点的方向导数，scale为绘图比例尺，比例尺越大，箭头越长 quiver(2,3,-3,-5,0.6) %(2,3)点开始向(-3,-5)方向画箭头，比例尺为0.6 quiver3(2,3
17、,5,-3,-5,-8,0.6) %(2,3,5)点开始向(-3,-5,-8)方向画箭头，比例尺为0.6 梯度和方向导数的关系：函数在某点的梯度是这样一个向量，它的方向与取得最大方向导数的方向一致，而它的模为方向导数的最大值.数值微分应用举例数值微分应用举例v研究偶极子研究偶极子(Dipole)的电势（的电势（Electric potential）和电场强度（）和电场强度（Electric field density）。设在）。设在(a,b)处有电荷处有电荷+q，在，在(-a,-b)处有电荷处有电荷-q。那么在电荷所在平。那么在电荷所在平面上任何一点的电势和场强分别为面上任何一点的电势和场强分
18、别为 。其中。其中 ： 。又设电荷。又设电荷q=2*10-6, a=1.5, b=-1.5。 clear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5; %清内存变量，清图形窗口，输入参数 x=-6:0.6:6 ;y=x; X,Y=meshgrid(x,y) %设置坐标网点 rp=sqrt(X-a).2+(Y-b).2);rm=sqrt(X+a).2+(Y+b).2);V=q*k*(1./rp-1./rm); Ex,Ey=gradient(-V,0.6,0.6);%计算电势V和场强 quiver(X,Y,Ex,Ey,0.7) %场强可视化 cv=linspace(min(min(
19、V),max(max(V),49); %设定画电势图等势线值 AE=sqrt(Ex.2+Ey.2);Ex=Ex./AE;Ey=Ey./AE; %场强归一化，使箭头等长 figure(2);contourf(X,Y,V,cv,k-) %另开2号图形窗口，用黑实线画填色等位线图 hold on; quiver(X,Y,Ex,Ey,0.7) %在等势线图上添加场强图 plot(a,b,wo,a,b,w+);plot(-a,-b,wo,-a,-b,w-) %用白线画正负电荷位置练习练习v求解线形方程组：v已知z(x,y)满足函数： 图视化表示函数z在x=-2,2，y=-2,2区间的大小和方向导数（网格
20、步长取0.2 ）习题答案习题答案a=16 2 3 4;5 11 10 8;12 11 6 9;8 4 6 5, b =190;302;309;194,x=ab x=-2:0.2:2 ;y=x; X,Y=meshgrid(x,y) %以0.2步长，-2 2区间划分计算网格点z=X.*exp(-X.2-Y.2) %计算网格点上z值的大小px,py=gradient(z,.2,.2) %0.2为网格步长，求z的导数 cz=linspace(min(min(z),max(max(z),8) %通过8条等值线描绘z的大小 contourf(X,Y,z,cz,k-) %contour(z)无格式hold
21、on %保持图形 AE=sqrt(px.2+py.2);px=px./AE;py=py./AE %箭头等长化，可略 quiver(X,Y,px,py,0.7) %方向导数可视化常见一元函数数值积分指令常见一元函数数值积分指令quad采用递推自适应Simpson法计算积分quadl采用递推自适应Lobatto法计算积分trapz采用梯形法求定积分，速度快精度差cumtrapz采用梯形法求区间上的积分曲线，速度快精度差sum等宽矩形法求定积分，速度快精度很差cumsum采用等宽矩形法求区间上的积分曲线，精度很差fnint利用样条函数求不定积分数值积分举例数值积分举例v例：求sin(x)在区间0,1
22、0上的积分。 x=0:0.1:10;y=sin(x); z1=cumsum(y)*0.1%矩形法求积分曲线 plot(x,y,r-,x,z1,g*) %矩形法积分曲线图 z2=trapz(y)*0.1 %梯形法求定积分 z3=quad(sin(x),0,10) %递推自适应Simpson法计算积分 z4=quadl(sin(x),0,10) %递推自适应Lobatto法计算积分数据显示格式数据显示格式vshort (缺省) 小数点后4位, long 15位数字，vshort e 5位科学计数， long e 15位科学计数v通过File下子菜单Preferences选择；用format加格式选
23、择矩阵的范数和条件数矩阵的范数和条件数v范数和条件数对求解过程中误差的放大现象进行定量描述v矩阵的范数确定了被该矩阵相乘向量长度的最大可能的放大倍数，2-范数的定义为：v矩阵的条件数描述方程Ax=b的解对b中误差的灵敏度，定义为： 条件数总大于等于1，正交矩阵的条件数为1，奇异矩阵的条件数是。条件数大到不可接受的矩阵称为病态矩阵，如Hilbert阵。a1=1 -2 2 3;a2=-2 4 -1 3;a3=-1 2 0 3;a4=0 6 2 3;a5=2 -6 3 4;A=a1 a2 a3 a4 a5例例：问向量组a1=（1 -2 2 3），a2=（-2 4 -1 3），a3=（-1 2 0 3），a4=（0 6 2 3），a5=（2 -6 3 4）是否线性相关？若线性相关求一个最大无关组解：构成矩阵format ratrank(A) %ans=3 A=1 2 2 1; 2 1 -2 -2; 1 -1 -4 -3; %输入系数矩阵 rank(A) %计算A的秩 format rat %指定有理化格式输出 B=null(A,r) %解线性齐次方程组的有理解 根据B的形式（由A的秩及未知数个数决定）确定变量数syms k1 k2 %定义符号变量-k1,k2分别对应的列X=k1*B(:,1)+k2*B(:,2) %给出方程组的通解pretty(X) %美化输出格式
