数学实验作业汇总

1、（1）产生一个5阶魔方矩阵M：M=magic(5)（2）将矩阵M的第3行4列元素赋值给变量t：t=M(3,4)（3）将由矩阵M第2，3，4行第2,5列构成的子矩阵赋给变N：N=M(2:4,2:3:5)（4）将由矩阵M的前3行赋给变量N： N=M(1:3,:)（5）将由矩阵M的后3列赋给变量N： N=M(:,end:-1:end-2)（6）提取M的主对角线元素，并以这些对角线元素构成对角矩阵N： N=diag(diag(M)或N=tril(triu(M)(7)随机产生1000个100以内的整数赋值给变量t： t=round(rand(1,1000)*100

2、)（8）随机产生100*5个100以内的实数赋值给变量M：M=rand(100,5)*100（1）删除矩阵M的第7个元素  M(7)=（2）将含有12个元素的向量t转换成3*4的矩阵：reshape(t,3,4)（3）产生和M同样大小的单位矩阵： eye(size(M)（4）寻找向量t中非零元素的下标：find(t)（5）逆序显示向量t中的元素：t(end:-1:1)（6）显示向量t偶数位置上的元素： t(2:2:end)（7）利用find函数，将向量t中小于10的整数置为0：t(find(t<10&rem(t,1)=0)=0（8）不用fi

3、nd函数，将向量t中小于10的整数置为0： t(t<10&rem(t,1)=0)=0（9）将向量t中的0元素用机器0（realmin）来代替： t(find(t=0)=realmin（10）将矩阵M中小于10的整数置为0： M(find(M<10)&rem(M,1)=0)=02、写出完成下列操作的命令及结果。（1）将150这50个整数按行优先存放到5*10的矩阵中，求该矩阵四周元素的和；>> t=1:10; >>  M=t;t+10;t+20;t+30;t+40M =1 &

4、#160;   2    3    4    5    6    7     8    9  10   11    12    13    14    15   

5、; 16    17    18    19    20   21    22    23    24    25    26    27    28    29    30 &#

6、160; 31    32    33    34    35    36    37    38    39    40   41    42    43    44    45 

7、;   46    47    48    49    50>>  N=M(2:4,2:9)N =12    13    14    15    16    17    18    19  22&#

8、160;   23    24    25    26    27    28    29  32    33    34    35    36    37    38  

9、60; 39  >> sum(sum(M)-sum(sum(n)ans =   663 2）n取100、1000、10000，求序列1、1/2、1/31/n的和。>> n=100;>> t=1:n;>> format rat>> M=t.-1;>> S=sum(M)S =2630/507>> n=1000;>> t=1:n;>> format rat>> M=t.-1;>> S=sum(M)S =1804

10、/241>> n=10000;>> t=1:n;>> format rat>> M=t.-1;>> S=sum(M)S =1106/1131. 在同一坐标系下绘制y1=sin(t),y2=sin(2t),y3=sin(3t),其中y1的数据点用星号，线形为黑色虚线，y2的数据点用方块，线形为红色实线，y3的数据点用小圆圈，线形为蓝色点线。（要求采用一次绘出和逐次填加两种方式完成绘图）>> t=linspace(0,2*pi,100);>> y1=sin(t); >> y2=sin

11、(2*t); >> y3=sin(3*t); >> plot(t,y1,*k:,t,y2,sr-,t,y3,ob-.)>> t=linspace(0,2*pi,100);>> y1=sin(t);>>  plot(t,y1,*k:)>> hold on>> y2=sin(2*t);>> plot(t,y2,sr-)>> hold on>> y3=sin(3*t);>> plot(t,y3,ob-.)>> hold off2

12、. 分别用plot和fplot函数绘制y=sin（1/x）的曲线，分析两曲线的差别>> x=linspace(0,1/(2*pi),100); >>  y=sin(x.-1);  >>  plot(x,y,*-) >> fplot(sin(x.-1),0,1/(2*pi),o-)两曲线的差别:plot曲线在确定自变量x的取值间隔时采用平均间隔，图像不是十分准确；fplot曲线自动取值，在函数值变化平稳时，它的数值点会自动相对稀疏一点，在函数值变化剧烈处，所取点会自动密集一点，所以曲线更加光滑准确。

13、6.已知曲面方程f（x,y）= ,x -1.5，1.5，y -2.5，2.5，用建立子窗口的方法在同一图形窗口绘制出三维线图，网线图，曲面图。>>  x=-1.5*pi:pi/50:1.5*pi;>> y=-2.5*pi:pi/50:2.5*pi;>> X,Y=meshgrid(x,y);>> Z=sin(sqrt(X.2+Y.2)./sqrt(1+X.2+Y.2);>> subplot(1,3,1);plot3(X,Y,Z);>> subplot(1,3,2);mesh(X,Y,Z);>> subp

14、lot(1,3,3);surf(X,Y,Z);8.将peaks函数生成的最高峰削去，并用色图矩阵“cool”修饰。>> x,y,z=peaks(30);>> x1=x(1,:);y1=y(:,1);>> i=find(y1>1&y1<3);>> j=find(x1>-1&x1<1);>> z(i,j)=NaN*z(i,j);>> surf(x,y,z)>> colormap(cool)3. 定义一个函数，函数的自变量为整数n，函数的功能是：随机产生n个三位整数，将其中小于

15、平均值的数用0代替。function mean,x=ff (n)  x=floor (100+899*rand (1,n);  m=length (x);  mean=sum (x)/m;  x (x<mean)=0;4. 编写函数，用来求下列函数的和，并给出n分别为100,1000,10000时，下列各式的值。 function y=s(n)y=1;for i=1:1:n x=4*i2/(4*i2-1); y=y*x;enddisp(y)s(100)=1.5669s(1000)=1.5704s(10000)=1.57085. 通过命令文件实现：随

16、机产生20个数，输出其中的最大数和最小数。通过函数文件实现：随机产生n个数，输出其中的最大最小数。命令文件>> t=rand(1,20);>> disp('max=');disp(max(t)max= 0.7942>> disp('min=');disp(min(t)min=0.0503函数文件function f3(n)t=rand(1,n);disp('max=');disp(max(t);disp('min=');disp(min(t);end3.求下列函数的一阶和二阶导数 >>

17、; syms x                                     >> diff(2/tan(x)+cos(x)/3,x,1) ans =- sin(x)/3 - (2*(tan(x)

18、2 + 1)/tan(x)2>> syms x                                     diff(2/tan(x)+cos(x)/3,x,2)4.求积分 >> syms x&

19、#160;int(sqrt(exp(x)+1),x) ans =2*(exp(x) + 1)(1/2) + 2*atan(exp(x) + 1)(1/2)*i)*i5.求下列级数的和>> syms n>> s=symsum(-1)(n+1)*1/n,1,inf)s =log(2)6.求函数在x=0处的泰勒展开式>> syms x>> taylor(exp(x)+exp(-x)/2,x,5,0)ans = x4/24 + x2/2 + 11. 利用randn函数声称符合正态分布的10*5随机矩阵A，进行以下操作：(1).A的各列

20、元素的均值和标准方差(2).A的最大元素及其所在位置(3).A的每行元素的和以及全部元素之和(4).分别对A的每行元素按升序排序(5).将A中的每行元素的总和按从大到小的顺序存入line_sum中，相应的行号存入line_num中>> A=randn(10,5);>> a1=mean(A)>> a2=std(A)>> AA=max(max(A)>> i j=find(A=AA)>> a3=sum(A,2)>> a4=sum(sum(A)>> a5=sort(A,2)>> line_sum

21、,line_num=sort(sum(A,2),'descend')2、补充题：利用导入向导（或借助函数imread）导入一幅单色图片存入变量ima_data中，然后依次完成下列操作：（1）用imshow函数显示图片；（2）删除图片前若干行（例如前100行）再次显示该图片。（3）将图片上、下翻转再次显示图片。先找到一个.bmp的文件，把它放入工作目录下，并修改名称为1.bmp,执行下列操作。ima_data=imread(1.bmp);(1)imshow(ima_data);(2)a=ima_data(101:end,:);imshow(a);(3)imshow(flipud(

22、ima_data);3.下表所示是090度内某些数的正弦近似值x度0153045607590Sinx00.25880.50.70710.8660.96591利用线性、样条差值求x=20、40、80度时正弦值，这两种方法哪个好？为什么实验步骤：利用inerp1函数先分别求出线性插值和三次样条插值所得到的y11和y12，再利用sin（x）函数得到准确的y1，比较y11和y1，y12和y12，不难得出结论。所用语句 clear;clc;x=0 15 30 45 60 75 90./180.*pi;y=sin(x);x1=20 40 80./180.*pi;y11=interp1(x,y,x1,lin

23、ear);y12=interp1(x,y,x1,spline);y1=sin(x1);主要结果 y11= 0.3392    0.6381    0.9773;y12=0.3420    0.6428    0.9849;y1=0.3420    0.6428    0.9848;4.已知某次实验测得数据如下：x11.41.82.22.633.43.84.24.65y0.870.525.213.5114.2

24、919.4314.1341.5313.9158.5614.99x5.45.86.26.677.47.88.28.699.4y130.4744.8221.2543.15281.25200.09177.93344.53509.84531.07260.49（1）请用3次多项式进行拟合，并给出拟合函数在0、0.5、1、1.59、9.5处的值（2）估计用几阶多项式拟合的效果较好，并说明理由。4.(1)clear;clc;x=1:0.4:9.4;y=0.87 0.52 5.21 3.51 14.29 19.43 14.13 41.53 13.91 58.56 14.99 130.47 44.82 21.2

25、5 43.15 281.25 200.09 177.93 344.53 509.84 531.07 260.49;x1=0:0.5:9.5;p=polyfit(x,y,3);y1=polyval(p,x1);主要结果:y1=50.55 33.03 18.91 8.38 1.61 -1.23 0.05 5.62 15.65 30.32 49.80 74.28 103.92 138.91 179.41 225.61 277.67 335.79 400.12 470.85(2) 19阶拟合效果最好。理由通过编写差方和函数(基于最小二乘原理）

26、f(n)f(n)函数如下：function tz=f(n)t=;x=1:0.4:9.4;y=0.87 0.52 5.21 3.51 14.29 19.43 14.13 41.53 13.91 58.56 14.99 130.47 44.82 21.25 43.15 281.25 200.09 177.93 344.53 509.84 531.07 260.49;for i=1:n    p=polyfit(x,y,i);    y1=polyval(p,x);    c=sum(y-y1).2,2);&#

27、160;   t=t c;endtz=find(t=min(t);令n=22（一共22组数据）f函数值最小时是19阶时所以得出结论19阶多项式拟合效果最好。再用拟合图像（p=polyfit(x,y,19)，plot(x,y,:o,x,polyval(p,x),-*)）也可以看出19阶多项式拟合效果最好。2、自行练习题。下列填空题是期中考试出错比较多的题目，请认真考虑并上机调试。（6）逆序显示向量t中的元素：（7）显示向量t偶数位置上的元素 ：（9）删除向量t中最小的5个数：（17）将150按列优先存放到5*10的矩阵M中：（18）求矩阵M最大值所在的位置：（19）

28、统计字符串S中小写字母的个数：（20）设A是n阶0、1方阵，A边界上1的个数：(6).t(end:-1:1)(7).t(2:2:end)(9).M=sort(t)     a=find(t<M(6)     t(a)=(17).t=1:5:46       M=t;t+1;t+2;t+3;t+4(18).i,j=find(M=max(max(M)(19).a=find(s>=a&s<=z)   num

29、=length(a)(20).B=A(2:end-1,2:end-1)   num=sum(sum(A)-sum(sum(B)1.分别用矩阵求逆、矩阵除法以及矩阵分解求线性方程组的解 矩阵求逆>> A=2,2,-1,1;4,3,-1,2;8,3,-3,4;3,3,-2,-2;>> b=4,6,12,6;>> inv(A)*b运用左除运算符>> A=2,2,-1,1;4,3,-1,2;8,3,-3,4;3,3,-2,-2;>> b=4,6,12,6;>> x=Ab运用矩阵分解>> A=2,2,

30、-1,1;4,3,-1,2;8,3,-3,4;3,3,-2,-2;>> b=4,6,12,6;>> Q,R=qr(A);>> x=R(Qb)4.在区间30,50内，求 的零点。>> f=5*sin(x)-2*(log(x)/log(3)+1.8;>> ezplot(f,30,50)>> fzero(f,33)ans =   32.5547>> fzero(f,34)ans =   33.3960>> fzero(f,38)ans =   3

31、9.0426>> fzero(f,39.4，39.5)ans =   39.4785则方程有四个零点6. 给出实验数据如下：x2345678910111213141516y6.248.209.589.609.6010.029.939.9910.4710.5910.6010.8010.6010.9010.75试分别用 做拟合形式，求出a和b及拟合曲线，并画图进行比较。>> x=2:16;>> y=6.24,8.20,9.58,9.60,9.60,10.02,9.93,9.99,10.47,10.59,10.60,10.80,10.60,10

32、.90,10.75;>> X=1./x;>> Y=log(y);>> P=polyfit(X,Y,1)P =   -1.1552    2.4629>> exp(2.4629)ans =   11.7388则a=11.7388     b=-1.1552作图：>> Y1=polyval(P,X)>> y1=exp(Y1);>> plot(x,y,:o,x,y1,-*)>> x=2:16;&g

33、t;> y=6.24,8.20,9.58,9.60,9.60,10.02,9.93,9.99,10.47,10.59,10.60,10.80,10.60,10.90,10.75;>> Y=1./y;>> X=1./x;>> P=polyfit(X,Y,1)P =    0.1384    0.0815则a=0.0815   b=0.1384 作图：>> Y1=polyval(P,X);>> y1=1./Y1;>> plot(x,

34、y,:o,x,y1,-*)3.求下列方程或方程的根在指定点的近似根 ，初值 function f=myFun(x)f(1)=sin(x(1)+x(2)2+log(x(3)-7;f(2)=3*x(1)+2x(2)-x(3)3+1;f(3)=x(1)+x(2)+x(3)-5;>> X=1,1,1;>> op=optimset(display,off);>> x=fsolve(myfun,X,op)x =    0.5991    2.3959    2.00502. 已知

35、，求y的单调增区间和y的极值>> fplot(2*sin(x)+cos(2*x),0,pi/2)>> syms x>> f=2*sin(x)+cos(2*x);>> s=diff(f)s = 2*cos(x) - 2*sin(2*x)>> fzero(2*cos(x) - 2*sin(2*x),0.5)ans =    0.5236由图知单调递增区间为0,0.5236；将ans的值代入原式中，得y的极值为1.5。3. 求解线性约束最优化问题 function f=fop(x)f=0.5*x(1)

36、2+x(2)2-x(1)*x(2)-2*x(1)-6*x(2);>> x0=0.5;0.5;>> A=1,1;-1,2;2,1;>> b=2;2;3;>> lb=0;0;>> options=optimset(display,off);>> x,f=fmincon(fop,x0,A,b,lb,options)x =    0.6667    1.3333f =   -8.22221、 请你构造一个生成素数的公式,并将你的工作与Euler的工作比

37、较。 采用素数生成公式p=n2-79*n+1601(1)编写函数f(x),用来计算素数多项式生成公式，在100以内和1000以内，产生素数的百分比,程序如下：function tz=f(x)n=0:x(1,3);t=n.2+x(1,1)*n+x(1,2);t1=find(isprime(t);tz=length(t1)/length(n);end(2)代入Euler公式系数x1=1 41 100,x2=1 41 1000与p=n2-79*n+1601系数y1=-79 1601 100,y2=-79 1601 1000比较得到结果f(x1)=0.8614;f(x2)=0.5814;f(y1)=0

38、.9505;f(y2)=0.6014;所以可得结论该公式比Eluer的公式生成素数的概率要高；2、 研究百万以内素数的间隔规律。a=primes(1000000);b=a;b(1)=;a(length(a)=; t=b-a; plot(a,t,.); t1=unique(t)           %求相邻素数间的间隔值t1 =  Columns 1 through 14     1   

39、  2     4     6     8    10    12    14    16    18    20    22    24    26  Columns 15 through

40、 28    28    30    32    34    36    38    40    42    44    46    48    50    52    54

41、60; Columns 29 through 42    56    58    60    62    64    66    68    70    72    74    76    78    80&#

42、160;   82  Columns 43 through 52    84    86    88    90    92    96    98   100   112   114s=zeros(2,length(t1);for i=1:length(t1)   &

43、#160; s(1,i)=t1(i);s(2,i)=length(find(t=t1(i);enddisp(s)                           %统计间隔重复的次数 Columns 1 through 7        

44、0;  1           2           4           6           8       

45、   10          12           1        8169        8143       13549     

46、   5569        7079        8005  Columns 8 through 14          14          16        

47、  18          20          22          24          26        4233  

48、0;     2881        4909        2401        2172        2682        1175  Columns 15 through 21 

49、         28          30          32          34          36   

50、60;      38          40        1234        1914         550         557 

51、60;       767         330         424  Columns 22 through 28          42          44  

52、;        46          48          50          52          54    &#

53、160;    476         202         155         196         106          77 

54、60;       140  Columns 29 through 35          56          58          60          62&

55、#160;         64          66          68          53          54  

56、0;       96          16          24          48          13  Columns 36 through 42

57、60;         70          72          74          76          78   

58、       80          82          22          13          12     

59、60;     6          13           3           5  Columns 43 through 49          84 

60、;         86          88          90          92          96   &#

61、160;      98           6           4           1           4  

62、60;        1           2           1  Columns 50 through 52         100        

63、112         114           2           1           1max(t1)        

64、0;          %求最大间隔值ans =114    间隔规律:百万以内相邻素数间隔值有52个,其中间隔值2,4,6,8,10，12重复的次数较多,最大间隔值为114;另外10000以内最大间隔值为36,100000以内最大间隔值为72，所以随着整数范围的扩大,最大间隔值也随着扩大。1、 若在构造Koch曲线的过程中将向量CE绕点C逆时针旋转90度，并作出迭代三次的分形图。function q=koch(p)     

65、    q=;         t=90*pi/180;         M=cos(t),-sin(t);sin(t),cos(t);         for i=1:length(p)-1           

66、60; A=p(:,i);B=p(:,i+1);             C=A/3*2+B/3;             E=A/3+B/3*2;             D=C+M*(E-C);   

67、60;         q=q,A,C,D,E,B;         end p=0,1;0,0; q=koch(koch(koch(p); plot(q(1,:),q(2,:) axis(0 1 0 0.6) title(迭代三次的koch曲线)2、修改Sierpinski三角形的生成元，使其不使用中点而用一个三等份点，黑色的三角形调整为随机颜色的三角形，并作出迭代四次的分形图。fun

68、ction q=sierpinsk(p)q=;for i=1:3:length(p)    A=p(:,i);B=p(:,i+1);C=p(:,i+2);    D=A/3*2+B/3;E=B/3*2+C/3;F=C/3*2+A/3;    q=q,A,D,F,B,E,D,C,F,E;endfunction viewsierpinsk(p)hold onfor i=1:3:length(p)    fill(p(1,i:i+2),p(2,i:i+2),rand();endho

69、ld off clf pol=-1,1,0;0,0,sqrt(3); q=sierpinsk(sierpinsk(sierpinsk(sierpinsk(pol); viewsierpinsk(q)3、参考图10-4，分析Minkowwski“香肠”的生成元，并作出迭代五次的分形图。function q=minkowwsk(p)q=;t=90*pi/180;M=cos(t),-sin(t);sin(t),cos(t);N=cos(-t),-sin(-t);sin(-t),cos(-t);for i=1:length(p)-1  

70、60; A=p(:,i);B=p(:,i+1);    C=A/4*3+B/4;    E=(A+B)/2;    G=A/4+B/4*3;    D=C+M*(E-C);    F=E+N*(G-E);    H=E+N*(C-E);    J=G+M*(E-G);    q=q,A,C,D,H,E,F,J,G,B;endp=0,1;0,0; q

71、=minkowwsk(minkowwsk(minkowwsk(minkowwsk(minkowwsk(p);plot(q(1,:),q(2,:)2.对于logistic映射，选取适当的a，使迭代序列进入3,4,5,6周期，并给出周期轨道所用函数：function y=logistic(a,x0,n)f=(x)a*x*(1-x);y=;for i=1:n    y=y,x0;    x0=f(x0);end x=;y=;for a=0:0.02:4    x0=0.2;f=(x)a*x*(1-x)

72、;    for i=1:50        x0=f(x0);    end    for i=1:50        x0=f(x0);    end    for i=1:100        x0=f(x0);x=x,a;y=y,x

73、0;    endendplot(x,y,.)所用方法：首先用logistic函数来生成迭代序列，其次构造函数生成feigenbaum图，然后通过调整a的取值范围来观察图中周期分布并取近似值并一一试行。所得结果：logistic(3.84,0.02,100)（即a=3.84可使迭代序列进入3周期）                  周期轨道： 0.4880  

74、0; 0.9595    0.1494                  logistic(3.46,0.02,100)（即a=3.46可使迭代序列进入4周期）                  周期轨道： 0.8389 

75、;   0.4675    0.8613    0.4132                  logistic(3.74,0.02,100)（即a=3.74可使迭代序列进入5周期）             &

76、#160;    周期轨道： 0.6572    0.8425    0.4962    0.9349    0.2275                  logistic(3.628,0.02,100)（即a=3.628可使迭代序列进入6周期）  

77、;                周期轨道： 0.7705    0.6415    0.8344    0.5014    0.9070    0.30602、对于1000之内的n，求Mersenne数Mn=2n-1是素数的最大的n及对应的Mersenne素数的位数。只给出结果对于1000之内的n，Mersenne数Mn=2n-1是素数的最大的n 是607；对应的Mersenne素数的位数是183。1、已知采用密钥为5的加法加密方案的密文为 N fr f xyzijsy!，求明文。function dd=jf(ss,n)dd=ss-n;k=find(isletter(ss);dd(k)=ss(k);k=find(ss>=a&dd<a);dd(k)