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

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:1x 10故方程 4sinx x 2 = 0 的近似根为：x1=1.8476,x2=0.7623,x3=-2.9117.8.作曲面z = x2 -y2的三维图形.程序: x=-10:0.0