数学实验作业汇总终审稿)

数学实验作业汇总

文稿归稿存档编号：[KKUY-KKIO69-OTM243-OLUI129-G00I-FDQS58-

（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)

（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）不用find函数，将向量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)=0

2、写出完成下列操作的命令及结果。

（1）将1~50这50个整数按行优先存放到5*10的矩阵中，求该矩阵四周元素的和；

>> t=[1:10];

>>M=[t;t+10;t+20;t+30;t+40]

M =

1 2 3 45 67 8 9 10

11 12 13 14 15 16 17 18 19 20

21 22 23 24 25 26 27 28 29 30

31 32 33 34 35 36 37 38 39 40

41 42 43 44 45 46 47 48 49 50

>>N=M(2:4,2:9)

N=

12 13 14 15 16 17 18 19

22 23 24 25 26 27 28 29

32 33 34 35 36 37 38 39

>> sum(sum(M))-sum(sum(n))

ans =

663

2）n取100、1000、10000，求序列1、1/2、1/3……1/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/241

>>?n=10000;

>> t=[1:n];>> format rat

>> M=t.^-1;

>> S=sum(M)

S =

1106/113

1.在同一坐标系下绘制y1=sin(t),y2=sin(2t),y3=sin(3t),其中y1的数据点用星号，线

形为黑色虚线，y2的数据点用方块，线形为红色实线，y3的数据点用小圆圈，线形为蓝色点线。（要求采用一次绘出和逐次填加两种方式完成绘图）

>> t=linspace(0,2*pi,100);

>> y1=sin(t);

>> y2=sin(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 off

2.分别用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曲线自动取值，在函数值变化平稳时，它的数值点会自动相对稀疏一点，在函数值变化剧烈处，所取点会自动密集一点，所以曲线更加光滑准确。

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);

>> subplot(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个三位整数，将其

中小于平均值的数用0代替。

function [mean,x]=ff (n) x=floor (100+899*rand (1,n)); m=length (x);

mean=sum (x)/m; x (x

4.编写函数，用来求下列函数的和，并给出n分别为100,1000,10000时，下列各式的

值。

function y=s(n)

y=1;

for i=1:1:n

x=4*i^2/(4*i^2-1);

y=y*x;

end

disp(y)

s(100)=1.5669

s(1000)=1.5704

s(10000)=1.5708

5.通过命令文件实现：随机产生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));

end

3.求下列函数的一阶和二阶导数

>> syms x>> diff(2/tan(x)+cos(x)/3,’x’,1)ans =- sin(x)/3 - (2*(tan(x)^2 + 1))/tan(x)^2

>> syms x?diff(2/tan(x)+cos(x)/3,’x’,2)

4.求积分

>> syms x

int(sqrt(exp(x)+1),x)

ans =

2*(exp(x) + 1)^(1/2) + 2*atan((exp(x) + 1)^(1/2)*i)*i

5.求下列级数的和

>> 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 =?x^4/24 + x^2/2 + 1

1.利用randn函数声称符合正态分布的10*5随机矩阵A，进行以下操作：

(1).A的各列元素的均值和标准方差

(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,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(ima_data));

3.下表所示是0~90度内某些数的正弦近似值

利用线性、样条差值求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,’linear’);

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.已知某次实验测得数据如下：

（1）请用3次多项式进行拟合，并给出拟合函数在0、0.5、1、1.5^9、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.25 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.381.61-1.23 0.05 5.62 15.65 30.3249.8074.28 103.92 138.91179.41225.61 277.67 335.79 400.12 470.85]

(2) 19阶拟合效果最好。理由通过编写差方和函数(基于最小二乘原理）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); 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）将1~50按列优先存放到5*10的矩阵M中：

（18）求矩阵M最大值所在的位置：

（19）统计字符串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

(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=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=A \b 运用矩阵分解

>> A=[2,2,-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\(Q\b)

4.在区间[30,50]内，求3()5sin()2log 1.8f x x x =-+ 的零点。 >> f=’5*sin(x)-2*(log(x)/log(3))+1.8’; >> ezplot(f,30,50) >> fzero(f,33) ans =

>> fzero(f,34)

ans =

33.3960

>> fzero(f,38)

ans =

39.0426

>> fzero(f,[39.4，39.5]) ans =

39.4785

则方程有四个零点

6.给出实验数据如下：

试分别用

b

x

b

y ae y a

x

==+

和做拟合形式，求出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. 90,10.75];>> X=1./x;>> Y=log(y);>> P=polyfit(X,Y,1)

P =

-1.1552 2.4629

>> exp(2.4629)

ans =

则a=11.7388 b=-1.1552 作图：>> Y1=polyval(P,X)

>> y1=exp(Y1);>> plot(x,y,’:o’,x,y1,’-*’) >> 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.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,y,’:o’,x,y1,’-*’)

3.求下列方程或方程的根在指定点的近似根

23

sin()ln 703210

50y x y z x z x y z ?++-=?+-+=??++-=?

，初值0001,1,1x y z === function f=myFun(x)f(1)=sin(x(1))+x(2)^2+log(x(3))-7;f(2)=3*x(1)+2^x(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.0050

2.

已知2sin cos 2(02)y x x x π=+≤≤ ，求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)^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.2222

1、 请你构造一个生成素数的公式,并将你的工作与

Euler 的工作比较。

采用素数生成公式p=n^2-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=n^2-79*n+1601系数y1=[-79 1601 100],y2=[-79 1601 1000]比较

得到结果

f(x1)=0.8614;f(x2)=0.5814;

f(y1)=0.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 2 4 6 8 10 12 14 16 18 20 22 24 26

Columns 15 through 28

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

