Matlab软件使用简介(72页).doc

-Matlab软件使用简介一、基本内容1. 向量的产生基本格式：start: step: end（step缺省时为1）x=1:5x=0:0.1:2*pidot 向量点集 cross 向量叉集 .* 对应元素乘积Matlab的每条命令后，若为逗号或无标点符号，则显示命令的结果；若命令后为分号，则禁止显示结果。linspace 线性等分向量 a=linspace(0,2*pi,100) logspace 对数等分向量 a=logspace(0,2*pi,100)100为插入点数（99等分），默认值为100.2. 矩阵的输入约定：（1）元素之间用空格或逗号隔开； （2）用中括号（ ）把所有元素括起来； （3）用分号（；）说明行结束。矩阵输入时，按Enter键表示开始输入新的一行，输入矩阵时，要求所有的行具有相同的列。例如：a=1 2 3; 4,5 6; 7 8,9该矩阵一直保存在工作空间，直至被修改。x=1:6y=sin(x)z=cos(x)b=x;y;z特殊矩阵：a= 产生一个空矩阵b=zeros(m,n)产生一个m行、n列的零矩阵c=ones(m,n)产生一个m行n列的元素全为1的矩阵d=eye(m,n)产生一个m行、n列的单位矩阵3. 大矩阵中抽取小矩阵c1=b(:,1:2); c2=b(:,5:6) ; d=c1,c2d=c1;c2或d=b(:,1:2,5:6)4. 固定变量(1) ans: 在没有定义变量名时，系统默认变量名为ans；(2) eps: 容许误差，非常小的数；(3) pi: 即圆周率；(4) i: 虚数单位；(5) inf: 表示正无穷大，由1/0运算产生；(6) NaN: 表示不定值，由inf/inf或0/0运算产生。5. 基本运算5.1 算算术运算符+ 加 - 减 * 矩阵乘法 .* 数组乘（对应元素相乘） 矩阵幂 . 数组幂（各个元素求幂） ./ 数组除（对应元素除） 左除或反斜杠 / 右除或斜杠 如果a为一个非奇异矩阵：ab = inv(a)*b： 表示a*x=b的解；a/b = b*inv(a)： 表示项x*a=b的解。例如：a=1 2 3; 4 2 6; 7 4 9 b=4;1;2 ab5.2 2.2 关系运算符= = 等号 = 不等号 < 小于 > 大于 <= 小于或等于 >= 大于或等于 5.3 逻辑运算符& 逻辑与 | 逻辑或 逻辑非 xor 异或 any 有非零元则为真 all 所有元素均非零则为真 6. 矩阵的基本操作' 转置inv 矩阵求逆det 行列式的值 v d=eig(a) 特征值和特征向量rank 秩 trace 迹size 矩阵的行数和列数diag 对角矩阵和矩阵对角线 fliplr 从左自右翻转矩阵 flipud 从上到下翻转矩阵 roy90 矩阵翻转90度 tril 矩阵的下三角 triu 矩阵的上三角 7. 常用函数（1）clc清除指令窗口（2）clear从内存中清除变量和函数（3）who列出工作内存中的变量名（4）whos列出工作内存中的变量详细信息（4）pause暂停二、绘图函数1. plot是绘制一维曲线的基本函数，但在使用此函数之前，我们需先定义曲线上每一点的x 及y座标。下例可画出一条正弦曲线：x=linspace(0, 2*pi, 100); % 100个点的x座标 y=sin(x); % 对应的y座标 plot(x,y); 若要画出多条曲线，只需将座标对依次放入plot函数即可plot(x, sin(x), x, cos(x); 若要改变颜色，在座标对后面加上相关字串即可： plot(x, sin(x), 'c', x, cos(x), 'g'); 若要同时改变颜色及图线型态（Line style），也是在座标对后面加上相关字串即可： plot(x, sin(x), 'co', x, cos(x), 'g*') plot绘图函数的参数符号颜色符号线形符号线形y黄色.点-虚线m洋红色o圆d菱形c青色×叉号>向右三角形r红色+加号<向左三角形g绿色*星号s正方形b蓝色-实线P正五角星w白色:点线h正六角星k黑色-.点划线图形完成后，我们可用axis(xmin,xmax,ymin,ymax)函数来调整图轴的范围： axis(0, 6, -1.2, 1.2); 此外，Matlab也可对图形加上各种注解与处理： xlabel('Input Value'); % x轴注解 ylabel('Function Value'); % y轴注解 title('Two Trigonometric Functions'); % 图形标题 legend('y = sin(x)','y = cos(x)'); % 图形注解 grid on; % 显示格线 注意：如果标题太长，需要分行，请用如下格式str='你需要多少行？''3行可以吗''试试就知道了'title(str)  我们可用subplot来同时画出数个小图形于同一个视窗之中： subplot(2,2,1); plot(x, sin(x);title('y=sinx') subplot(2,2,2); plot(x, sin(2*x); title('y=sin2x')subplot(2,2,3); plot(x, sin(3*x); title('y=sin3x')subplot(2,2,4); plot(x, sin(4*x); title('y=sin4x')Matlab还有其他各种二维绘图函数，以适合不同的应用。bar 长条图errorbar 图形加上误差范围 fplot 较精确的函数图形 polar 极座标图hist 直方图rose 绘制角度直方图（玫瑰图）stairs 阶梯图stem 针状图 fill 多边形填充图feather 羽毛图compass 罗盘图quiver 向量场图 以下我们针对每个函数举例。 （1）当数据点数量不多时，长条图是很适合的表示方式： close all; % 关闭所有的图形视窗 x=1:10; y=rand(size(x); bar(x,y); rand 均匀颁随机数和数组（2）如果已知数据的误差量，就可用errorbar来表示。下例以单位标准差来做资的误差量： x = linspace(0,2*pi,30); y = sin(x); e = std(y)*ones(size(x); % std求标准差errorbar(x,y,e) （3）对於变化剧烈的函数，可用fplot来进行较精确的绘图，会对剧烈变化处进行较密集的取样，如下例： fplot('sin(1/x)', 0.02 0.2); % 0.02 0.2是绘图范围 （4）若要产生极座标图形，可用polar： theta=linspace(0, 2*pi); r=cos(4*theta); polar(theta, r); （5）对於大量的数据，我们可用hist来显示资料的分布情况和统计特性。下面几个命令可用来验证由randn产生的正态分布： x=randn(5000, 1); % 产生5000个均值为0，方差为1的正态分布随机数 hist(x,20); % 20代表长条的个数   （6）rose和hist很接近，只不过是将数据大小视为角度，数据个数视为距离，并用极座标绘制表示： x=randn(1000, 1); rose(x); （7）stairs可画出阶梯图： x=linspace(0,10,50); y=sin(x).*exp(-x/3); stairs(x,y); （8）stems可产生针状图，常被用来绘制数位讯号： x=linspace(0,10,50); y=sin(x).*exp(-x/3); stem(x,y); （9）stairs将资料点视为多边行顶点，并将此多边行涂上颜色： x=linspace(0,10,50); y=sin(x).*exp(-x/3); fill(x,y,'b'); % 'b'为蓝色 （10）feather将每一个数据点视复数，并以箭号画出： theta=linspace(0, 2*pi, 20); z = cos(theta)+i*sin(theta); feather(z); （10）compass和feather很接近，但每个箭号的起点都在圆点： theta=linspace(0, 2*pi, 20); z = cos(theta)+i*sin(theta); compass(z); 2. plot3用来绘制空间曲线格式：plot(x,y,z, cs) 其中若x,y,z为长度相等的向量，则在空间表示一条曲线；若x,y,z为m*n的矩阵，则每一列对应一条曲线。Cs表示曲线的颜色、线形等性质。例螺旋线：t=0:pi/50:10*pi; x=sin(t); y=cos(t); plot3(x,y,t)例多条曲线：x=-3:0.1:3;y=1:0.1:5; X,Y=meshgrid(x,y); Z=(X+Y).2; plot3(X,Y,Z)（meshgrid(x,y)的作用是产生一个以向量x为行、向量y为列的矩阵）3. surf用来绘制空间曲面格式：surf(x,y,z)x,y,z为矩阵，显式函数的图形。例：显示的图形。n=-8:0.5:8;m=n' x y=meshgrid(n,m);t=sqrt(x.2+y.2)+eps;z=sin(t)./t;surf(x,y,z)4. mesh用来绘制网格曲面格式同上5. scatter离散点绘图格式：scatter(x,y,cs)例: x=0:pi/10:2*pi;y=sin(x);scatter(x,y,'+')6. 其它函数（1）hold on: 保持当前图形, 以便继续在当前图上画图； hold off：释放当前窗口（2）figure(n): 打开或创建图形窗口n三、微积分内容1. limit求极限（必须首先声明变量）格式：limit(F,x,a)求limit(F,a) 求limit(F) 求limit(F,x,a,'right')的右极限，即limit(F,x,a,'left')的右极限，即例：syms x; %声明x为符号变量（cosn c）limit(1+1/x)x,x,inf)；limit(sin(x)/x,x,0); limit(sin(x)/x,0); limit(sin(x)/x);syms h; limit(sin(x+h)-sin(x)/h,h,0) 2. diff求导数格式：diff(F)求函数的一阶导数diff(F, v)求函数对变量一阶偏导数diff(F,n)求函数的阶导数diff(F, v, n)求函数对变量的阶偏导数例：diff(x*y*z); diff(x*y*z,2) diff(diff(x*y*z,x),y) diff(x2*cos(x),x,2)3. 积分函数（1）符号积分int(f) 传回f对预设独立变数的积分值int(f,'t') 传回f对独立变数t的积分值int(f,a,b) 传回f对预设独立变数的积分值，积分区间为a,b，a和b为数值式int(f,'t',a,b) 传回f对独立变数t的积分值，积分区间为a,b，a和b为数值式int(f,'m','n') 传回f对预设变数的积分值，积分区间为m,n，m和n为符号式例 syms x y int(x) int(x,x) int(x*y,x)Int(x,1,2) int(x,x,1,2) 注意：int(x,x,a,b) int(x,x,'a','b')的区别。例：计算二重积分,其中由轴，轴及围成。syms x y c=int(int(x+y,y,0,1-x),x,0,1)vpa(c,10) %给出数值型符号结果（2）数值积分trapz(x,y) 梯形法沿列方向求函数Y关于自变量X的积分cumtrapz(x,y) 梯形法沿列方向求函数Y关于自变量X的累计积分quad(fun,a,b) 采用递推自适应Simpson法计算积分（低阶方法）quadl(fun,a,b) 采用递推自适应Lobatto法求数值积分（高阶方法）dblquad(fun,xmin,xmax,ymin,ymax,tol) 矩形区域二重数值积分triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol) 长方体区域三重数值积分例：计算定积分的值。符号积分：syms x; int(exp(-sin(x),x,0,4) 数值积分：参看ex01.m例：计算积分。quad('sin(x)./x',0,2) quadl('sin(x)./x',0,2)注意：quad函数和quadl函数求积分时，其中求解变量是以向量来计算的，所以这里的sin(x)和x都是向量，因此要用点除；同时被积函数用' '括起来。例：计算二重积分，dblquad('x*y',0,1,0,1)例：triplequad(' y*sin(x)+z*cos(x)',0,pi,0,1,-1,1)注意：最新版的2009a有关于一般区域二重积分的计算函数quad2d，如果用的不是2009a，那么你可以利用NIT工具箱里的quad2dggen函数。quad2dggen(fun,xmin,xmax,ymin,ymax,tol) 任意区域上二重数值积分如果既没有NIT工具箱用的也不是2009a，可以利用两次quadl函数。例：计算二重积分，quadl(x) arrayfun(xx) quadl(y) xx+y,0, xx),x),0,1)% 创建函数 4. 无穷级数格式：r = symsum(s) r = symsum(s,v) r = symsum(s,a,b) r = symsum(s,v,a,b) syms k n x; symsum(k2)symsum(k2,1,2) symsum(k2,1,inf)5. 方程求解基本函数:solve(eq)solve(eq,var)solve(eq1,eq2,.,eqn)g = solve(eq1,eq2,.,eqn,var1,var2,.,varn)例 solve('2*x+1=0','x')或syms x; solve('2*x+1=0') S = solve('x + y = 1','x - 11*y = 5') S.y S.xroots(vec) vec为系数组成的向量例 roots(2 1)线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断：若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；若系数矩阵的秩r<n，则可能有无穷解；线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解；其特解的求法属于解的第一类问题，通解部分属第二类问题。5.1 求线性方程组的唯一解或特解（第一类问题）这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 直接法；另一类是解大型稀疏矩阵 迭代法。5.1.1 利用矩阵除法求线性方程组的特解（或一个解）方程：AX=b解法：X=Ab例 求方程组的解。解：>>A=5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5; B=1 0 0 0 1'R_A=rank(A) %求秩X=AB %求解运行后结果如下R_A = 5X = 2.2662 -1.7218 1.0571 -0.5940 0.3188这就是方程组的解。或用函数rref求解：C=A,B %由系数矩阵和常数列构成增广矩阵CR=rref(C) %将C化成行最简行R = 1.0000 0 0 0 0 2.2662 0 1.0000 0 0 0 -1.7218 0 0 1.0000 0 0 1.0571 0 0 0 1.0000 0 -0.5940 0 0 0 0 1.0000 0.3188则R的最后一列元素就是所求之解。例 求方程组的一个特解。解：>>A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;>>B=1 4 0'>>X=AB %由于系数矩阵不满秩，该解法可能存在误差。X = 0 0 -0.5333 0.6000（一个特解近似值）。若用rref求解，则比较精确：>> A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;B=1 4 0'>> C=A,B; %构成增广矩阵>> R=rref(C)R = 1.0000 0 -1.5000 0.7500 1.2500 0 1.0000 -1.5000 -1.7500 -0.2500 0 0 0 0 0由此得解向量X=1.2500 0.2500 0 0（一个特解）。5.1.2 利用矩阵的LU、QR和cholesky分解求方程组的解（1）LU分解：LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。则：A*X=b 变成L*U*X=b所以X=U(Lb) 这样可以大大提高运算速度。命令 L，U=lu (A)例 求方程组的一个特解。解：A=4 2 -1;3 -1 2;11 3 -1;B=5 7 14'D=det(A)L,U=lu(A)X=U(LB)见ex04.m例 求方程组的一个特解。解：>>A=4 2 -1;3 -1 2;11 3 0;>>B=2 10 8'>>D=det(A)>>L,U=lu(A)>>X=U(LB)显示结果如下：D = 0L = 0.3636 -0.5000 1.0000 0.2727 1.0000 0 1.0000 0 0U = 11.0000 3.0000 0 0 -1.8182 2.0000 0 0 0.0000Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.018587e-017.> In D:Matlabpujunlx0720.m at line 4X = 1.0e+016 * -0.4053 1.4862 1.3511说明：结果中的警告是由于系数行列式为零产生的。可以通过A*X验证其正确性。（2）Cholesky分解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即： 其中R为上三角阵。方程 A*X=b 变成 所以 R = chol(X)（3）QR分解对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR方程 A*X=b 变形成 QRX=b所以 X=R(Qb)上例中 Q, R=qr(A)X=R(QB)说明：这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。5.2 求线性齐次方程组的通解在Matlab中，函数null用来求解零空间，即满足A·X=0的解空间，实际上是求出解空间的一组基（基础解系）。格式 z = null % z的列向量为方程组的正交规范基，满足。 % z的列向量是方程AX=0的有理基例 求解方程组的通解：解：>>A=1 2 2 1;2 1 -2 -2;1 -1 -4 -3;>>format rat %指定有理式格式输出>>B=null(A,'r') %求解空间的有理基运行后显示结果如下：B = 2 5/3 -2 -4/3 1 0 0 1或通过行最简行得到基：>> B=rref(A)B = 1.0000 0 -2.0000 -1.6667 0 1.0000 2.0000 1.3333 0 0 0 0即可写出其基础解系（与上面结果一致）。写出通解：syms k1 k2X=k1*B(:,1)+k2*B(:,2) %写出方程组的通解pretty(X) %让通解表达式更加精美运行后结果如下：X = 2*k1+5/3*k2 -2*k1-4/3*k2 k1 k2% 下面是其简化形式 2k1 + 5/3k2 -2k1 - 4/3k2 k1 k2 5.3 求非齐次线性方程组的通解非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。因此，步骤为：第一步：判断AX=b是否有解，若有解则进行第二步第二步：求AX=b的一个特解第三步：求AX=0的通解第四步：AX=b的通解= AX=0的通解+AX=b的一个特解。例 求解方程组解：见ex06.m例 求解方程组的通解：解法一：在Matlab编辑器中建立M文件如下：A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0'B=A b;n=4;R_A=rank(A)R_B=rank(B)format ratif R_A=R_B&R_A=n X=Abelseif R_A=R_B&R_A<n X=Ab C=null(A,'r')else X='Equation has no solves'end运行后结果显示为：R_A = 2R_B = 2Warning: Rank deficient, rank = 2 tol = 8.8373e-015.> In ex07 at 13X = 0 0 -8/15 3/5 C = 3/2 -3/4 3/2 7/4 1 0 0 1 所以原方程组的通解为X=k1+k2+解法二：用rref求解A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0'B=A b;C=rref(B) %求增广矩阵的行最简形，可得最简同解方程组。运行后结果显示为：C = 1 0 -3/2 3/4 5/4 0 1 -3/2 -7/4 -1/4 0 0 0 0 0 对应齐次方程组的基础解系为：， 非齐次方程组的特解为：所以，原方程组的通解为：X=k1+k2+。5.4 非线性方程的解非线性方程的标准形式为f(x)=0函数 fzero格式 x = fzero (fun,x0) %用fun定义表达式f(x)，x0为初始解。x = fzero (fun,x0,options)x,fval = fzero() %fval=f(x)x,fval,exitflag = fzero()x,fval,exitflag,output = fzero()说明 该函数采用数值解求方程f(x)=0的根。例 求的根解： fun='x3-2*x-5'z=fzero(fun,2) %初始估计值为25.5 非线性方程组的解非线性方程组的标准形式为：F(x) = 0其中：x为向量，F(x)为函数向量。函数 fsolve格式 x = fsolve(fun,x0) %用fun定义向量函数，其定义方式为：先定义方程函数function F = myfun (x)。F =表达式1；表达式2；表达式m %保存为myfun.m，并用下面方式调用：x = fsolve(myfun,x0)，x0为初始估计值。x = fsolve(fun,x0,options)x,fval = fsolve() %fval=F(x)，即函数值向量x,fval,exitflag = fsolve()x,fval,exitflag,output = fsolve()例 求下列方程组的根解：化为标准形式 设初值点为x0=-5 -5。先建立方程函数文件，并保存为myfun.m：function F = myfun(x)F = 2*x(1) - x(2) - exp(-x(1); -x(1) + 2*x(2) - exp(-x(2);然后调用优化程序x0 = -5; -5; % 初始点options=optimset('Display','iter'); % 显示输出信息x,fval = fsolve(myfun,x0,options) 6. 常微分方程求解求微分方程的解析解：r = dsolve('eq1,eq2,.', 'cond1,cond2,.', 'v')r = dsolve('eq1','eq2',.,'cond1','cond2',.,'v')dsolve('eq1,eq2,.', 'cond1,cond2,.', 'v')注意：在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分。任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省。例如，微分方程 应表达为：D2y=0.dsolve('Dy=x','x')例:求微分方程的特解.y=dsolve('D2y+4*Dy+29*y=0' ,'x')y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')微分方程数值解法：t,x=ode23(xfun,t0 tf,y0,tol) t,x=ode45(xfun, t0 tf,y0,tol)其中：xfun为以待解方程或方程组写成的M函数 t0 tf:自变量初值和终值 yo:函数的初值 tol:设置容许误差（相对误差和绝对误差）求解过程：options = odeset('RelTol',1e-4,'AbsTol',1e-4 1e-4 1e-5);T,Y = ode45(rigid,0 12,0 1 1,options);plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')四、数据处理插值和拟合都是数据优化的一种方法，当实验数据不够多时经常需要用到这种方法来画图。在matlab中都有特定的函数来完成这些功能。这两种方法的确别在于：当测量值是准确的，没有误差时，一般用插值；当测量值与真实值有误差时，一般用数据拟合。1. 插值对于一维曲线的插值，一般用到的函数yi=interp1(X,Y,xi,method) ，其中method包括nearst，linear，spline，cubic。命令1 interp1功能 一维数据插值。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系如下图。图1 数据点与插值点关系示意图格式 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。yi = interp1(x,Y,xi,method) %用指定的算法计算插值：nearest：最近邻点插值，直接完成计算；linear：线性插值（缺省方式），直接完成计算；spline：三次样条函数插值；pchip：分段三次Hermite插值；cubic：与pchip操作相同；v5cubic：在MATLAB 5.0中的三次插值。对于超出x范围的xi的分量，使用方法nearest、linear、v5cubic的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。例 x = 0:10; y = x.*sin(x); xx = 0:.25:10; yy = interp1(x,y,xx); plot(x,y,'kd',xx,yy)例 year = 1900:10:2010;product = 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)对于二维曲面的插值，一般用到的函数zi=interp2(X,Y,Z,xi,yi,method)，其中method也和上面一样，常用的是cubic。命令2 interp2功能 二维数据内插值格式 ZI = interp2(X,Y,Z,XI,YI) %返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)Xi(i,j),yi(i,j)。ZI = interp2(Z,XI,YI) %缺省地，X=1:n、Y=1:m，其中m,n=size(Z)。再按第一种情形进行计算。ZI = interp2(Z,n) %作n次递归计算，在Z的每两个元素之间插入它们的二维插值，这样，Z的阶数将不断增加。interp2(Z)等价于interp2(z,1)。ZI = interp2(X,Y,Z,XI,YI,method) %用指定的算法method计算二维插值：linear：双线性插值算法（缺省算法）；nearest：最临近插值；spline：三次样条插值；cubic：双三次插值。例 X,Y =