MATLAB实验指导书加程序+上机实例

1、MATLAB 语言实验指导书华东交通大学电气学院张永贤2006 年 2 月实验一 MATLAB工作环境熟悉及简单命令的执行一、实验目的：熟悉MATLAB勺工作环境，学会使用 MATLAB®行一些简单的运算。二、实验内容：MATLAB勺启动和退出，熟悉MATLAB勺桌面(Desktop ),包括菜单(Menu)、 工具条 (Toolbar )、命令窗口(CommandWindow)、历史命令窗口、工作空间(Workspace) 等；完成一些基本的矩阵操作；学习使用在线帮助系统。三、实验步骤：1、启动 MATLAB熟悉 MATLAB勺桌面。2、在命令窗口执行命令完成以下运算，观察 wor

2、kspace的变化，记录运算结果。(1) (365-52 2-70 )3(2) >>area=pi*A2(3)已知 x=3, y=4,在 MATLAB求 z：2 3z xy z2x y(4)将下面的矩阵赋值给变量ml,在workspace中察看ml在内存中占用的字节数。16 23135 11 10 8m1 =976124 14 15 1执行以下命令>>m1( 2,3 )>>m1( 11 )>>m1( :, 3 )>>m1( 2 : 3,1 : 3 )>>m1( 1 ,4 ) + m1( 2 ,3 ) + m1( 3 ,2

3、) + m1( 4 ,1)(5)执行命令>>help abs查看函数abs的用法及用途，计算 abs( 3 + 4i )(6)执行命令>>x=0:6*pi;>>y=5*sin(x);>>plot(x,y)(6)运行MATLAB勺演示程序，>>demo以便对 MATLABT一个总体了解。四、思考题1、以下变量名是否合法为什么(1) x2(2) 3col(3) _row(4) for2、求以下变量的值，并在 MATLA呻验证。(1) a = 1 : 2 : 5 ;(2) b = a' a' a'(3) c = a

4、+ b ( 2 ,:)实验二MATLAB语言矩阵运算掌握基本的矩阵运及常用的函数。、实验内容：123 a4562411351147c 0 d 85223601、下列运是否合法，为什么如合法，结果是多少(1) result1 = a'(2) result2 = a * b(3) result3 = a + b(4) result4 = b * d(5) result5 = b ; c' * d(6) result6 = a . * b(7) result7 = a . / b(8) result8 = a . * c(9) result9 = a . b(10) resultIO

5、 = a . A2(11) resultll = a 人2(12) result11 = 2 . A a2、用MATLABt下面的的方程组。7212x1491532x27(1)22 115x3113213x40 xyz1 x 2y z w 82x y 3w 33x 3y 5z 6w 572123、已知A9153222 115132 13求矩阵A的秩(rank)求矩阵A的行列式(determinant)(3)求矩阵A的逆(inverse)(4) 求矩阵 A 的特征值及特征向量(eigenvalue and eigenvector) 4、关系运与逻辑运已知 a=20,b=-2,c=0,d=1(1)

6、 r1 = a > b(2) r2 = a > b & c > d(3) r3 = a = b* (-10)(4) r4 = b | c三、思考题10Qni 109 999q10y 22222,求 y=(用 format long 查看 y 的值)n 10实验三程序的编辑及调试一、实验目的：掌握MATLA翼序编辑、运行及调试方法。二、实验内容：1、启动 MATLAB后，点击File|New|M-File ,启动 MATLAB的程序 编辑及 调试器(Editor/Debugger )，编辑以下程序，点击File|Save 保存程序，注意文件名最好用英 文字符。点击 De

7、bug|Run运行程序，在命令窗口查看运行结果，程序如有错误则改正。注：数论中一个有趣的题目：任意一个正整数，若为偶数，则用2除之，若为奇数，则与3相乘再加上1。重复此过程，最终得到的结果为1。如：2 1310 516 842163 105 168421运行下面的程序，按程序提示输入n=1,2,3,5,7等数来验证这一结论。%classic "3n+1" problem from number theory.while 1n=input( 'Enter n,negative quits:' );if n<=0breakenda=n;while n>

