Matlab.doc

1数值计算和符号计算功能例如，求解线性方程组：在MATLAB命令窗口输入命令：a=2,3,-1;3,-5,3;6,3,-8; b=7;8;9; x=inv(a)*b也可以通过符号计算来解此方程syms x1 x2 x3x1,x2,x3=solve(2*x1+3*x2-x3-7,3*x1-5*x2+3*x3-8,6*x1+3*x2-8*x3-9)2绘图功能例如，分别绘制函数y=300sinx/x和y=x2的曲线x=-20:0.1:20;plot(x,300*sin(x)./x,:,x,x.2);2设置搜索路径（1）用path命令设置搜索路径例如，将用户目录c:mydir加到搜索路径下：path(path,e:matlabwork)（2）用对话框设置搜索路径 在MATLAB的File菜单中选择Set Path命令或在命令窗口执行pathtool命令1变量命名在MATLAB 7.X中，变量名是以字母开头，后接字母、数字或下画线的字符序列，最多63个字符。【例1.1】当 时，计算表达式 的值，并将结果赋给变量y，然后显示出结果。在MATLAB命令窗口分别输入命令：x=sqrt(1+pi);y=(exp(x)+log(abs(sin(x)2-sin(x*x)/(x-5*i) y = 0.5690 + 1.3980i其中，pi和i都是MATLAB定义的变量，分别代表圆周率和虚数单位。【例1.2】利用M文件建立mydata矩阵。（1）启动有关文本编辑程序或MATLAB的M-file编辑器（见第2章），并输入待建矩阵：mydata=1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9; 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9; 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9;（2）把输入的内容存盘（设文件名为matfund.m）。（3）在MATLAB命令窗口中输入matfund，即运行该M文件，就会自动建立一个名为mydata的矩阵，可供以后使用。【例1.3】建立矩阵A，然后找出在10，20区间的元素的位置。A=4,15,-45,10,6;56,0,17,-45,0A = 4 15 -45 10 6 56 0 17 -45 0 find(A=10 & A=a&ch=z); %找小写字母的位置ch(k)=ch(k)-(a-A); %将小写字母变成相应的大写字母char(ch) length(k) %统计小写字母的个数ans = 14把字符串的内容作为对应的MATLAB语句来执行。 eval(t)其中t为字符串。例如：x=3;y=4;m=sqrt(x*x+y*y),x+y;y=eval(m)y = 5 7（1）若字符串中的字符含有单撇号，则该单撇号字符应用两个单撇号来表示。例如：disp(Its a book.)将输出：Its a book.（2）可以用字符串向量的形式连接多个字符串，即用中括号括起来。例如：x=8; y=10;disp(num2str(x), + ,num2str(y), = ,num2str(x+y)其中disp函数的自变量是一个长字符串。输出为8 + 10 = 18【例2.3】计算分段函数 ：程序如下：x=input(请输入x的值:);if x=A & c=a& c=0& c=9 disp(str2num(c)2);else disp(c);end【例2.5】将例2.4改用switch语句实现。c=input(请输入一个字符：,s);cc=abs(c);switch(cc) case num2cell(abs(A):abs(Z) disp(lower(c); case num2cell(abs(a):abs(z) disp(upper(c); case num2cell(abs(0):abs(9) disp(abs(c)-abs(0)2); otherwise disp(c);endnum2cell函数是将数值矩阵转化为单元矩阵，num2cell（1:5）等价于1,2,3,4,5。【例2.7】一个3位整数各位数字的立方和等于该数本身则称该数为水仙花数。输出全部水仙花数。程序如下：shu=; %用于存放结果，先赋空值for m=100:999 m1=fix(m/100); %求m的百位数字 m2=rem(fix(m/10),10); %求m的十位数字 m3=rem(m,10); %求m的个位数字 if m=m1*m1*m1+m2*m2*m2+m3*m3*m3 shu=shu,m; %存入结果 endendshu【例2.10】求使 最小的n。程序如下：y=0;n=0;while (y=1.5)n=n+1;y=y+1/n/n;enddisp(满足条件的n是：,num2str(n)【例2.19】计算 ，当n=100时，求y的值。用循环结构实现y=0; 采用向量求和的方法n=100; n=100;for i=1:n i=1:n;y=y+1/i/i; f=1./i2;end y=sum(f);disp(y) disp(y)【例3.1】绘制曲线t=0:0.1:8*pi;x=cos(t)+t.*sin(t);y=sin(t)-t.*cos(t);plot(x,y);【例3.2】用不同标度在同一坐标内绘制曲线y1=0.2e0.5xcos(4px)和y2=1.5e0.5x cos(px)。程序如下：x=0:pi/100:2*pi;y1=0.2*exp(-0.5*x).*cos(4*pi*x);y2=1.5*exp(-0.5*x).*cos(pi*x);plotyy(x,y1,x,y2);【例3.3】在同一坐标内，分别用不同线型和颜色绘制曲线y1=0.2e0.5xcos(4px)和y2=1.5e0.5x cos(px)。标记两曲线交叉点。程序如下：x=linspace(0,2*pi,1000);y1=0.2*exp(-0.5*x).*cos(4*pi*x);y2=1.5*exp(-0.5*x).*cos(pi*x);k=find(abs(y1-y2)1e-2); %查找y1与y2相等点(近似相等)的下标x1=x(k);%取y1与y2相等点的x坐标y3=0.2*exp(-0.5*x1).*cos(4*pi*x1); %求y1与y2值相等点的y坐标plot(x,y1,x,y2,k:,x1,y3,bp);【例3.4】在0x2p区间内，绘制曲线y1=e0.5x和y2=cos(4px)e0.5x，并添加图形标注。程序如下：x=0:pi/100:2*pi;y1=exp(-0.5*x);y2=exp(-0.5*x).*sin(2*x);plot(x,y1,x,y2)title(x from 0 to 2pi);%加图形标题xlabel(Variable X);%加X轴说明ylabel(Variable Y);%加Y轴说明text(1.5,0.5,曲线y1=e-0.5x);%在指定位置添加图形说明text(3,-0.1,曲线y2=cos(4pix)e-0.5x); legend(y1,y2)%加图例【例3.5】绘制曲线及其包络线。程序如下：t=(0:pi/100:pi);y1=sin(t)*1,-1;%包络线函数值y2=sin(t).*sin(9*t);plot(t,y1,y2)grid on;%加网格线box on;%加坐标边框axis equal%坐标轴采用等刻度【例3.7】在图形窗口中，以子图形式同时绘制多根曲线。程序如下：subplot(2,2,1 3);fplot(x-cos(x3)-sin(2*x2),-3,3); xlabel(a);x=-3:0.1:3;subplot(2,2,2);y2 = sin(2.*x.2); plot(x,y2);xlabel(b); axis(-3 3 -1.2 1.2)subplot(2,2,4);y3 = cos(x.3); plot(x,y3);xlabel(c); axis(-3 3 -1.2 1.2); grid on;【例3.8】绘制y=ex的对数坐标图并与直角线性坐标图进行比较。程序如下：x=0:0.1:10;y=exp(-x);subplot(2,2,1);plot(x,y);title(plot(x,y);grid on;subplot(2,2,2);semilogx(x,y);title(semilogx(x,y);grid on;subplot(2,2,3);semilogy(x,y);title(semilogy(x,y);grid on;subplot(2,2,4);loglog(x,y);title(loglog(x,y);grid on;【例3.9】已知t 0,6p，绘制阿基米德螺线r=a+bt图，并标记数据点。程序如下：t=0:pi/20:6*pi;a=2;b=3;r=a+b*t;polar(t,r,-*);x=51,82,34,47;67,78,68,90;78,85,65,50;subplot(1,2,1);bar(x,group);title(Group);axis(0 5 0 100);subplot(1,2,2);barh(x,stack);title(Stack);【例3.12】绘制一个蓝色的六边形。程序如下：n=6; dt=2*pi/n; st=0:dt:2*pi;t=st,st(1); %数据向量的首尾重合，使图形封闭。x=sin(t);y=cos(t);fill(x,y,b);axis(-1.5 1.5 -1.5 1.5)【例3.13】绘制三维曲线。程序如下：t=0:pi/10:10*pi;x=sin(t)+t.*cos(t);y=cos(t)-t.*sin(t);z=t;plot3(x,y,z);axis(-30 30 -30 30 0 35)title(Line in 3-D Space);xlabel(X);ylabel(Y);zlabel(Z);grid on;【例3.14】绘制三维曲面图z=sinx2+cosy2，x0，p，y0，p/2。程序如下：x,y=meshgrid(0:pi/100:pi, 0:pi/100:pi/2);z=sin(x.2)+cos(y.2);mesh(x,y,z);axis(0 4 0 1.8 -1.5 1.5);【例3.18】绘制多峰函数的瀑布图和等高线图。程序如下：subplot(1,2,1);X,Y,Z=peaks(30);waterfall(X,Y,Z)xlabel(X-axis),ylabel(Y-axis),zlabel(Z-axis);subplot(1,2,2);contour3(X,Y,Z,12,k); %其中12代表高度的等级数xlabel(X-axis),ylabel(Y-axis),zlabel(Z-axis);【例3.22】已知，绘制三维曲面图，并裁掉图中x和y都小于1.5的部分。程序如下：x,y=meshgrid(-5:0.2:5);z=sin(sqrt(x.2+y.2)./(sqrt(x.2+y.2)+eps);subplot(1,2,1);mesh(x,y,z);i=find(x-1.5 & y-1.5);z1=z;z1(i)=NaN;subplot(1,2,2);mesh(x,y,z1);【例3.23】隐函数绘图应用举例。程序如下：subplot(2,2,1);ezplot(cos(tan(pi*x),0,1)subplot(2,2,3);ezplot(x2-y4)subplot(2,2,2 4);ezplot(5*cos(5*t),4*sqrt(2*t) ,0,2*pi);2ezplot3函数ezplot3函数的调用格式为ezplot3(x,y,z,a,b)对于参数方程在区间a tb绘制x=funx(t)，y=funy(t)和z=funz(t)的图形。未指定区间时，默认在区间0t2p绘制图形。例如，绘制例3.13的三维曲线也可以采用以下命令：ezplot3(sin(t)+t*cos(t),cos(t)-t*sin(t),t,0,10*pi);【例3.25】绘制一个水平放置的瓶状柱面，并且将它绕z轴旋转。程序如下：t=0:pi/20:2*pi;x,y,z= cylinder(2+sin(t),30);mesh(z,y,x)axis off;shading interp;colormap(hsv);for k=1:20 view(-37.5+18*(k-1),30)%改变视点 M(k) =getframe;%将图形保存到m矩阵endmovie(M,2);%播放画面2次【例3.26】生成例3.13的图形轨迹。t=0:pi/250:10*pi;x=sin(t)+t.*cos(t);y=cos(t)-t.*sin(t);comet3(x,y,t);【例4.1】分别建立33、32和与矩阵A同样大小的零矩阵。（1）建立一个33的零矩阵。zeros(3)ans= 0 0 0 0 0 0 0 0 0（2）建立一个23的零矩阵。zeros(2,3)（3）设A为23矩阵，则可以用zeros(size(A)建立一个与矩阵A同样大小的零矩阵。A=1 2 3;4 5 6; %产生一个23阶矩阵Azeros(size(A) %产生一个与矩阵A同样大小的零矩阵【例4.2】建立随机矩阵：（1）在区间10,30内均匀分布的4阶随机矩阵。（2）均值为0.6、方差为0.1的4阶正态分布随机矩阵。产生(0,1)区间均匀分布随机矩阵使用rand函数，假设得到了一组满足(0,1)区间均匀分布的随机数xi，则若想得到在任意a,b区间上均匀分布的随机数，只需用yi=a+(ba)xi计算即可。产生均值为0、方差为1的标准正态分布随机矩阵使用randn函数，假设已经得到了一组标准正态分布随机数xi，如果想更一般地得到均值为、方差为2的随机数，可用yi=+xi计算出来。针对本例，命令如下：a=10;b=30;x=a+(b-a)*rand(4)y=0.6+sqrt(0.1)*randn(4)【例4.3】将101125等25个数填入一个5行5列的表格中，使其每行、每列及对角线的和均为565。一个5阶魔方矩阵的每行、每列及对角线的和均为65，对其每个元素都加100后这些和变为565。完成其功能的命令如下：M=100+magic(5)【例4.5】求(x+y)4的展开式。在MATLAB命令窗口，输入命令：pascal(5)ans= 1 1 1 1 1 1 2 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70矩阵次对角线上的元素1，4，6，4，1即为展开式的系数，即(x+y)4=x4+4x3y+6x2y2+4xy3+y4构造对角矩阵设V为具有m个元素的向量，diag(V,k)的功能是产生一个nn（n=m+|k|）对角阵，其第k条对角线的元素即为向量V的元素。 例如：diag(1:3,-1)ans = 0 0 0 0 1 0 0 0 0 2 0 0 0 0 3 0省略k时，相当于k为0，其主对角线元素即为向量V的元素。 【例4.6】先建立55矩阵A，然后将A的第一行元素乘以1，第二行乘以2，第五行乘以5。用一个对角矩阵左乘一个矩阵时，相当于用对角阵的第一个元素乘以该矩阵的第一行，用对角阵的第二个元素乘以该矩阵的第二行依此类推，因此，只需按要求构造一个对角矩阵D，并用D左乘A即可。命令如下：A=1:5;2:6;3:7;4:8;5:9D=diag(1:5);D*A %用D左乘A,对A的每行乘以一个指定常数【例4.7】用求特征值的方法解方程：3x57x4+5x2+2x18=0先构造与方程对应的多项式的伴随矩阵A，再求A的特征值。A的特征值即为方程的根。命令如下：p=3,-7,0,5,2,-18;A=compan(p); %构造伴随矩阵Ax1=eig(A) %求A的特征值x2=roots(p) %直接求多项式p的零点【例4.9】利用矩阵求逆方法解线性方程组A=1,-2,3;3,-1,5;2,1,5;b=1;2;3;x=inv(A)*b【例4.10】用左除运算符求解下列相同系数矩阵的两个线性代数方程组的解。解法1：分别解线性方程组。A=1,-1,1;5,-4,3;2,1,1; b1=2;-3;1;b2=3;4;-5;x=Ab1y=Ab2解法2：将两个线性方程组连在一起求解。A=1,-1,1;5,-4,3;2,1,1;b=2,3;-3,4;1,-5;xy=Ab利用第2种格式对矩阵A进行LU分解：L,U,P=lu(A)clearA=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;ticx2=Ab; %用左除运算求解tocticx1=inv(A)*b; %用求逆运算求解tocticL,U=lu(A); %LU分解x3=U(Lb); %用LU分解求解tocx1=x1x2=x2x3=x3其中tic和toc两个函数配合使用用于计算程序的执行时间，tic记录当前时间，toc记录或显示使用tic函数以来所花费的时间。运行结果说明，x1、x2、x3的值相同，通过LU分解求值所化运行时间最少。 【例5.2】已知 ，求矩阵x、y所有同一位置上的较大元素构成的新矩阵p。x=443,45,43;67,34,-43;y=65,73,34;61,84,326;p=max(x,y)a=7,4,9,8,10,7,8,7,8,7;7,6,10,5,9,8,10,9,5,6;mean(a)std(a)两人成绩的平均值相同，但小明的成绩的标准方差较小，说明小明的成绩波动较小，成绩更稳定。【例5.9】表5.3所示为我国06个月婴儿的体重、身长参考标准，用3次样条插值分别求得婴儿出生后半个月到5个半月每隔1个月的身长、体重参考值。表5.3我国婴儿体重、身长计量表出生1月2月3月4月5月6月身长（cm）50.656.559.662.364.665.968.1体重（kg）3.274.975.956.737.327.708.22设时间变量t为一行向量，c为一个两列矩阵，其中第1列存放身长，第2列储存体重。命令如下：h =0:1:6;t=50.6,56.5,59.6,62.3,64.6,65.9,68.1;3.27,4.97,5.95,6.73,7.32,7.70,8.22;XI =0.5:1:5.5XI = 0.5000 1.5000 2.5000 3.5000 4.5000 5.5000YI=interp1(h,t,XI,spline) %用3次样条插值计算【例5.11】已知数据x，y如表5.5所示，试求2次拟合多项式f(x)，然后求x=0.05，0.25，0.45，0.65，0.85，1.05各点的函数近似值。x=0:0.1:1;y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2;f=polyfit(x,y,2) %计算2次拟合多项式的系数f = -9.8108 20.1293 -0.0317故f(x)=-9.8108 x2+20.1293x-0.0317。可以用polyval函数按所得的多项式计算xi各点上多项式的值，polyval函数将在5.5.3节中详细介绍。下面先用polyval求得xi各点上的函数近似值：xi=0.05:0.2:1.05;yi=polyval(f,xi)plot(x,y,:o,xi,yi,-*)【例5.12】给定数学函数x(t)= 3cos(220t+/3)+5sin(250t)，取N=128，试对t从01s采样，用fft作快速傅立叶变换，绘制相应的振幅频率图。N=128; % 采样点数T=1; % 采样时间终点t=linspace(0,T,N); % 给出N个采样时间ti(I=1:N)x=3*cos(2*pi*20*t+pi/3)+5*sin(2*pi*50*t); % 求各采样点样本值xdt=t(2)-t(1); % 采样周期f=1/dt; % 采样频率(Hz)X=fft(x); % 计算x的快速傅立叶变换XF=X(1:N/2+1); % F(k)=X(k)(k=1:N/2+1)f=f*(0:N/2)/N; % 使频率轴f从零开始plot(f,abs(F),-*) % 绘制振幅频率图xlabel(Frequency);ylabel(|F(k)|)A=1,8,0,0,-10;B=2,-1,3;C=conv(A,B)Q,r=deconv(A,B)从上面的运行可知，两个多项式的乘积是一个6次多项式：2x6+ 15x5- 5x4+24x3-20x2+10x-30。多项式A除以多项式B的商式P为0.5x2+4.25x+1.375，余式r为-11.375x-14.125。【例6.1】设f(x)=sinx，用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f(x)的图像。为确定计算数值导数的点，假设在0，20区间内以为步长求数值导数。下面用3种方法求f(x)在这些点的导数。首先用一个5次多项式p(x)拟合函数f(x)，并对p(x)求一般意义下的导数dp(x)，求出dp(x)在假设点的值；第2种方法用diff函数直接求f(x)在假设点的数值导数；第3种方法先求出导函数f(x)=cosx，然后直接求f(x)在假设点的导数。x=0:pi/24:pi;%用5次多项式p拟合f(x),并对拟合多项式p求导数dp在假设点的函数值p=polyfit(x,sin(x),5);dp=polyder(p);dpx=polyval(dp,x);%直接对sin(x)求数值导数dx=diff(sin(x,pi+pi/24)/(pi/24);%求函数f的导函数g在假设点的导数gx=cos(x); %作图plot(x,dpx,x,dx,o,x,gx,+); 【例6.2】生成一个5阶魔方矩阵，按列进行差分运算。M=magic(5)M= 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9DM=diff(M) %计算V的一阶差分DM= 6 19 6 6 1 19 1 6 6 6 6 6 6 1 19 1 6 6 19 6可以看出，diff函数对矩阵的每一列都进行差分运算，因而结果矩阵的列数是不变的，只有行数减1。矩阵DM第3列值相同，表明原矩阵第3列是等间距的。【例6.3】求S=xsin(x)/(1+|cos(x)|)在0,pi的积分。（1）建立被积函数文件fe.m。function f=fe(x)f=x.*sin(x)./(1+abs(cos(x);（2）调用数值积分函数quad求定积分。quad(fe,0,pi)【例6.4】分别用quad函数和quadl函数求椭圆积分 的近似值，并在相同的积分精度下，比较函数的调用次数。调用函数quad求定积分：format long; fx=inline(1./sqrt(1+X.4); %定义一个语句函数I,n=quad(fx,0,1,1e-10)%注意函数名不加号调用函数quadl求定积分：I,n=quadl(fx,0,1,1e-10)【例6.5】求定积分（1）建立被积函数文件feln.m。function y=feln(x) y=exp(x).*log(x);（2）调用数值积分函数quadgk求定积分。format long;I=quadgk(feln,0,1)I= 1.317902162414081也可使用匿名函数方法求解：I=quadgk(x)(exp(x).*log(x),0,1);匿名函数是构造简单函数的一种方法，是函数指针，第1个括号里面是自变量，第2个括号里面是代表函数运算的表达式。【例6.7】计算二重定积分（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.5745Ki %Ki=1050 如果使用inline函数，则命令如下：f=inline(exp(-x.2/2).*sin(x.2+y),x,y);I=dblquad(f,2,2,1,1)（1）建立函数文件funst.m。function yp=funst(x,y)yp=sec(x)+y*tan(x);（2）求解微分方程。x0=0;xf=1;y0=pi/2;x,y=ode23(funst,x0,xf,y0);%求数值解y1=(x+pi/2)./cos(x); %求精确解x,y,y1（1）建立函数文件verderpol.m。function xprime=verderpol(t,x)global mu;xprime=x(2);mu*(1x(1)2)*x(2)x(1);（2）求解微分方程。global mu;mu=2;y0=1;0;t,x=ode45(verderpol,0,20,y0);（3）绘制解的曲线。 subplot(1,2,1);plot(t,x); %系统时间响应曲线 subplot(1,2,2);plot(x(:,1),x(:,2) %系统相平面曲线（1）建立模型的函数文件lorenz.m。function xdot=lorenz(t,x)xdot=8/3,0,x(2);0,10,10;x(2),28,1*x;（2）解微分方程组。t,x=ode23(lorenz,0,80,0;0;eps);（3）绘制系统相平面图plot3(x(:,1),x(:,2),x(:,3);axis(10,45,-15,20,-30,25);作业题（1）syms x y %符号运算求带变量上下限二重积分f=x*y;a=int(f,x,y.2,y+2);b=int(a,y,-1,2);eval(b)（2）x=0.3:0.2:1.5;y=0.3895,0.6598,0.9147,1.1611,1.3971,1.6212,1.8325;f=polyfit(x,y,2); %2次多项式拟合fx=inline(poly2sym(f); %转换成函数I=quad(fx,0.5,1.5) %求积分 第7章 符号计算【例1】syms x y;s1=x3-y3;factor(s1) %对s分解因式ans =(x - y)*(x2 + x*y + y2)s2=(-7*x2-8*y2)*(-x2+3*y2);expand(s2) %对s展开ans =7*x4 - 13*x2*y2 - 24*y4s3=(x+y)*(x2+y2+1)collect(s3,y) %对s按变量x合并同类项ans =y3 + x*y2 + (x2 + 1)*y + x*(x2 + 1)factor(sym(630) %对符号整数分解因式ans =2*32*5*7极限3：f=sym(x*(sqrt(x2+1)-x);limit(f,sym(x),inf,left)ans =1/2极限4：f= sym(cot(x) (1/log(x);limit(f,x,0,right)ans =1/exp(1)极限1：f=sym(sin(x+h)-sin(x)/h);limit(f,sym(h),0)ans =cos(x)极限2：f=sym(1+t/x) x);limit(f,inf)ans =exp(t)%求(1)f1=cos(x*x);diff(f1) %未指定求导变量和阶数，按默认规则处理ans =(-2)*x*sin(x2)diff(f1,x,2) %求f1对x的二阶导数ans =- 2*sin(x2) - 4*x2*cos(x2)diff(f1,x,3) %求f1对x的三阶导数ans =8*x3*sin(x2) - 12*x*cos(x2)%求(2)f21=a*(t-sin(t);f22=b*(1-cos(t);diff(f22)/diff(f21) %求y对x的一阶导数 ans =-(b*sin(t)/(a*(cos(t) - 1)%求(3)f3= x6-3*y4+2*x2*y2;diff(f3,x) %求