专题2简单的程序设计(数学实验)

《专题2简单的程序设计(数学实验)》由会员分享，可在线阅读，更多相关《专题2简单的程序设计(数学实验)（55页珍藏版）》请在装配图网上搜索。

1、数学建模实验报告专题2 MATLAB 基本使用2018年10月25日目录 3MatLab电子版作业 32.1 练习1 :基本运算 32.2 练习数组及矩阵运算 62.3 程序设计 192.4 练习4 绘图 30-预备实验作业参考解答 36说明：专题 2 作业： Matlab 基本使用：电子版作业（基本运算+数组+程序设计+作图）+数学实验P35-实验内容专题2丁书上命令练习:Diary filename20+30*rand(5) diary off二、MatLab电子版作业2.1 练习1 :基本运算注意：(1),请按要求作题，写出程序并把结果附上，并保存文件1, (5*2+1.3-0.8)*1

2、0A2/25答案：(5*2+1.3-0.8)*10A2/25ans =422, y =sin(10n) e(m3+42) + log423答案：sin(10*pi)*exp(-0.3+4A2)+log(23)/log(4)ans =2.26183, z = 10sin(冗/3)cos(二 /3)答案：z = 10*sin(pi/3)* .cos(pi/3)4.33014, x = sin(223/3), y = xA2 , z = y*10 ;求 x+2y-5z答案：x = sin(223/3); y = xA2;z = y*10 ; x+2*y-5*zans =-37.60255, 4 =2

3、+7i,z2 =2i,z3 =5e2.计算 z = 2-z2 z3答案： z1=2+7*iz1 =2.0000 + 7.0000i z2=2*i;z3=5*exp(2*pi*i)z3 =5.0000 - 0.0000i z=z1*z2/(z2+z3)z =-2.1379 + 1.6552i6,建立起始值=3,增量值=5.5,终止值=44的一维数组x(填写程序语句)答案：3: 5.5: 447,建立等差一维数组x:首项为0,末项为冗，项数为15 (填写程序语句)答案:linspace(0,pi,15)3 2 38,计算行列式的值A = 4 2 67 8 1答案： A=3 2 3;4 2 6;7

4、8 13 234 267 81 det(A)ans =-83231119,矩阵A=426，矩阵B=222 ;分别求出A x B及A与B中对应元素781 -333-之间的乘积的程序语句。答案：A=3 2 3;4 2 6;7 8 1;B=1 1 1;2 2 2;3 3 3; A*B,A*B ans =16 16 1626 26 2626 26 26 ans =32384 1221 24310,写出计算Sin(36 o)的程序语句.答案：sin(36*pi/180)ans =0.58782.2练习数组及矩阵运算34-1191065074-161-47168A =2-45-6123- 36-78111

5、.输入矩阵4)a4 =7-612-16-815. 改变 A 的元素 A( 4， 5)为0 值。答案： A(4,5)=0A =34 -11 -9 1065074 -161 -47 -16 -82 -45 -60 -8-36-78-118-4913016.改变A的第6列为b答案：A(:,6)=bA =34-11-916507431 -47-1652 -45-607-36-78-198-49131117.建立三角函数表阵S。第一列是自变量X的取值，从0-2冗共取20个值。第 二列是对应Y=cos X的值。答案：X=linspace(0,2*pi,20); Y=cos(X); S=X;Y0.3307

6、0.94580.6614 0.78910.9921 0.54691.3228 0.24551.6535 -0.08261.9842 -0.40172.3149 -0.67732.6456 -0.87952.9762 -0.98643.3069 -0.98643.6376 -0.87953.9683 -0.67734.2990 -0.40174.6297 -0.08264.9604 0.24555.2911 0.54695.6218 0.78915.9525 0.94586.2832 1.000018. 以红色的+号 画 y 的图形答案：plot(X,Y,r+)19. 分别取 a2 的 1、 2