8、;1if rem(n,2)=0n=n/2;elsen=3*n+1;enda=a,n;endaendm2、编程求满足i 1210000的最小mt。三、思考题用对分法求解方程 2e xsinx在0, 1内的解，并验证，在程序中统计出对分次数。提示：先将原方程转化成f (x) 2e x sin x 0的形式。对分法的基本思想是：一个一元方程f(x)=0，若f(x1)*f(x2)<0 ，则在x1,x2区间内有实数解。取该区间的中点 xm=(x1+x2)/2，判定f(x1)和f(x2)二者中哪一个与f(xm)异号,若f(x1)*f(xm)<0,则解存在的区间缩小为x1,xm,否则解存在的区间

9、缩小为 xm,x2。重复这样的步骤，直到区间的长度小于一个可以接受的小数(比如 原方程的解。1e-10),则认为中点即是实验四 函数的编写及调试一、实验目的：掌握MATLAESi数的编写及调试方法。二、实验内容：1、编写一个函数，计算下面函数的值，给出标量x的值，调用该函数后，返回 y的值。function y=myfun1(x)选择一些数据测试你编写的函数。2、编写一个函数求向量 x中元素的平均值、最大值、最小值、均方根值。function m_x,max_x,min_x,rms_x=myfun2(x)方均根值(RootMean Square)的计算公式为：rmsNx21用下面数据测试你写的

10、函数：(1) x=sin(0:6*pi)(2) x=rand(1,200),得到的 x为 200 个(0, 1)之间均匀分布的随机数。3、编写一个函数，给出一个向量 x%,function v=myvander(x)x2 ,xn ,生成如111x1x2xn222x1x2xnn 1n 1n，x1x2xn卜范德蒙矩阵。1例如：>>v=myvander(2 3 4 得v=5)16252764125生成一些数据测试你写的函数。三、思考题编写程序,用如下迭代公式求a的值分别为:3,17, 113。迭代的终止条件为10 5,迭代初值x01.0 ,迭代次数不超过100次。分别对迭代结果和222x

11、 a 2x x aa2xn 1实验五MATLAB勺绘图1、在同一坐标系下绘制下面三个函数在t0,4的图象。y 1 ty 2，ty 34 e 0.1t sin( t)2、编写程序，选择合适的步距，绘制下面函数在区间-6 , 6中的图象。3、用compass函数画下面相量图ua =1; ub =cos(-2*pi/3)+sin(-2*pi/3)*iuc=cos(2*pi/3)+sin(2*pi/3)*i;compass(ua,ub,uc,ua-ub,ub-uc,uc-ua)4、三维空间曲线绘制 z=0:4*pi;x=cos(z);y=sin(z);plot3(x,y,z)5、用mesh或surf函

12、数，绘制下面方程所表示的三维空间曲面，x和y的取值范围设为-3 ,3。22z 王L10 10三、思考题在同一坐标系下，用不同颜色和线型绘制以下两个函数在t -2, 2 范围内的图象。051t0.2120第y2 2e实验六MATLA吸值运算一、实验目的：掌握MATLAB1的数值运算函数。二、实验内容：1、求代数方程3x5 4x4 7x3 2x2 9x 12 0的5个根，并将其用星号(*)标记在复平面图上。(用roots和plot函数)。52、求代数方程x 1 0的5个根，并将其用星号(*)标记在复平面图上。(用roots和 plot函数)。3、求下面函数在,4区间内的过零点。(用fzero函)3

13、2 . ,、/、1f (x) x 2x sin(x) 5xcos(x) 一4、50已知R=50欧姆,OU=4V二极管D正向电流与电压的关系为:U dqKFIdIse其中：Ud为二极管正向电压Is为反向饱合电流，取 10-12AK为玻尔茨曼常数，*10-23T为绝对温度，取 300开尔文(27摄氏度)q为电子电荷*10-19C求此电路中的电流Id和二极管正向电压 Ud (要求用fsolve 函数求解)5、实验数据处理：已知某压力传感器的测试数据如下表pu101113141718222429343932p为压力值，u为电压值，试用多项式u(p) ap bp cp d来拟合其特性函数，求出a,b,c