Columns 29 through 42

56 58 60 62 64 66 68 70 72 74 76 78 80 82

Columns 43 through 52

84 86 88 90 92 96 98 100 112 114

s=zeros(2,length(t1));

for i=1:length(t1)

s(1,i)=t1(i);s(2,i)=length(find(t==t1(i)));

end

disp(s) %统计间隔重复的次数

Columns 1 through 7

1 2 4 6 8 10 12 1 8169 8143 13549 5569 7079 8005

Columns 8 through 14

14 16 18 20 22 24 26 4233 2881 4909 2401 2172 2682 1175 Columns 15 through 21

28 30 32 34 36 38 40 1234 1914 550 557 767 330 424

Columns 22 through 28

42 44 46 48 50 52 54 476 202 155 196 106 77 140

Columns 29 through 35

56 58 60 62 64 66 68 53 54 96 16 24 48 13

Columns 36 through 42

70 72 74 76 78 80 82 22 13 12 6 13 3 5

Columns 43 through 49

84 86 88 90 92 96 98 6 4 1 4 1 2 1

Columns 50 through 52

100 112 114 2 1 1

max(t1)%求最大间隔值

ans =114

间隔规律:百万以内相邻素数间隔值有52个,其中间隔值2,4,6,8,10，12重复的次数较多,最大间隔值为114;另外10000以内最大间隔值为36,100000以内最大间隔值为72，所以随着整数范围的扩大,最大间隔值也随着扩大。

1、若在构造Koch曲线的过程中将向量CE绕点C逆时针旋转90度，并作出迭代三次的分

形图。

function q=koch(p) q=[]; t=90*pi/180; M=[cos(t),-sin(t);sin(t),cos(t)]; for

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

q=[q,A,C,D,E,B];

endp=[0,1;0,0];q=koch(koch(koch(p)));plot(q(1,:),q(2,:))axis([0 1 0

0.6])title(’迭代三次的koch曲线’)

2、修改Sierpinski三角形的生成元，使其不使用中点而用一个三等份点，黑色的三角形调整为随机颜色的三角形，并作出迭代四次的分形图。

function 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];end

function viewsierpinsk(p)hold onfor i=1:3:length(p)

fill(p(1,i:i+2),p(2,i:i+2),rand());

endhold off

clfpol=[-

1,1,0;0,0,sqrt(3)];q=sierpinsk(sierpinsk(sierpinsk(sierpinsk(pol))));viewsier pinsk(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 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];end

p=[0,1;0,0];q=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);

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,x0]; endendplot(x,y,’.’)

所用方法：首先用logistic函数来生成迭代序列，其次构造函数生成feigenbaum图，然后通过调整a的取值范围来观察图中周期分布并取近似值并一一试行。

所得结果：logistic(3.84,0.02,100)（即a=3.84可使迭代序列进入3周期）

周期轨道：0.4880 0.9595 0.1494

logistic(3.46,0.02,100)（即a=3.46可使迭代序列进入4周期）

周期轨道：0.8389 0.4675 0.8613 0.4132

logistic(3.74,0.02,100)（即a=3.74可使迭代序列进入5周期）

周期轨道：0.6572 0.8425 0.4962 0.9349 0.2275

logistic(3.628,0.02,100)（即a=3.628可使迭代序列进入6周期）

周期轨道：0.7705 0.6415 0.8344 0.5014 0.9070 0.3060

=2n-1是素数的最大的n及对应的Mersenne素数2、对于1000之内的n，求Mersenne数M

n

的位数。只给出结果

对于1000之内的n，Mersenne数M

=2n-1是素数的最大的n 是607；对应的Mersenne素

n

数的位数是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)=dd(k)+26 ;k=find(ss>=’A’&dd<’A’);dd(k)=ss(k)+26;dd=char(dd);

步骤：先在M文件创建jf.m文件，然后在matlab程序中输入?dd=jf(’N fr f

xyzijsy!’,5)

结果：dd =I am a student!

2、已知采用密钥为“good”维吉尼亚加密方案的密文为 Nck gfs eci!，求明文。function dd=wjf(ss,key)

m=length(key);

for i=1:length(ss) c=mod(i,m); if c==0 c=m; end dd(i)=jf(ss(i),key(c)-

’a’);end

dd=wjf(’Nck gfs eci!’,’good’)

结果：ans =How are you!

3、A收到与之有秘密通信往来的B的一个密文信息,密文内容:

WOWUYSBACPGZSAVCOVKPEWCPADKPPABUJCQLYXQEZAACPP

按照双方的约定,采用Hill密码,密钥为 a={{1,2},{0,3}}，A~Z与整数0~25对应如下:A~Y对应1~25,Z对应0

求其原文

步骤：首先由A=[1，2；0，3]可得|A|=3，其次由命题条件可知3的逆为9

然后在Matlab程序中输入C=mod(9*[3,-2;0,1],26)，得A逆矩阵C=[1,8;0,9]

再者,输入A=[1,2;0,3];

B=[23,23,25,2,3,7,19,22,15,11,5,3,1,11,16,2,10,17,25,17,0,1,16;

15,21,19,1,16,0,1,3,22,16,23,16,4,16,1,21,3,12,24,5,1,3,16];

mod((C*B),26)

结果:ans =

Columns 1 through 14

13 9 21 10 1 7 1 20 9 9 7 1 7 9 5 7 15 9 14 0 9 1 16 14 25 14 10 14

Columns 15 through 23

24 14 8 9 9 5 8 25 14 9 7 1 4 8 19 9 1 14

最后根据明文字母表可得出原文为MEIGUOJIANGZAITAIPINGYANGJINXINGHAIDIHESHIYANN