1、Matlab 的建模与仿真第三章：matlab 的数值计算例 3-1 创建矩阵举例x=1 3 5;2 4 6y=1,3,5;2,4,6z=x;7 8 9t=1:10例 3-2 矩阵元素的访问举例C=1 2 3 4 5;6 7 8 9 10;11 12 13 14 15;16 17 18 19 20C(2,3)C(:,3)C(3,:)C(1:3,3:4) 例 3-3 特殊矩阵生成示例i=3;j=4;zeros(i,j)ones(i,j)eye(j)rand(j)magic(i)例 3-4 矩阵的四则运算示例A=1 2 3;4 5 6;7 8 9B=9 8 7;6 5 4;3 2 1A+BB-AA
2、*B例 3-5 help matfun 示例help matfun例 3-6 help rank 示例help rank例 3-7 求矩阵 A 的有关特征函数，并计算矩阵 A 的特征值与奇异值A=1 2 3;4 5 6;7 8 9x,y=eig(A) %x 为特征向量，y 为特征值svd(A) %求 A 的奇异值例 3-8 数组的创建与访问示例A=9 8 7;6 5 4;3 2 1B=1 2 3;4 5 6;7 8 9A(2,3)B(:,3)B(3,:)A(1:2,2:3)例 3-9 数组内元素的更换示例A=9 8 7;6 5 4;3 2 1A(3,3)=10例 3-10 数组相加减示例A=9
3、 8 7;6 5 4;3 2 1B=1 2 3;4 5 6;7 8 9C1=A+BC2=A-B例 3-11 数组相乘除示例A=9 8 7;6 5 4;3 2 1B=1 2 3;4 5 6;7 8 9C1=A.*BC2=A./BC3=A.B例 3-12 如果用户希望了解基本数学函数的全部命令，可以先进入 C:Proam FileMATLAB704toolboxmatlabelfum,然后在 MATLAB 命令窗口中输入以下命令语句:type content例 3-13 求解齐次线性代数方程组 x350yz420xyz40xyzA=1 3 5;4 -1 2;1 1 -4;R=rank(A) %这里
4、 A 的秩等于 A 的列数，所以方程组有零解例 3-14 求解齐次方程组4x20318yzzA=4 4 2;-3 12 3;8 -2 1;2 12 4;R=rank(A) %这里 A 的秩小于 A 的列数 3，所以方程组有解，可用 null 求解x1=null(A)x2=null(sym(A)例 3-15 求解265438170xyz%首先确定系数矩阵和扩展矩阵的秩A=1 -2 6;3 8 1;18 -1 3;B=54 -6 70;RA=rank(A)RB=rank(A B) %此时，RA=RB=3,刚好与矩阵的列数相同det(A) %det(A)不等于零，则方程组有唯一解，可以通过左除法 A
5、B 求解方程组xx=AB 例 3-16 求解231456789xyz%首先确定系数矩阵和扩展矩阵的秩A=1 2 3;4 5 6;7 8 9B=-1 -2 -3RA=rank(A)RB=rank(A B) %由上式可知， RA=RB=2epsfor i=k+1:nm=(i,k)=A1(i,k)/A1(k,k);for j=k+1:nA1(i,j)=A1(I,j)-m(I,k)*A1(k,j);endB1(i)=B1(i)-m(i,k)*B1(k);end %loop ielsedisp(Nave Gaussian Elimination fails)returnend % ifend % loo
6、p k%获得结果%x(n)=B1(n)/A1(n,n);for i=n-1:-1:1s=B1(i);for j=i+1:ns=s-A1(i,j)*x(j);end %loop jx(i)=s/A1(i,i);endx例 3-18 求方阵 A=1 3 5;2 4 6;7 8 9的特征多项式、特征值和特征向量%输入方阵 AA=1 3 5;2 4 6;7 8 9;%f1 为方阵 A 特征多项式的系数 f1=poly(A)%f2 为特征多项式的表达式f2=poly2str(f1,x)f2=x3-14x2-40x-4.2437e-014%用 roots 求特征多项式 f1 的根x1=roots(f1)%
7、用 eig 求方阵 A 的特征向量x2=eig(A)%使用 eig 求方阵 A 的特征值和特征向量d x3=eig(A)例 3-19 用符号法求解二阶常微分方程 的通解及满足 =0.4, =0.7221dytt(0)y()的特解%用符号法求解时，首先需将常微分方程符号化，可用Dmy表示函数的 m 阶导数%小写字母t为系统默认自变量而且可以缺省，综上，例子即可写成 D2y+y=1-t2%求解常微分方程通解y=dsolve(D2y+y=1-t2)%求解常微分方程的特解y=dslove(D2y+y=1-t2,y(0)=0.4,Dy(0)=0.7)例 3-20 求解常微分方程 的通解23452dxyz
8、tdzxyzt%使用 dslove 函数求解常微分方程%设计程序求解x,y,z=dslove(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t);x=simple(x) %将 x 简化y=simple(y) %将 y 简化z=simple(z) %将 z 简化%运行结果如下例 3-21 用 2/3 阶龙格库塔法求解二阶常微分方程 ，满足 ，221dxytt(0)4y，t-10 10(0)7y%首先需要用 function 建立微分方程的 M 函数文件，然后以 ex320 存盘function dy=ex320(t,y) %定义函数文件名，输入输出变量
9、dy=zeros(2,1);dy(1)=dy(2);dy(2)=-y(1)+1-t2;%在命令窗口中输入 ode23 命令ode23ex320,-10 10,4 7),grid%则输出图形如图 3-1 所示例 3-22 已知微分方程组 ，绘制 y1,y2,y3 曲线123312130.8(),(),(0)yyy%首先需要用 function 建立微分方程的 M 函数文件，然后以 ex321 存盘function dy=ex321(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.8*y(1)*y(2);%在命令窗口中输入 o
10、de23 命令t,y=ode45(ex321,0 20,1 2 3);plot(t,y(:,1),-,t,y(:,2),*,t,y(:,3),+)legend(y1,y2,y3)例3-23依据泰勒展开法设计程序求解下面的联立常微分方程式 23()(,(0)1,01)4xtyttxyty%设计程序求解a=0;b=1;m=10;X0=1;0;n=length(X0);rs=zeros(n,m);d1=zeros(n,1);d2=d1;d3=d1;d4=d1;t=a;X=X0;h=(b-a)/m;for k=1:m%计算X,X”,X”and X”对t 的导数d1(1)=X(1)-X(2)+t*(2-
11、t*(1+t); % d1 is Xd1(2)=X(1)+X(2)+t2*(-4+t); d2(1)=d1(1)-d1(2)+2-t*(2+3*t); % d2 is X”d2(2)=d1(1)+d1(2)+t*(-8+3*t); d3(1)=d2(1)-d2(2)-2-6*t; %d3 isX” d3(2)=d2(1)+d2(2)-8+6*t; d4(1)=d3(1)-d3(2)-6; %d3 is X4 d4(2)=d3(1)+d3(2)+6; %计算X(t+h)for i=1:nX(i)=X(i)+h*(d1(i)+h*(d2(i)+h*(d3(i)+h*d4(i)/4)/3)/2);r
12、s(i,k)=X(i);end %for it=t+h;end %for kRs=X0 rs;%添增初值的数据Rs 第四章MATLAB的图形绘制例4-1 点的图形生成示例%在命令窗口中键入plot(4,5,bd) %在(4,5)坐标处家一个蓝色的菱形%程序运行结果显示如图4-1所示例4-2 线的图形生成示例%在命令窗口中键入x=1:10; %定义变量x的数值y=3*x; %定义函数yplot(x,y,b-*) %生成坐标点为* 形成的蓝色实线型曲线%程序运行结果显示如图4-2所示例4-3 单图单曲线图形生成示例%在命令窗口中键入t=0:2*pi/100:2*pi; %定义弧度行向量y=sin(
13、t); %计算正弦函数向量plot(t,y,g) %生成绿颜色的正弦曲线%程序运行结果显示如图4-3所示例4-4 单图双曲线图形生成示例%在命令窗口中键入t=0:2*pi/100:2*pi; %定义弧度行向量y1=sin(2*t); %计算正弦函数向量y2=cos(3*t); %计算余弦函数向量plot(t,y1,g,t,y2.r) %生成绿颜色的正弦和红颜色的余弦曲线%程序运行结果显示如图 4-4 所示亦或%在命令窗口中键入t=0:2*pi/100:2*pi; %定义弧度行向量y=sin(2*t);cos(3*t); %计算正弦和余弦函数向量plot(t,y) %生成绿颜色的正弦和红颜色的余
14、弦曲线%程序运行结果显示如图 4-5 所示例 4-5 多窗口曲线图形生成示例%在命令窗口中键入t=0:2*pi/100:2*pi; %定义弧度行向量figure(1) %创建序号为 1 的绘图窗口y1=sin(2*t); %计算正弦函数向量plot(t,y1,g) %在1号窗口中生成绿颜色的正弦曲线figure(2) %创建序号为 2 的绘图窗口y2=cos(3*t); %计算余弦函数向量plot(t,y2,g) %在2号窗口中生成绿颜色的正弦曲线例4-6 单窗口子图曲线图形生成示例%在命令窗口键入t=linspace(0,4*pi,100); %定义弧度行向量y1=sin(t); %计算正弦
15、函数向量y2=cos(t); %计算余弦函数向量y3=sin(t)./(cos(t)+eps); %计算正切函数向量y4=cos(t)./(sin(t)+eps); %计算正切函数向量subplot(2,2,1) %画左上角的图形plot(t,y1)subplot(2,2,2) %画右上角的图形plot(t,y2)subplot(2,2,3) %画左下角的图形plot(t,y3)subplot(2,2,4) %画右下角的图形plot(t,y4)%程序运行结果显示如图 4-7 所示例4-7 设置坐标轴名称及图形标题示例%在命令窗口键入t=0:2*pi/100:2*pi; %定义弧度行向量y1=s
16、in(2*t); %计算正弦函数向量y2=cos(3*t); %计算余弦函数向量plot(t,y1,g-,t,y2.r-) %生成绿颜色的正弦和红颜色的余弦线xlabel(x=0to2pi,FontSize,12) %标注 x 轴名称并设置名称为 12 号字体ylabel(幅值,FontSize,12) %标注y轴名称并设置名称为12号字体title(例4-7 图形输出) %定义全图名称%程序运行结果显示如图 4-8 示例4-8 对曲线进行文本注释示例%在命令窗口键入t=0:2*pi/100:2*pi; %定义弧度行向量y1=sin(2*t); %计算正弦函数向量y2=cos(3*t); %计
17、算余弦函数向量plot(t,y1,g-,t,y2,r-.) %生成绿颜色的正弦和红颜色的余弦线text(pi/2,sin(2*pi/10),y1leftarrow,FontSize,12);text(pi,cos(3*pi/10),y2rightarrow,FontSize,12);gtext(单击鼠标,FontSize,12); %在所需注释的位置单击鼠标放入文本字符串例4-9调整坐标轴和标定图例示例%在命令窗口中键入x=linspace(0,3*pi,60);y1=sin(x);y2=cos(x);subplot(221)plot(x,y1,x,y2)axis off; %取消横纵轴axi
18、s(square,equal); %设置横纵轴比例legend(y1,y2); %标定y1 和y2的图例grid on; %打开网格线subplot(222)plot(x,y1,x,y2)axis off; %取消横纵轴axis(xy,normal); %以预设置画横纵轴legend(y1,y2); %标定y1 和y2的图例grid on; %打开网格线subplot(223)plot(x,y1,x,y2)axis on; %恢复横纵轴axis(square,equal); %以预设置画横纵轴subplot(224)plot(x,y1,x,y2)axis on; %恢复横纵轴axis(xy,n
19、ormal); %设置横纵轴比例例4-10 调整坐标轴和标定图例示例%在命令窗口中键入t=0:0.05*pi:2*pi; %定义弧度行向量y1=sin(2*t); %计算正弦函数向量y2=cos(3*t); %计算余弦函数向量plot(t,y1,g-,t,y2,r-.) %生成绿颜色的正弦和红颜色的余弦线text(pi/2,sin(2*pi/10),y1leftarrow,FontSize,12);text(pi,cos(3*pi/10),y2rightarrow,FontSize,12);hold on %打开图形保持，若关闭图形保持，则为hold offy3=sin(2*t).*cos(3
20、*t);plot(t,y3,b-.)legend(y1,y2,y3)zoom on %启动图形缩放，通过单击鼠标左右键调整大小%程序运行结果显示如图 4-11 示例4-11复数向量绘图示例%在命令窗口中键入t=0:pi/90:4*pi;x=t.*exp(i*t); %t点乘负指数函数，形成等距螺旋线plot(x) %绘制复数向量图axis(image) %修饰图形，是得曲线居中%程序运行结果显示如图4-12所示例4-12分别绘制y=|500cos(3x)|+2的单对数和双对数坐标图%在命令窗口中键入x=0:0.02:2*pi;y=abs(500*cos(3*x)+2；figure(1)semi
21、logx(x,y); %单对数对x 轴绘图figure(2)semilogy(x,y); %单对数对y轴绘图figure(3)loglog(x,y); %双对数坐标绘图%程序运行结果显示如图4-13所示例4-13绘制cos(x)*sin(x) 的极坐标图%在命令窗口中键入x=0:0.02:2*pi;y=cos(x).*sin(x);polar(x,y); %绘制极坐标title(绘制极坐标);%程序运行结果显示如图4-14所示例4-14求方程 的根1505x%在命令窗口中键入x=1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1;y=roots(x); %求解方程的根r=abs(y
22、);t=angle(y);polar(t,r,bd); %将方程的根绘制出来title(方程根的分布图);%程序运行结果显示如图4-15所示例4-15绘制直方图和火柴杆图%在命令窗口中键入t=0:0.1:10;y=1-exp(-t).*cos(2*t);figure(1)bar(t,y) %绘制直方图title(直方图);figure(2)stem(t,y)title(火柴杆图)；%程序运行结果显示如图4-16所示例4-16绘制阶梯图，饼图，彗星曲线图，填充区域图，复数矢量图和极坐标累计图%在命令窗口中键入t=-20:2:20;y1=cos(t).*exp(-t/4);y2=exp(-1/3*
23、t).*cos(t);y3=sin(t);y4=cos(t);figure(1);stairs(t,y1); %绘制阶梯图title(阶梯图);figure(2);pie(t); %绘制饼图title(饼图);figure(3)aera(t,y3); %绘制填充区域图title(填充区域图);figure(4);compass(t,y3); %绘制复数矢量图title(复数矢量图);figure(5);rose(t,y4); %绘制极坐标累计图title(极坐标累计图);%程序运行结果显示如图4-17所示例4-17绘制三维螺旋线t=0:pi/100:5*pi;plot3(sin(3*t),co
24、s(3*t),t)grid%程序运行结果显示如图4-17所示例4-18绘制两条空间折线，第一条线上的空间坐标为 (1,4,7),(2,10,1),(10,6,9),第二条线上的空间节点坐标为(4,10,1),(5,1,2),(3,1,3)%将空间坐标点依次按照X,Y,Z 建立三个矩阵%在命令窗口中键入X=1 4;2 5;10 3;Y=4 10;10 1;6 1;Z=7 1;1 2;9 3;plot3(X Y Z)gridbox %添加三维箱体方框%程序运行结果显示如图4-19所示例4-19绘制三维随机网线图%在命令窗口中键入z=rand(10);mesh(z)%程序运行结果显示如图4-19所示
25、例4-20绘制例4-19的三维网线%在命令窗口中键入X=1 4;2 5;10 3;Y=4 10;10 1;6 1;Z=7 1;1 2;9 3;mesh(X,Y,Z)grid onbox %添加三维箱体框图colormap(0 0 1); %设置网线颜色，0 0 1为蓝色%程序运行结果显示如图4-21所示例4-21绘制4-20的三维曲面%在命令窗口中键入X=1 4;2 5;10 3;Y=4 10;10 1;6 1;Z=7 1;1 2;9 3;surf(X,Y,Z)grid onbox %添加三维箱体框图colormap(0 1 0 ); %设置网线颜色，0 1 0为绿色%程序运行结果显示如图4-
26、22所示例4-22设置图形颜色示例%在命令窗口中键入peak(35) %生成三维示范曲面图形colormap(1 0 0); %生成红色图形colormap(hot); %生成暖色色图%程序运行结果显示如图4-23所示例4-23设置图形视觉角度示例%在命令窗口中键入t=0:0.01:0.95;t1=0:0.01:0.95;t2=1:0.05:3;T=t1,-t2+3;x,y,z=cylinder(T,40);mesh(x,y,z);subplot(221);mesh(x,y,z);view(0,0);subplot(222);mesh(x,y,z);view(20,20);subplot(22
27、3);mesh(x,y,z);view(40,40);subplot(224);mesh(x,y,z);view(80,80);%程序运行结果显示如图4-24所示例4-24特殊三维绘图示例%在命令窗口中键入x,y,z=sphere(20);k=abs(z);subplot(221);bar3(3);title(三维直方图)subplot(222);pie3(30);title(三维饼图)subplot(223);stem3(5);title(三维离散杆图)subplot(224);surf(x,y,z,k);title(球面图)%程序运行结果显示如图4-25所示例4-25画一条直线，并将线的颜
28、色改为绿色，并且获得当前窗口属性%在命令窗口中键入subplot(121)L=line(0:5,0:5); %创建一个线对象subplot(122)L=line(0:5,0:5); set(L,Color,0 1 0) %设置线对象颜色为绿色%程序运行结果显示如图4-26所示%获得当前窗口图形属性get(gcf)例4-26句柄操作示例%在命令窗口中键入t=0:pi/270:2*pi;y1=sin(2*t);y2=cos(2*t);plot(t,y1,t,y2)%程序运行结果显示如图4-27所示%获得所有图形对象的句柄值h=findobj第五章：MATLAB程序设计及其仿真例5-1 已知一元二次
29、方程 ,试编写程序求解2y=x51023(1)()yy%首先建立一个M 函数：m5_1.mfunction y=m5_1(x)y=2*x2+5*x+10%在命令窗口输入变量1，2，3以及y(1)+y2(2)+y3(3)x1=1;x2=2;x3=3;m5_1(x1)+m5_1(x2)*m5_1(x2)+m5_1(x3)*m5_1(x3)*m5_1(x3)%运行结果如下例5-2 使用条件语句if编写程序示例%首先建立一个M 函数：m5_2.mfunction y=m5_2(x)if x=0y=x2+4*x+4elsey=x2-4*x+4end%在命令窗口输入变量x1=5x1=5;y1=m5_2(x
30、1);%运行结果如下%在命令窗口输入变量x2=-5x2=-5;y2=m5_2(x2);%运行结果如下例5-3 使用条件语句switch编写程序示例y=Acase x=90sum=m5_1(n);sum%运行结果如下例5-5 使用循环结构语句for编写程序示例%用for语句求1到任意自然函数的和%首先建立一个M 的函数：m5_5.mfunction m,sum=m5_5(n)sum=0;m=1;for m=1:nsum=sum+m;m=m+1;endm=m-1;%在命令窗口输入变量n=600n=600;n,sum=m5_5(n)%运行结果如下例5-6 使用 for编写嵌套循环程序示例%首先建立一
31、个M 的函数：m5_6.mm=1:2:20;n=length(m);for i=1:nfor j=1:nf(i)=cos(m(i)*pi/20)*sin(m(i)*pi/20);m(j)=m(j)*4;endendnmf%运行结果如下例5-7 求解所有3位整数中每位数字的立方和等于该数本身的数%本程序关键是如何确定3位整数的百位数，十位数和个位数%在MATLAB里可以通过使用函数fix(n/100)求3位整数n的百位数%使用函数rem(fix(n/10),10)求3位整数n的十位数%使用函数rem(n,10)求3位整数n 的个位数%此时建立一个M 文本文件：m5_7.mfor n=100:1:
32、999n1=fix(n/100);n2=rem(fix(n/10),10);n3=rem(n,10);if n=n1*n1*n1+n2*n2*n2+n3*n3*n3disp(n)endend%运行结果如下例5-8 随意输入一个字符串，字符串中的字符若为数字字符，则输出所对应的数值；若为小写字母，则输出所对应的大写字母；若为大写字母，则输出对应的小写字母；否则，则原样输出%首先建立一个M 文本文件：m5_8.mzf=input(pleas input a random string:,s);n=length(zf);for i=1:nif zf(i)=0elsedisp(zf(i);endend
33、zf%运行结果如下例5-9 随机输入3个自然数，然后按照从大到小的顺序排列%首先建立一个M 文本文件：m5_9.m%输入三个自然数clear all;x1=input(x1=);x2=input(x2=);x3=input(x3=);x=x1;x2;x3;s=length(x);%定义一个1行3 列的零矩阵sx=zeros(1,3);for i=1:sif x(i)=max(x) %求出三个自然数中最大的数sx(1)=x(i);elseif x(i)=min(x) %求出三个自然数中最小的数sx(3)=x(i);elsesx(2)=x(i);endendsx%运行结果如下例5-10猩猩吃香蕉的
34、问题：有一堆不知数目的香蕉，猩猩第一天吃掉一半，觉得没吃够，又多吃了一个，第二天依旧如此吃香蕉，即吃掉剩下香蕉的一半再加一个，以后天天如此，直至第十天早上发现只剩下一个香蕉了。问这堆香蕉原来的数量和每天香蕉的数量？%此题初看起来感觉无从下手，其实这是一个典型的递推。即先假设第1天共有X1个香蕉，第2天有x2个香蕉，,第9天有x9个，第10天有x10 个。从题干中可以知道x10=1,而且可以看出从x1x10之间存在xi=2*(xi+1+1)的关系，其中i=9,8,7,1%首先建立一个M 文本文件：m5_10.m%定义一个1x10 的零矩阵x=zeros(1,10);%给定初始值x(10)=1;i
35、=9;for n=i:-1:1x(n)=2*(x(n+1)+1);Total=x(n);end%输出香蕉总数量Total%输出每天香蕉的数量x%运行结果如下例5-11编写程序判断某一年是否为闰年众所周知，判断闰年的条件是能被4整除不能被100整除或者能被400整除的年份，那么接下来就以此为条件编写程序%首先建立一个M 文本文件：m5_11.m%输入四位数的某一年year=input(please input year:);sg=0;if rem(year,4)=0sg=sg+1;endif rem(year,100)=0sg=sg-1;endif rem(year,400)=0sg=sg+1;
36、endif sg=1fprintf(这是闰年 n,year)elsefprintf(这不是闰年 n,year)%运行结果如下%此时，在M文本编辑窗口中多次运行，每次运行时在MATLAB 命令窗口中分别输入1000,2000,2008,和2010%运行结果如下例5-12某商城对商品实行消费折扣促销，具体标准如表 5-2所示消费额度 折扣若总消费小于300元 无折扣若总消费大于等于300元并小于500元 1%折扣若总消费大于等于500元并小于800元 3%折扣若总消费大于等于800元并小于1500元 5%折扣若总消费大于等于1500元并小于3000元 7%折扣若总消费大于等于3000元并小于100
37、00元 9%折扣若总消费大于等于10000元 12%折扣表5-2促销折扣%首先建立一个M 文本文件：m5_12.mxf=input(:);switch fix(xf/100)case0,1,2zk=0;case3,4zk=1/100;case5,6,7zk=3/100;case8:14zk=5/100;case15:29zk=7/100;case30:99zk=9/100;otherwisezk=12/100;end%折扣后的实际消费额度xf=xf*(1-zk)%运行结果如下例 5-13 设计一个程序管理 10 个学生的成绩单%首先建立一个M 文本文件：m5_13.m%首先划分等级:F(60分
38、以下 )，D(60-69),C(70-79),B(80-89),A(90-99),A+(100)for i=1:9ai=60+i;bi=70+i;ci=80+i;di=90+i;end%输入学生名单以及对应的成绩name=甲,乙,丙,丁 ,戊, 己, 庚 ,辛,任, 癸;grade=50,77,88,66,99,100,70,90,82,55;%创建一个10元素的元胞数组data=cell(1,10);result=struct(name,name,grade,grade,data,data);%使用switch语句求出学生分数所对应的等级for j=1:10switch result(j).
39、gradecase aresult(j).data=D;case bresult(j).data=C;case cresult(j).data=B;case dresult(j).data=A;case 100result(j).data=A+;otherwiseresult(j).data=F;endend%打印学生姓名，分数以及等级信息disp(姓名,blanks(6), 分数 ,blanks(6),等级 );disp()for k=1:10disp(result(k).name,blanks(8),num2str(result(k).grade),blanks(8),result(k).data);end%运行结果如下