7、 行作成两个向量A1、 A2答案：A1=a2(1,:),A2=a2(2,:)A1 =74 -16A2 =-6 12 -820 .作向量的点积C1=A1- A2答案： C1=dot(A1,A2)C1 =134或者： C1=A1*A2C1 =13421 .作向量的叉积C2=A A2C2=cross(A1,A2)答案：C2 =160 152 10822 3程序设计1,编辑函数文件hanshu.m,该函数是f (x) = x3 -3x2 -x + 3 ,并用它来计算f(0.34)、f(3)、f (6.87) f (3) - f (2) f2 (-4)答案：函数文件hanshu.m如下：function

8、 y=hanshu(x)y=xA3-3*xA2-1*x+3;计算程序及结果如下：x=0.34 3 6.87x =0.3400 3.0000 6.8700 hanshu(x) ans =2.35250 178.7820 hanshu(3)-hanshu(2)*(hanshu(-4)A2ans =330752,编程求一个四位数abca,使abca = (5c + 1)2。答案：(填写程序语句和结果)3,利用二分法编程求方程x3-3x2-x + 3 = 0在0,3内的根.答案：(填写程序语句和结果)先编辑函数文件hanshu.m如下：function y=hanshu(x)y=x.A3-3*x.A2

9、-x+3 ;然后编写m文件如下：cleara=0; b=3;while abs(a-b)10A(-6)if hanshu2(a)*hanshu2(b)0 b=(a+b)/2;else a=(a+b)/2;end end gen123=a 结果：gen123 =3.0000*4,利用二分法编程求方程tan(x)-x-0.5=0在0,1.5内的根.答案：(填写程序语句和结果)先定义函数：fin.mfunction y=fin(x) y=tan(x)-x-0.5;然后编辑m文件qiugen.m:fplot(fin,0 1.5)a=0;b=1.5;while abs(a-b)10A(-6)if fin

10、(a)*fin(a+b)/2)0 a=(a+b)/2;else b=(a+b)/2;endenda结果：0.97505,列出所有的水仙花数，水仙花数是一个三位数，其各位数字立方和等于该数本身. 例如：153 =13 +53 +33答案：(填写程序语句和结果)程序1:for a=1:9for b=0:9for c=0:9abc=100*a+10*b+c;if abc=a.A3+b.A3+c.A3 disp(abc)endendendend程序2:k=1;for a=1:9for b=0:9for c=0:9abc=100*a+10*b+c;if abc=a.A3+b.A3+c.A3shuixia

11、nhua(k)=abc;k=k+1;endendendendshuixianhua结果：shuixianhua =153 370 371 407程序3:for m=100:999m1=fix(m/100);% 求m的百位数字 m2=rem(fix(m/10),10);% 求 m的 10 位数字 m3=rem(m,10);if m=m1*m1*m1+m2*m2*m2+m3*m3*m3 disp(m)endend6.利用-&1-1 +1+ 公式求n的近似值,直到误差小于10”。43 5 7k=1;s=0;while 4/(2*k-1)10A(-6)s=s+(1)A(k+1)/(2*k1)；k=k+

12、1;endpai=4*s结果：pai =3.14159215358972冗 tan = 1【 中学时我们就已经知道： 4，从而n=4arctan1 ,如果我们应用泰勒3572n 1公式将 展开，就可以得到 arctan x = x - -+ -+(-1)na- 3572n-1把 x =1 代入，就得到一=1 + +*( 1)2 + 43 5 72n -1这个简单的公式就可以用循环实现来实现圆周率的近似计算(虽然它并不真正适用于计算圆周率)】7,写出小于5000的、立方的末四位是8888的所有自然数的程序答案：(填写程序语句和结果)clear;clcn=5000;ii=1;for kj=1:50

13、00if rem(kjA3,10000)=8888shu(ii)=kj;ii=ii+1;endendshu结果：shu =19424442*8,先思考如何用推理的方法(1)求Fibonacc数列中第2007个数除以6的余 数(2)若将Fibonacc数列按如下方式分组：(1)、(1,2)、(3,5,8)、(13,21, 34, 55)、问2007组数的和除以6的余数是几？然后再用 matlab编程 验证你的结果。答案：(填写你的思考过程和你的结果以及 matlab程序语句)提示：(1)可以发现Fibonacc数列中的数除以6的余数呈周期性的变化(2)第2007组数开始是几？结尾是几？考虑一个周