14、,d ,并把拟合曲线和各个测试数据点画在同一幅图上。实验七MATLAB1、以原点为奇对称中心的方波y(wt),可以用相应频率的基波及其奇次谐波合成。，、4 .， y(wt) - sin wt1 . C ,1 . L , sin3wt sin5wt 351(2n 1)sin(2n1)wtn 1,2,3,取的阶数越多，越接近方波，但总消除不了边缘上的尖峰，这称为吉布斯效应。设方波频率为50Hz,时间 t 取0秒(f=50Hz,w=2*pi*f,h=1e-5,tf=40e-3,t= 0:h:tf),编写程序，画出如下用1次谐波、1,3次谐波、1,3,5,7,9次谐波，1,3,5，,19次谐波合成的近

15、似方波。(产生方波的函数为：square)2、用Simulink求解下图所示电路0100微秒内的响应。已知R=6*104欧，C=1700 微法，L=6*10-9 享，uc(0)=15kV。模块参数设置：Integrator 的 Initial condition:15kV在命令窗口为R,L,C赋值。仿真参数设置如下：Start time:0Stop time:100e-6Solver Type:Variable-stepSolver:ode45Max step size:1e-7Min step size:autoInitial step size:autoRelative tolerance

16、:1e-3Absolute tolerance:1e-6MATLA戚验程序实验1第1题.x=2*sin(85*pi/180)/(1+exp(2)x =(2) .x=2 1+2i; 5;y=*log(x+sqrt(1+xA(2)y =- +-或x=2 1+2i; 5;d=*log(x+sqrt(1+x*x) d =- +-或x=2 1+2*i; 5;y=*log(x+sqrt(1+xA(2)y =- +-(3) .a=:;g=(exp*a)-exp*a).*sin(a+/2+log(+a)/2)结果略(4)> > t=0:;> > f1=t.A2;> > f2

17、=t.A2-1;> > f3=t.A2-2*t+1;> > z=(t>=0&t<1).*f1+(t>=1&t<2).*f2+(t>=2&t<3).*f3 z =00第 2题(1)> > A=12 34 -4;34 7 87;3 65 7;> > B=1 3 -1;2 0 3;3 -2 7;> > A+6*Bans =18 52 -1046710521 53 49> > A-B+eye(3)ans =12 32-2338851 681( 2)>> A*B

18、ans =68 4462309 -72596154 -5241>> A.*Bans =12 10246802619 -13049(3)>> AA3ans =372262338244860424737014918860076678688454142118820>> A.A3ans =172839304-643930434365850327274625343(4)>> A/Bans =>> BA ans =(5)>> A,Bans = 12 34334765-487730-2-137>> A(1,3,:);BA2a

19、ns = 12 3 4 11 20346550-5-4711940第 3题> > A=1 2 3 4 5;6 7 8 9 10;11 12 13 14 15;16 17 18 19 20;21 22 23 24 25;> > B=3 0 16;17 -6 9;0 23 -4;9 7 0;4 13 11;> > C=A*B9315077258335237423520397588705557753890717>> D=C(3:5,2:3)D = 520 705 890第 4题(1)397557717>> a=100:999;>>

20、; k=find(rem(a,21)=0);素位置%找出能杯21 整除的元素位置,find() 函数找出不为0 的元>> x=length(k)%获得向量k 的元素个数并赋值给变量xx =43>> k=find(rem(a,21)=0) % 显示能杯21 整除的元素位置k =Columns 1 through 24627927 48 69 90 111 132 153 174 195 216 237 258300 321 342 363 384 405 426 447 468 489Columns 25 through 43510 531 552 573 594 615

21、 636783 804 825 846 867 888657 678 699720741762>> y=100+k-1%显示能杯21整除的元素y =Columns 1 through 23105 126 147 168 189 210 231252 273 294315336357378 399 420 441 462 483 504 525 546 567Columns 24 through 43588 609 630 651 672 693 714 735 756 777798819840>> sh(k)=;>> x=sh(1:end)x =e345fg

22、69%删除大写字母%显示处理后的字符实验2第 1题a=1 2 3;4 5 6;b=2 4 -1;1 3 5;c=1;0;-2;d=1 4 7;8 5 2;3 6 0;>> result1=a'result1 =1 42 53 6>> result2=a*b>> result3=a+b result3 =36258 11>> result4=b*d result4 =31 22 2240 49 13>> result5=b;c'*d% a 的转置%error 应采用点乘%求两个矩阵的和%矩阵相乘result5 = 31

23、22 2240 49 13-5 -87>> result6=a.*b result6 =28 -34 15 30%矩阵点乘>> result7=a./b result7 =%矩阵右点除861 882 903 924 945 966 987(2)sh='CDe345Efg69K'>> k=find(sh>='A'&sh<='Z');% 找出大写字母的位置>> result8=a.*c>> result9=a.b result9 =%error%矩阵左点除a和c维数不同&

24、gt;> result10=a.A2 result10=1 4916 25 36>> result11=aA2>> result12=2.Aa result12 =2 4816 32 64第 2题(1)%error等价于 a*a» A=7 2 1 -2;9 15 3-2;-2 -2 11 5;1 3 2 13;» B=4;7;-1;0;A/B 等价于 A*inv(B)» X=inv(A)*B%人旧等价于 inv(A)*B,X =» X1=ABX1 =(2)» a=1 1 1 0;1 2 1 -1;2 -1 0 -3

25、;3 3 5-6;>> b=1;8;3;5;» x=inv(a)*b x =第3题» A=7 2 1 -2;9 15 3-2;-2 -2 11 5;1 3 2 13;» a1=rank(A)a1 =4» a2=det(A)a2 =12568» a3=inv(A)a3 =%“向量A的特征向量，D为特征值>> V,D=eig(A) V =+ -+ -D =00+00000第4题» a=20;» b=-2;» c=0;» d=1;» r1 =a>br1 =1»

26、r2=a>b&c>dr2 =» y=myfun1(4)0» r3 = a = b* (-10) r3 =1» r4 = -b | cr4 =05» y=0;for k=-10:10y=y+pow2(k);endformat long»yy = +003实验3方法一a=0;for i=1:20a=a+pow2(i);ifa>10000m=i;breakendend方法二a=0;i=1;while (a<10000)a=a+pow2(i);i=i+1;endm=i-1;m实验4第1题function y=myfun1

27、(x)if x<=0y=sin(x);elseif x>0&x<=3y=x;elseif x>3y=-x+6;end运行结果：» y=myfun1(-pi/2)y =-1» y=myfun1(0)y =0» y=myfun1(2)y =2y =2第 2题function m_x,max_x,min_x,rms_x=myfun2(x)%求平均值sum_x=sum(x);%向量元素求和m,n=size(x);%最好用n=length(x);m_x=sum_x/n;%求最大值采用逐个比较方式if x(1)>x(2)max_x=x(1

28、);elsemax_x=x(2);endfor k=3:nif max_x<x(k)max_x=x(k);elsemax_x=max_x;%可省略endend%求最小值if x(1)<x(2)min_x=x(1);elsemin_x=x(2);endfor k=3:nif min_x>x(k)min_x=x(k);elsemin_x=min_x;%可省略endend%求均方根值sum_x2=0;for k=1:nsum_x2=sum_x2+x(k2;rms_x=sqrt(sum_x2/n);endm_x;max_x;min_x;rms_x;%按照函数值行参顺序输出结果运行结果

29、:>> m_x,max_x,min_x,rms_x=myfun2(sin(0:6*pi) m_x =max_x =min_x =rms_x = >> m_x,max_x,min_x,rms_x=myfun2(rand(1,200)m_x =max_x = min_x =rms_x =第 3题function v=myvander(x)v1=vander(x);%生成范德蒙矩阵v2=v1'v=flipud(v2);%实现矩阵上下翻转运行结果:>> v=myvander(2:5)v =1 1112 34549162582764125思考题function

30、 x,n=sqrt_a(a)x=;for k=1:100m=x;x=x/2+a/(2*x);if abs(x-m)<=10A(-5) breakendendx;n=k;s=(x-sqrt(a);if s<=10A(-5)disp('正确');elsedisp('错误');end运行结果:>> x,n=sqrt_a(3)正确x =n =5>> x,n=sqrt_a(17)正确x =n =6>> x,n=sqrt_a(113)正确x =8实验5第 1题 .方法1>> t=linspace(0,4*pi,20

31、0);y1=t;y2=sqrt(t);y3=4*pi*exp*t).*sin(t);plot(t,y1,'b',t,y2,'g',t,y3,'r')方法 2>> t=linspace(0,4*pi,200);y1=t;y2=sqrt(t);y3=4*pi*exp*t).*sin(t);t=t,t,t;y=y1,y2,y3;plot(t,y)第 2题>> x=linspace(-6,6,100);y=;for x0=xif x0<=0 y=y,sin(x0);elseif x0<=3y=y,x0;elsey=y,

32、-x0+6;endendplot(x,y);axis(-7 7 -2 4);title(' 分段函数曲线');text(-3*pi/2,1,'y=sin(x)');text(2,2,'y=x');text(4,2,'y=-x+6');第 3题ua=1;ub=cos(-2*pi/3)+sin(-2*pi/3)*i;uc=cos(2*pi/3)+sin(2*pi/3)*i;compass(ua,ub,uc,ua-ub,ub-uc,uc-ua) title( 相量图);第 4题>> z=0:4*pi;x=cos(z);y=s

33、in(z);plot3(x,y,z,'rp');title(' 三维空间曲线');text(0,0,0,'origin');xlabel('X'),ylabel('Y'),zlabel('Z');grid;第 5题( 1)>>x=-3:3;x,y=meshgrid(x);z=-x.A2/10+y.A2/10;mesh(x,y,z);xlabel('X'),ylabel('Y'),zlabel('Z');title(' 立体网状图&#

34、39;);( 2)>>x=-3:3;x,y=meshgrid(x);z=-x.A2/10+y.A2/10;surf(x,y,z);xlabel('X'),ylabel('Y'),zlabel('Z');title(' 立体曲面图');思考题( 3) -2*pi:pi/100:2*pi;y1=2.A*abs(t);y2=2*exp*t);plot(t,y1,'-g');hold on;plot(t,y2,':r');legend(' 曲线 y1',' 曲线 y2&#

35、39;); hold off;grid;实验 6第 1题>> A=3,4,7,2,9,12;x=roots(A)plot(x,'*');grid;x =+ -+第 2题>> A=1,0,0,0,0,-1; x=roots(A) plot(x,'*');grid;x =+-+第 3题%估计零点fplot('xA3+1/x',4);hold on;fplot('2*xA2*sin(x)-5*x*cos(x)',4);hold off;m,n=ginput%建立函数function y=f(x)y=xA3-2*x

36、A2*sin(x)+5*x*cos(x)+1/x;%调用函数>> y1=fzero('fz',y1 =>> y2=fzero('fz',y2 =第 4题%估计零点axis(0,1,0,);fplot('10A(-12)*exp(Ud*10A(-19)/*10人(-23)*300)-1)',0,4);hold on;fplot('(Ud-4)/50',0,4);hold off;m,n=ginput%建立函数function f=myfun(X)Id=X(1);Ud=X(2);f(1)=Id-10A(-12)

37、*exp(Ud*10A(-19)/*10A(-23)*300)-1);f(2)=50*Id-Ud-4;%调用函数>> x=fsolve('myfun',0,optimset('Display','off') x =%验证结果>> K=myfun(x)K =*0第 5题>> p=,;u=10,11,13,14,17,18,22,24,29,34,39;x=polyfit(p,u,3)%得多项式系数t=linspace(0,10,100);y=polyval(x,t);%求多项式得值plot(p,u,'*&

38、#39;,t,y,'r')%画拟和曲线x =实验 7第 1题>>f=50;w=2*pi*f;h=1e-5;tf=40e-3;t=0:h:tf;wt=w*t;y1=4/pi*sin(wt);y2=4/pi*(sin(wt)+1/3*sin(3*wt);y3=4/pi*(sin(wt)+1/3*sin(3*wt)+1/5*sin(5*wt)+1/7*sin(7*wt)+1/9*sin(9*wt);y4=4/pi*(sin(wt)+1/3*sin(3*wt)+1/5*sin(5*wt)+1/7*sin(7*wt)+1/9*sin(9*wt)+1/11*s in(11*wt

39、)+1/13*sin(13*wt)+1/15*sin(15*wt)+1/17*sin(17*wt)+1/19*sin(19*wt);y=square(wt);subplot(2,2,1);plot(wt,y1,wt,y);title('1 次谐波 ');subplot(2,2,2);plot(wt,y2,wt,y);title('1,3 次谐波 ');subplot(2,2,3);plot(wt,y3,wt,y);title('1,3,5,7,9 次谐波 ');subplot(2,2,4);plot(wt,y4,wt,y);title('

40、1,3,5，19 次谐波)；第 2题参数如下:>> R=6e-4;C=17e-4;L=6e-9;模块参数设置:Integrator1 的 Initial condition:15kV在命令窗口为R,L,C 赋值。仿真参数设置如下：Start time:0Stop time:100e-6Solver Type:Variable-stepSolver:ode45Max step size:1e-7Min step size:autoInitial step size:autoRelative tolerance:1e-3Absolute tolerance:1e-6实验六用 matla

41、b 求解常微分方程1 微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为若上式中的系数均与无关，称之为常系数。2 常微分方程的解析解有些微分方程可直接通过积分求解. 例如 , 一解常系数常微分方程可化为, 两边积分可得通解为 . 其中为任意常数. 有些常微分方程可用一些技

42、巧, 如分离变量法, 积分因子法, 常数变异法 , 降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理, 从而他们的求解可归结为求一个特解和相应齐次微分方程的通解. 一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，已给一个阶方程设，可将上式化为一阶方程组反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。3微分方程的数值解法除常系数线性微分方

43、程可用特征根法求解，少数特殊方程可用初等积分法求解外，大 部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题其中所谓数值解法，就是寻求在一系列离散节点上的近似值称为步长，通常取为常量。最简单的数值解法是 Euler法。Euler法的思路极其简单：在节点出用差商近似代替导数这样导出计算公式（称为 Euler格式）他能求解各种形式的微分方程。Euler法也称折线法。Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程（组）初值问题。边值问题

44、采用不同方法，如差分法、有限元法等。数值算法的主要缺点是它缺乏物理理解。4 ,解微分方程的 MATLA瑜令MATLAB3主要用dsolve求符号解析解，ode45,ode23,ode15s求数值解。s=dsolve（'方程 1 ,'方程2' ,'初始条件1','初始条件2'，自 变量）用字符串方程表示，自变量缺省值为t。导数用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形，s为一个符号结构。tout,yout=ode45（ 'yprime ' ,t0,tf,y0）采用变步长四阶Runge-Kutta 法和

45、五阶Runge-Kutta-Felhberg法求数值解，yprime是用以表示 f（t,y） 的M文件名，t0表示自变量的初始值，tf表示自变量的终值，y0表示初始向量值。 输出向量tout表示 节点（t 0,t 1,t n）T,输出矩阵yout表示数值解，每一列对应y的一个分量。若无输出参数，则自动作出图形。ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23与ode45类似，只是精度低一些。ode12s用来求解刚性方程组， 是用格式同ode45。可以用help dsolve, help ode45 查阅有关这些命令的详细信息 例1求下列微分方程的解析解(1)(2)(3)方程(1)求解的MATLAB弋码为:>>clear;>>s