14、期内的和除以6的余数。9.应用e=1+1十1+.+工近似求e之值，要求误差小于10： 2! n!程序：k=1;s=1;while 1/prod(1:k)10A(-6)s=s+1/prod(1:k);k=k+1;ende=s结果：e =2.7182815255731910,定义一个函数shjiech计算双阶乘，要求对输入的一个整数N,该函数能立即计算出其双阶乘。并用它来计算 200! !答案：(填写程序语句和结果)function jiech=shjiech(n)if rem(n,2)=0as=2:2:n;else as=1:2:n;endjiech=prod(as);11,定义多元函数 f(x

15、,y) = (x2+y2)eYx2*2),并用它来计算 f(-2,3)、f(3,4)、f(0,0),并画出它的图形。答案：(填写程序语句和结果)函数文件hanshu4.mfunction fxy=hanshu4(a,b)fxy=(aA2+bA2)*exp(-(aA2+bA2);计算：hanshu4(-2,3)ans =2.938428229075371e-005 hanshu4(3,4),hanshu4(0,0)ans =3.471985966241005e-010ans =0绘图程序：x=-2:0.01:2;y=x;xx,yy=meshgrid(x,y);zz=(xx.A2+yy.A2).*

16、exp(-(xx.A2+yy.A2); mesh(zz)结果：*12,编写函数文件SQRT1.M迭代法求 石之值，迭代公式是Xn书(Xn+2),2 Xn误差小于10-6 0答案：(填写程序语句和结果)function genhaoa=SQRT1(a)an=1;while 1/2*abs(a/an-an)10A(-6) an=(an+a/an)/2;endgenhaoa=an;比如：SQRT1(2)ans =1.41421356237469SQRT1(2456)ans =49.55804687847163补充1.求100,200 之间的素数I.for循环程序：for m=1012200k=fix

17、(sqrt(m);for i=2:k+1if rem(m,i)=0break;endif i=k+1% 除遍了所有（ 2到 K+1 ）disp(int2str(m)endend或者： (以数组的形式显示)u=0;for m=101:2:200k=fix(sqrt(m);for i=2:k+1if rem(m,i)=0 break;endendif i=k+1u=u+1;a(u)=m;endendb=aII.while 循环程序：m=101;while m200k=fix(sqrt(m);i=2;while i=k+1disp(int2str(m)endm=m+2;end或：i=1;x=;for

18、 j=100:200if isprime(j)=1%此命令可直接判断j 是否为素数x(i)=j;endi=i+1;end或：N=input( 请输入数值： m=)x=100:N;for j=2:sqrt(N);i=find(rem(x,j)=0);x(i)=;endx2. 若一个数等于它各个真因子之和，则称该数为完数，如 6=1+2+3 ，所以 6 就是完数，求1,500 之间的全部完数。程序： for m=1:500s=0;for k=1:m/2if rem(m,k)=0s=s+k;endendif m=sdisp(m)endend3. 求鸡兔同笼问题，头共36，脚共100.求鸡、兔各多少？

19、程序：i=1;while 1if rem(100-i*2,4)=0&(i+(100-i*2)/4)=36 break;endi=i+1;enda1=ia2=(100-2*i)/4a1 =22a2 =144 .用起泡法对10 个数由小到大排序，即相邻两个比较，将小的调到前头。程序 1 ：x=input(x=)n=length(x);for i=1:n% 比较 n 轮for j=1:n-i % 两个比较第一轮比较n-1 次， 。 。if x(j)x(j+1)a=x(j);x(j)=x(j+1);x(j+1)=a;endendendy=x程序 2.function y=qp(x)m,n=size(x

20、);for i=1:n% 比较 n 轮for j=1:n-i % 两个比较第一轮比较n-1 次， 。 。if x(j)x(j+1)a=x(j);x(j)=x(j+1);x(j+1)=a;endendendy=x5 .有一个4X5矩阵，编程求出最大值及其所处的位置程序 1 ：x=input(x=)c,t=max(x);y,i=max(c);t(i);disp(最大值为：,num2str(y)disp(位置为：,num2str(t(i),行,num2str(i),列)程序 2 ：function max-1(x)c,t=max(x);y,i=max(c);t(i);disp(最大值为：,num2s

21、tr(y)disp(位置为：,num2str(t(i),行,num2str(i),列)2.4练习4绘图1,在0,4pi画sin(x),cos(x)(在同一个图象中)； 其中cos(x)图象用红色小圆圈画.并在函数图上标注“y=sin(x) ，“y=cos(x) ” ,x轴,y轴，标题为“正弦余弦函数图象”.答案：函数文件如下：t=0:0.02:4*pi;plot(t,sin(t),m*,t,cos(t),rO)text(pi/3,sin(pi/3),y=sin(x)text(pi/3,cos(pi/3),y=cos(x)title(正线函数和余弦函数)xlabel(x 轴)ylabel( y

22、轴)结果显示：正线函数和余弦函数102460101214X轴060.40.2*- 0.2- 0 4- 06- 0.3222,画出椭圆人+工=1。169答案：函数文件如下： 刎椭圆 tt=0:0.02:2*pi;xx=4*cos(tt);yy=3*sin(tt);axis equal plot(xx,yy,m*)title( 椭圆 xA2/16+yA2/9=1)答案2:函数文件：xx=-4:0.02:4;yy=3*sqrt(1-xx.A2/16);plot(xx,yy,m*)axis equal hold on plot(xx,-yy,r+)title( 椭圆 xA2/16+yA2/9=1)椭圆

23、/16+尸尼=1-4-3-2-101234x = (1 2cos u)3,绘制螺旋线.BO.60、20.210四叶玫瑰线9015,用subplot分别在不同的坐标系下作出下列图形，并为每幅图形加上标题。第一副：标题：”概率曲线2-xy =e第二副：标题“双扭线”=4cos 2t3t第三副：标题“叶形线”x 二1t33t2y二门第四副：标题：“曳物线”xdn1-；1-y2+ 可y答案：x=-5:0.1:5;y=exp(-x.A2);subplot(2,2,1)plot(x,y)title(概率曲线)tt=-pi/4:0.03:pi/4;pp=2*sqrt(cos(tt*2);subplot(2,

24、2,2)polar(tt,pp)hold onpolar(tt,-pp)title( 双扭线)subplot(2,2,3)t=0:0.01:20;x=3*t./(1+t.A3);y=3*t.A2./(1+t.A3);plot(x,y)title( 叶形线)subplot(2,2,4)yyy=0:0.01:1;xxx=log(1+sqrt(1-yyy.A2)./(yyy+eps)-sqrt(1-yyy.A2); plot(xxx,yyy)hold onxxx=log(1-sqrt(1-yyy.A2)./(yyy+eps)+sqrt(1-yyy.A2); plot(xxx,yyy)title( 曳

25、物线)概率曲线三、31.设有分块矩阵A =10 2用和对角矩阵，试通过数值计算验证程序：E=eye(3);R=rand(3,2);O=zeros(2,3);S=diag(1,2);A=E,R;O,S;B=A*A,C=E R+R*S;O S*S,a=B=CR32。，其中E, R, 0, S分别为单位矩阵、随机矩阵、零矩阵S22 JE R+ RSS2结果：B =1.0000000.82051.058601.000001.78732.4395001.00000.11580.02960001.0000000004.0000C =1.0000000.82051.058601.000001.78732.4

26、395001.00000.11580.02960001.0000000004.0000a =11111111111111111111111112.用命令magic(n)生成幻方矩阵，通过计算研究它的性质，如行和、歹U和、两条对角 线和等(可以利用命令 diag, sum, fliplr, flipud ,其用法可以查阅 MATLAB帮助系统). 程序：A=magic(4),s1=sum(A),s2=sum(A), s3=sum(diag(A),s4=sum(diag(fliplr(A)结果：A =16231351110897612414151s1 =34343434s2 =34343434s3

27、 =34s4 =343 .设y1二-1/(1+x2),y2 = exp(-x2/2) ,y3 =sin(2x) ,y4=Y4-xx 在22 2内适当离故化,计算 y1 +y2, y1y2, y3/y2, (5y4 -y1)/y2.程序：x=-2:0.5:2;y1=-1./(1+x.A2);y2=exp(-x.A2 /2);y3=sin(2*x);y4=sqrt(4-x.A2); a=y1+y2,b=y1.*y2,c=y3./y1,d=(5*y4-y1)./y2.A2结果：-0.06470.0170-0.06470.10650.082500.08250.10650.0170b =-0.0271-

28、0.0999-0.0271c =-3.78400.45863.7840-0.3033-0.7060-1.0000-0.7060-0.3033-0.09991.81861.05180-1.0518-1.8186-0.4586d =10.919665.674810.919624.900213.459711.000013.459724.900265.67484 .自己选择一非负单调递减序列an a2,，an定0, a1远大于an,用从1到n和从n至U 1两种顺序计算至3k,观察哪个更准确些，分析原因。i=1程序与结果： syms x symsum(1/x,1,100) ans =1446663627

29、9520351160221518043104131447711 /2788815009188499086581352357412492142 272 symsum(1/x,100,1) ans =-291248328005999177069358804052937859600407 /69720375229712477164533808935312303556800另外：s=0; for k=1:100000s=s+1/k;ends(1,2,3,，100000)=,ss1=0;for k=100000:-1:1s1=s1+1/k;ends1(100000,99999,99998,1)=,s1a

30、ns =s(1,2,3,100000)=s =12.09014612986334 ans =s1(100000,99999,99998,，1)=s1 =12.09014612986341for k=1:100000s=s-1/k;s1=s1-1/k;ends,s1 s =-7.305902946150560e-014s1 =-2.283990460953331e-0165.对 In = f0xnex 1dx ( n = 0, 1,2, )证明如下递推公式：1Io = 1 e ) In = 1 nIn1, n = 1,2,用递推公式计算Ii, I2,，In,观察n多大时结果就不对了(考虑一个简单

31、的判断结果错误的标准)，为什么会出现这种情况.如果将递推公式反过来，即In 1 = ( 1 -In )/n .从In倒过来计算In-1, , I1, Io,而In由下式估计x 1.1,(1),(2)/ -x 1,1(0mxn_1e ) n+1= In In =exp(-1)/(n+1).& (in0In(n+1)=1-n*In(n);i=i+1;endIn,i反向推导：IIn(19)=(exp(-1)+1)/2;for n=18:-1:1IIn(n)=(1-IIn(n+1)/nendIInclc,clearI0=1-exp(-1);In(1)=I0;i=1;while In0In(i+1)=1

32、-i*In(i);i=i+1;endIn,iIIn(19)=(exp(-1)+1)/2;for n=18:-1:1IIn(n)=(1-IIn(n+1)/n endIInduibizhen=In,IIn 结果：duibizhen =0.632120558828560.367879441171440.264241117657120.207276647028650.170893411885380.145532940573080.126802356561520.112383504069360.100931967445090.091612292994170.083877070058290.0773522

33、29358780.071773247694640.066947779969720.062731080423870.059033793641900.055459301729570.05719187059731-0.029453670751540.632120558828560.367879441171440.264241117657120.207276647028650.170893411885390.145532940573070.126802356561600.112383504068800.100931967449580.091612292953740.083877070462640.07

34、7352224911000.071773301068060.066947086115270.062740794386170.058888084207400.057790652681650.017558904411900.683939720585726 用几种方法作x2， x3， x4， x5 的图形，如一个图上画几条曲线， 用 subplot 作多幅图形等，考虑如何画上x 轴， y 轴，并在图上加各种标注程序 1： x=1:0.01:3;y1=x.A2;y2=x.A3;y3=x.A4;y4=x.A5;plot(x,y1,x,y2,x,y3,x,y4) gtext(xA2),gtext(xA3)

35、,gtext(xA4),gtext(xA5)结果 1：250程序2:x=1:0.01:3;y1=x.A2;y2=x.A3;y3=x.A4;y4=x.A5;subplot(2,2,1),plot(x,y1),title(xA2)subplot(2,2,2),plot(x,y2),title(xA3)subplot(2,2,3),plot(x,y3),title(xA4)subplot(2,2,4),plot(x,y4),title(xA5)gtext(xA2),gtext(xA3),gtext(xA4),gtext(xA5)结果：,2,3x=-2:0.05:2;y=x.A2;x.A3;x.A4;

36、x.A5;z=0*x;plot(x,y,x,z);grid on;legend(xA2,xA3,xA4,xA5)7.用作图法求x2 = 81nx和4sinx -x -2 = 0的根的近似值.程序1:x=1:0.01:3;y1=x.A2;y2=8*log(x);plot(x,y1,x,y2),gtext(xA2),gtext(8Inx)结果1：9故方程 x2 = 81nx 的近似根为 x1=1.1957,x2=2.9348程序2:x=-4:0.01:4;y=4*sin(x)-x-2;z=0*x;p1ot(x,y,x,z),axis(-4 4 0 1),gtext(4sinx i -2)结果2:1

37、x 10故方程 4sinx x 2 = 0 的近似根为：x1=1.8476,x2=0.7623,x3=-2.9117.8.作曲面z = x2 -y2的三维图形.程序: x=-10:0.01:10;y=x;X,Y=meshgrid(x,y); Z=X.A2-Y.A2; mesh(X,YZ)结果：9.图20中电阻R = (R1+R2)为分压器，Rl为负载电阻，试将分压比y = Ul/U表示为x = R2/R 和2 = Rl/R的函数，并以a为参数(a =10, 1,0.1),作函数y(x, n)的图形，对结果作出解释.Rp1-1u Uk 41分析：U/Ul=_Rl - R2/(Rl+R2)RL(R

38、-R2)+R l - R2/(R l+R2)丫图20 R2=Rx,RL=R.化 简得：U/Ul=R/R 2+R/R L-R2/R L=(a+x-x2)/xa由此得：y=xa/(a+x-x 2)程序： x=0:0.01:1;a=10;y1=a*x./(a+x-x.A2); a=1;y2=a*x./(a+x-x.A2); a=0.1;y3=a*x./(a+x-x.A2); plot(x,y1,x,y2,x,y3),gtext(a=10),gtext(a=1),gtext(a=0.1)结果：化简彳导:y=x-x2(1-x)/x(1-x)+a,故当 a越大日X( x2(1-x)/x(1-x)+a越趋近

39、于 0,此时 y=x。10.建立M文件作以下计算：(1)自然数n的阶乘；(2) n中取m的组合；(3)已知任意两个多项式(不一定同阶)的系数，求它们的和.程序1:建立函数 M文件如下：function m=jiecheng(n)p=1；for i=1:np=p*i；endm=p;结果1：m=jiecheng(5)120程序2:建立函数 M文件如下:function q=zuhe(n,m) p=1;if (n=m & n=0 & m=0) m=min(n-m,m);for i=1:mp=p*(n-i+1)/i;endendq=p;结果2：q=zuhe(5,3)q =10程序3：建立函数M 文件如

40、下：function p=qiuhe(p1,p2) a=length(p1);b=length(p2);if abp2=zeros(1,a-b),p2;endp=p1+p2;结果3：p1=1 0 9 7 6; p2=2 4; p=qiuhe(p1,p2)p =109910P36,第 10 题：10(2)function y=zhuhe(n,m) y1=1;y2=1;y3=1;for i=1:ny1=y1*i;endfor j=1:my2=y2*j;endfor k=1:(n-m)y3=y3*k;endy=y1/y2/y3;10(3):function f=sumpoly(p1,p2)n=abs(length(p1)-length(p2);if length(p1)length(p2)p22=p2,zeros(1,n); p12=p1;elsep12=p1;p22=p2;endf=p12+p22;