MATLAB实验程序

实验一
要求： 查看matlab 的版本
代码：version
结果：ans =
7.1.0.246 (R14) Service Pack 3
1.先求下列表达式的值，然后显示MATLAB 工作空间的使用情况并保存全部变量。

（1）z 1=2185sin 2e
+
z1=(2*sin(85*pi/180))/(1+exp(2))
z =
0.2375
（2）z 2=⎥⎦⎤⎢⎣⎡+=++545
.0-i 212其中),1ln(212x x x x=[2,1+2*i;-0.45,5]
x =
2.0000 1.0000 + 2.0000i
-0.4500 5.0000
>> z2=1/2*log(x+sqrt(1+x*x))
z2 =
0.7114 - 0.0253i 0.8968 + 0.3658i
0.2139 + 0.9343i 1.1541 - 0.0044i
（3）z 3=2
3.0ln )3.0sin(23.03.0a a e e a
a +++--，a=-3.0,-2.9,-2.8, …3.0 提示：用冒号表达式生成a 向量，求各点的函数值时用点乘计算。

a=[-3.0:0.1:3.0]
z3=1/2.*((exp(0.3.*a)-exp(-0.3.*a))).*sin(a+0.3)+log((0.3+a).*1/2)
（4）z 4=⎪⎩⎪⎨⎧<≤+-<≤-<≤32,12
,21,110,222t t t t t t t 其中t=0:0.5:2.5
提示：用逻辑表达式求分段函数值。

t=0:0.5:2.5
t =
0 0.5000 1.0000 1.5000 2.0000 2.5000
>> z4=t.^2.*(t>=0&t<1)+(t.^2-1).*(t>=1&t<2)+(t.^2-2*t+1).*(t>=2&t<3) z4 =
0 0.2500 0 1.2500 1.0000 2.2500
2.已知：A=
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
7
65
3
87
7
34
4
-
34
12
，B=
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎣
⎡
7
2
-
3
3
2
1
-
3
1
求下列表达式的值：
（1）A+6*B和A-B+I（其中I为单位矩阵） A=[12,34,-4;34,7,87;3,65,7]
A =
12 34 -4
34 7 87
3 65 7
>> B=[1,3,-1;2,0,3;3,-2,7]
B =
1 3 -1
2 0 3
3 -2 7
>> A+6*B
ans =
18 52 -10
46 7 105
21 53 49
>> I=eye(size(A))
I =
1 0 0
0 1 0
0 0 1
>> A-B+I
ans =
12 31 -3
32 8 84
0 67 1
（2）A*B和A.*B
A*B
ans =
68 44 62
309 -72 596
154 -5 241
>> A.*B
ans =
12 102 4
68 0 261
9 -130 49
（3）A^3和A.^3
A^3
ans =
37226 233824 48604
247370 149188 600766
78688 454142 118820
>> A.^3
ans =
1728 39304 -64
39304 343 658503
27 274625 343
（4）A/B和B\A
A/B
ans =
16.4000 -13.6000 7.6000
35.8000 -76.2000 50.2000
67.0000 -134.0000 68.0000
>> B/A
ans =
0.1027 -0.0062 -0.0069
0.0617 0.0403 -0.0366
0.0205 0.0855 -0.0507
（5）[A,B]和[A([1,3],:);B^2]
[A,B]
ans =
12 34 -4 1 3 -1
34 7 87 2 0 3
3 65 7 3 -2 7
>> [A([1,3],:);B^2]
ans =
12 34 -4
3 65 7
4 5 1
11 0 19
20 -5 40
3.设有矩阵A和B
A=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡2524
2322212019181716151413121110987654321，B=⎥⎥⎥⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡111340794-23096-171603 （1） 求它们的乘积C 。

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]
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 6;17 -6 9;0 23 -4;9 7 0; 4 13 11]
B =
3 0 6
17 -6 9
0 23 -4
9 7 0
4 13 11
C=A*B
C =
93 150 67
258 335 177
423 520 287
588 705 397
753 890 507
（2） 将矩阵C 的左下角3×2子矩阵赋给D 。

(2)D=C([3,4,5],2:end)
D =
520 287
705 397
890 507
（3） 查看MATLAB 工作空间的使用情况。

4.完成下列操作：
（1）求[100，999]之间能被21整除的数的个数。

提示：先利用冒号表达式，再利用find 和length 函数。

A=100:999
P=rem(A,21)==0
Q=find(P)
>> length(Q)
ans =
43
（2）建立一个字符串向量，删除其中的大写字母。

提示：利用find 函数和空矩阵 。

ch='MmAaTtLlAaBb'
ch =
MmAaTtLlAaBb
>> k=find(ch>='A'&ch<='Z')
k =
1 3 5 7 9 11
ch(k)=[]
ch =
matlab
实验二、
1. 设有分块矩阵A=⎥⎦
⎤⎢⎣⎡⨯⨯⨯⨯23322333E S O R ，其中E 、R 、O 、S 分别为单位矩阵、随机矩阵、零矩阵和对角矩阵，试通过数值计算验证A=⎥⎦
⎤⎢⎣⎡+2E S O RS R 。

E=eye(3), R=rand([3,2]),O=zeros(2,3),S=diag([1,1]);A=[E,R;O,S]
E =
1 0 0
0 1 0
0 0 1
R =
0.9218 0.4057
0.7382 0.9355
0.1763 0.9169
O =
0 0 0
0 0 0
A =
1.0000 0 0 0.9218 0.4057
0 1.0000 0 0.7382 0.9355
0 0 1.0000 0.1763 0.9169
0 0 0 1.0000 0
0 0 0 0 1.0000
2.产生5阶希尔伯特矩阵H和5阶帕斯卡矩阵P，且求其行列式的值Hh
和Hp以及它们的条件数Th和Tp，判断哪个矩阵的性能更好。

为什么？H=hilb(5),P=pascal(5)
H =
1.0000 0.5000 0.3333 0.2500 0.2000
0.5000 0.3333 0.2500 0.2000 0.1667
0.3333 0.2500 0.2000 0.1667 0.1429
0.2500 0.2000 0.1667 0.1429 0.1250
0.2000 0.1667 0.1429 0.1250 0.1111
P =
1 1 1 1 1
1 2 3 4 5
1 3 6 10 15
1 4 10 20 35
1 5 15 35 70
>> Hh=det(H),Pp=det(P)
Hh =
3.7493e-012
Pp =
1
>> Th=cond(H),Tp=cond(P)
Th =
4.7661e+005
Tp =
8.5175e+003
矩阵T 条件数比矩阵H 的条件数更接近1，因此，矩阵T 的性能要好于矩阵H 。

3. 建立一个5×5矩阵，求它的行列式值、迹、秩和范数。

A=[1 4 7 10 10;2 5 8 11 7;3 6 9 4 8;11 12 3 2 1;5 6 9 8 7]
A =
1 4 7 10 10
2 5 8 11 7
3 6 9
4 8
11 12 3 2 1
5 6 9 8 7
>> a1=det(A),a2=trace(A),a3=rank(A),a4=norm(A)
a1 =
5376
a2 =
24
a3 =
5
a4 =
32.3947
4. 已知A=⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡588-1252018629-
求A 的特征值及特征向量，并分析其数学意义。

A=[-29 6 18;20 5 12;-8 8 5]
A =
-29 6 18
20 5 12
-8 8 5
>> [V,D]=eig(A)
V =
0.7130 0.2803 0.2733
-0.6084 -0.7867 0.8725
0.3487 0.5501 0.4050
D =
-25.3169 0 0
0 -10.5182 0
0 0 16.8351
数学意义：求得的3个特征值是-25.3169、-10.5182、16.8351，各特征值对应的特征向量为V 的各列构成的向量。

5.
下面是一个线性方程组：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡52.067.095.06/15/14/15/14/13/14/13/12/1321x x x (1) 求方程的解。

A=[1/2 1/3 1/4;1/3 1/4 1/5;1/4 1/5 1/6]
A =
0.5000 0.3333 0.2500
0.3333 0.2500 0.2000
0.2500 0.2000 0.1667
B=[0.95;0.67;0.52]
B =
0.9500
0.6700
0.5200
x=inv(A)*B
x =
1.2000
0.6000
0.6000
(2)
将方程右边向量元素b 3改为0.53，再求解，并比较b 3的变化和解的相对
变化。

B=[0.95;0.67;0.53]
B =
0.9500
0.6700
0.5300
>> x=inv(A)*B
x =
3.0000
-6.6000
6.6000
（4） 计算系数矩阵A 的条件数并分析结论。

cond(A)
ans =
1.3533e+003
6.建立A 矩阵，试比较sqrtm(A)和sqrt(A)，并分析它们的区别。

A=rand(3)
A =
0.4103 0.3529 0.1389
0.8936 0.8132 0.2028
0.0579 0.0099 0.1987
>> sqrtm(A)
ans =
0.4248 0.2963 0.1314
0.7418 0.7709 0.0873
0.0764 -0.0106 0.4354
>> sqrt(A)
ans =
0.6405 0.5940 0.3727
0.9453 0.9018 0.4503
0.2406 0.0993 0.4458
实验三、
1. 求分段函数的值。

其他503且0,1,65,6222<≤-≠<⎪⎩
⎪⎨⎧--+--+=x x x x x x x x x y
用if 语句实现，分别输出x=-5.0,-3.0,1.0,2.0,2.5,3.0,5.0时的y 值。

2. 输入一个百分制成绩，要求输出成绩等级A 、B 、C 、D 、E 。

其中90分～100分为A ，80分～89分为B ，70分～79分为C ，60分～69分为D ，60分以下为E 。

要求：
（1） 分别用if 语句和switch 语句实现。

（2） 输入百分制成绩后要判断该成绩的合理性，对不合理的成绩应输出出
错信息。

用IF 语句：x=input('请输入X 的值：')
if x<=100&x>=90
disp('A');
elseif x>=80&x<=89
disp('B');
elseif x>=70&x<=79
disp('C');
elseif x>=60&x<=69
disp('D');
elseif x>=0&x<=59
disp('E');
else x<0&x>100
disp('输入有误');
end
运行结果：请输入X 的值：97
x =
97
A
请输入X的值：132
x =
132
ans =
输入有误:
用swich语句：x=input('请输入x的值：');
swich floor(x/10)
case(9)
disp('A');
case(8)
disp('B');
case(7)
disp('C');
case(6)
disp('D');
case num2cell(0:5)
disp('E');
otherwise
disp('输入有误：');
end
3.硅谷公司员工的工资计算方法方法如下：
（1）工作时数超过120小时者，超过部分加发15%。

（2）工作时数低于60小时者，扣发700元。

（3）其余按每小时84小时计发。

试编程按输入的工号和该号员工的工时数，计算应付工资。

n=input('请输入工号：')
x=input('请输入该员工的工时数：')
if x>=120
y=84.*x+(x-120)*84*0.15;
elseif x<120&x>=60
y=84*x;
else x<60
y=84*x-700;
end
disp([n,'请输入该员工的工资：',x]);
y
请输入工号：001
n =
1
请输入该员工的工时数：100
x =
100
请输入员工的工资：d
y =
8400
请输入工号：002
n =
2
请输入该员工的工时数：59
x =
59
请输入员工的工资：;
y =
4256
请输入工号：003
n =
3
请输入该员工的工时数：122
x =
122
请输入员工的工资：z
y =
1.0273e+004
4.设计程序，完成两位数的加、减、乘、除四则运算，即产生两个两位随机整
数，再输入一个运算符号，做相应的运算，并显示相应的结果。

程序如下：
clear all
close all
clc
a=input('请输入第一个数¸öÊý');
b=input('请输入第二数¸öÊý');
c=input('请输入第运算符,'s');
if c=='+'
s=a+b;
else
if c=='-'
s=a-b;
else
if c=='*'
s=a*b;
else
if c=='/'
s=a/b;
else
disp('error');
end
end
end
end
disp(s);
请输入第一个数68
请输入第二5
请输入运算符/
1.0462
请输入第一个数2
请输入第二个数3
请输入运算符-
-1
>>
请输入第一个数2
请输入第二个数3
请输入运算符*
6
请输入第一个数2
请输入第二个数3
请输入运算符+
5
5.建立5×6矩阵，要求输出矩阵第n行元素。

当n值超过矩阵的行数时，自
动转为输出矩阵最后一行元素，并给出出错信息。

程序如下：
clear all
close all
clc
A=1:30;
B=reshape(A,5,6);
n=input('input');
try
disp(B(n,:))
catch
disp(B(5,:))
end
lasterr
结果：
input4
4 9 14 19 24 29
input8
5 10 15 20 25 30
Attempted to access B(8,:); index out of bounds because size(B)=[5,6].
实验四、
例题
例题4.7
程序如下：
for m=100:999
m1=fix(m/100);
m2=rem(fix(m/10),10);
m3=rem(m,10);
if m==m1*m1*m1+m2*m2*m2+m3*m3*m3
disp(m)
end
end
输出结果为：
153
370
371
407
例4.8
程序如下：
y=0;n=100;
for i=1:n
y=y+1/i/i;
end
y
输出结果为：
y =
1.6350
例4.9
程序如下：
a=0;b=3*pi;
n=1000;h=(b-a)/n;
x=a;s=0;
f0=exp(-0.5*x)*sin(x+pi/6);
for i=1:n
x=x+h;
f1=exp(-0.5*x)*sin(x+pi/6);
s=s+(f0+f1)*h/2;
f0=f1;
end
s
运行结果：
0.9008
方法二：
a=0;b=3*pi;
n=1000;h=(b-a)/n;
x=a:h:b;
f=exp(-0.5*x)*sin(x+pi/6);
for i=1:n
s(i)=(f(i)+f(i+1))*h/2;
end
s=sum(s)
例4.10
程序如下：
s=0;
a=[12,13,14;15,16,17;18,19,20;21,22,23]; for k=a
s=s+k;
end
disp(s);
运行结果：
39
48
57
66
例4.11
程序如下：
sum=0;
n=0;
x=input('Enter a number(end in 0):'); while x~=0
sum=sum+x;
n=n+1;
x=input('Enter a number(end in 0):'); end
if(n>0)
sum
mean=sum/n
end
运行结果：
Enter a number(end in 0):67
Enter a number(end in 0):89
Enter a number(end in 0):93
Enter a number(end in 0):70
Enter a number(end in 0):0
319
mean =
79.7500
例4.12
程序如下：
X=input('Enter X:');
E=zeros(size(X));
F=eye(size(X));
n=1;
while norm(F,1)>0
E=E+F;
F=F*X/n;
n=n+1;
end
E
expm(X)
结果如下：
Enter X:[0.5,2,0;1,-1,-0.5;0.9,1,0.75]
E =
2.6126 2.0579 -0.6376
0.7420 0.7504 -0.5942
2.5678 2.3359 1.5549
ans =
2.6126 2.0579 -0.6376
0.7420 0.7504 -0.5942
2.5678 2.3359 1.5549
例4.13
程序如下：
for n=100:200
if rem(n,21)~=0
continue
end
break
end
n
运行结果：
n =
105
例4.14
程序如下：
for m=1:500
s=0;
for k=1:m/2
if rem(m,k)==0
s=s+k;
end
end
if m==s
disp(m);
end
end
运行结果：
n =
105
6
28
496
例4.15
程序如下：
m=input('m=');
p=1:m;
p(1)=0;
for i=2:sqrt(m)
for j=2*i:i:m
p(j)=0;
end
end
n=find(p~=0);
p(n)
运行结果：
m=20
ans =
2 3 5 7 11 13 17 19
1.根据+=22
1
16π2221...3121n +++ ，求π的近似值。

当n 分别取100、1000、10000时，结果是多少？
要求：分别用循环结构和向量运算（使用sum 函数）来实现。

1. 程序如下：
n=input('Enter n:');
sum=0;
i=1;
while (i<n)
d=1/(i^2);
sum=sum+d;
i=i+1;
end
p=sqrt(6*sum);
disp(p)
disp(abs(p-pi))
0.0096
Enter n:1000
3.1406
9.5555e-004
运行结果：
Enter n:100
3.1320
0.0096
Enter n:1000
3.1406
9.5555e-004
Enter n:10000
3.1415
9.5499e-005
2.根据y=1+1
21...715131+++++n ，求： （1）y<3时的最大值。

（2）与（1）的n 值对应的y 值。

程序如下：
clear all
n=1;
y=0;
while (y<3)
x=1/(2*n-1);
n=1+n;
y=y+x;
end
disp(['y<3时的最大n 值是：',num2str(n-2)])
disp(['相应的y 值是：',num2str(y-x)])
运行结果：
y<3时的最大n 值是：56
相应的y 值是：2.9944
3.考虑以下迭代公式： n n x b a x +=
+1
其中a 、b 为正的常数。

（1） 编写程序求迭代的结果，迭代的终止条件是5110-+≤-n n x x ，迭代初值
0x =1.0，迭代次数不超过500次。

clear all
close all
clc
x=1;
a=1;
b=1;
n=1;
while (abs(a/(b+x)-x)>10e-5&n<=500)
运行结果：0.6180
10
（2） 如果迭代过程收敛于r ，那么r 的准确值是2
4-2a b b +±，当（a,b ）的值取（1,1）、（8,3）、（10,0.1）时，分别对迭代结果和准确值进行比较。

clear all
close all
clc
x=1;
a=1;
b=1;
n=1;
while(abs(a/(b+x)-x)>10e-5&n<=500)
x=a/(b+x);
n=n+1;
end
x
r1=(-b+sqrt(b^2+4*a))/2;
r2=(-b-sqrt(b^2+4*a))/2;
r1
r2
结果：x =
0.6180
r1 =
0.6180
r2 =
-1.6180
4.已知
⎪⎪⎩⎪⎪⎨⎧+-====---,2,1,0,13214321n n n f f f f f f f 3321>===n n n n
求1f ～100f 中：
（1） 最大值、最小值、各数之和。

（2） 正数、零、负数的个数。

5.若两个连续自然数的乘积减1是素数，则称这两个连续自然数是亲密数对，该素数是亲密素数。

例如，2×3-1=5，由于5是素数，所以2和3是亲密数对，5是亲密素数。

求[2,50]区间内：
（1） 亲密数对的对数。

（2） 与上述亲密数对对应的所有亲密素数之和。

clear all
close all
clc
p=0;
l=0;
H=[];
for m=2:50
a=m;
b=m-1;
c=a*b-1;
for k=1:c
if rem(c,k)==0
l=l+1;
end
end
if l==2
disp(['亲密对数是',num2str(a),',',num2str(b)])
p=p+1;
H=[H,c];
end
l=0;
end
sum=sum(H);
disp(['亲密对数个数是：',num2str(p)])
disp(['亲密对数和是：',num2str(sum)])
运行结果是：
亲密数对是3,2
亲密数对是4,3
亲密数对是5,4
亲密数对是6,5
亲密数对是7,6
亲密数对是9,8
亲密数对是10,9
亲密数对是11,10
亲密数对是12,11
亲密数对是14,13
亲密数对是16,15
亲密数对是17,16
亲密数对是20,19
亲密数对是21,20
亲密数对是22,21
亲密数对是25,24
亲密数对是27,26
亲密数对是29,28
亲密数对是31,30
亲密数对是32,31
亲密数对是36,35
亲密数对是39,38
亲密数对是40,39
亲密数对是42,41
亲密数对是45,44
亲密数对是46,45
亲密数对是47,46
亲密数对是49,48
亲密数对的个数是：28 亲密数对的和是：21066 >>
实验五、
例4.18
程序如下：
function f=factor(n) if n<=1
f=1;
else
f=n*factor(n-1) end
运行结果：>> factor(5) ans =
120
求s=1!+2!+3!+4!+5!
程序如下：
s=0;
for i=1:5
s=s+factor(i);
end
运行结果：
s =
153
例4.19
任意排列问题。

编写一个函数实现randperm（n）函数的功能。

function Y=rndprm1(X)
[m,n]=size(X)
if m>1
error('RNDPRM1 accepts as inputs only vecyors'); end
Y=[];
l=n;
for i=1:n
k=1+fix(l*rand);
x=X(k);
Y=[Y,x];
X(k)=[];
l=l-1;
end
rndprm1([34,6,3,54,2,5,454])
m =
1
n =
7
ans =
5 454 34 2 3
6 54
function Y=rndprm2(X)
[m,n]=size(X);
l=n;
if m>1
error('RNDPRM2 accepts as inputs only vectors')
end
if n<=1
Y=X;
else
k=1+fix(l*rand);
x=X(k);
X(k)=[];
Z=rndprm2(X);
Y=[Z,x];
l=l-1;
end
rndprm2('apple') ans = epapl >>
1. 定义一个函数文件，求给定复数的指数、对数、正弦和余弦，并在命令文件中调用该函数文件。

函数文件： complexnum.m
function [expval,logval,sinval,cosval]=complexnum(x) %COMPLEXNUM compute complex number value expval=exp(x); logval=log(x); sinval=sin(x); cosval=cos(x); 命令文件：exp0501. clear all; close all; clc;
x=input('请输入一个复数：');
[expval,logval,sinval,cosval]=complexnum(x); disp(['给定复数的指数是：', num2str(expval)]) disp(['给定复数的对数是：', num2str(logval)]) disp(['给定复数的正弦是：', num2str(sinval)]) disp(['给定复数的余弦是：', num2str(cosval)]) 运行结果：
请输入一个复数：1+2i
给定复数的指数是：-1.1312+2.4717i 给定复数的对数是：0.80472+1.1071i 给定复数的正弦是：3.1658+1.9596i 给定复数的余弦是：2.0327-3.0519i
2. 一物理系统可用下列方程组来表示：
⎥⎥⎥⎥⎦
⎤⎢⎢⎢
⎢⎣⎡=⎥
⎥
⎥⎥
⎦⎤
⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣
⎡----g m g m N N a a m m m m 2121212
111001cos 0
0sin 00cos 0
sin 0sin cos θ
θθθθθ 从键盘输入θ和、21m m 的值，求θ，输入8.9取的值。

其中和、2121g N N a a 时以角度为单位。

要求：定义一个求解线性方程组AX=B 的函数文件，然后在命令文件中调用该函数文件。

函数文件：equationsolution.m
function J= equationsolution( m1,m2,theta )
H=[m1*cos(theta*pi/180), -m1, -sin(theta*pi/180), 0;
m1*sin(theta*pi/180), 0, cos(theta*pi/180), 0;
0, m2, -sin(theta*pi/180), 0;
0, 0, -cos(theta*pi/180), 1];
K=[0;m1*9.8;0;m2*9.8];
J=inv(H)*K;
End
命令文件：exp0502.m
clear all;
close all;
clc; m1=input('m1=');
m2=input('m2=');
theta=input('theta=');
J=equationsolution(m1, m2, theta);
disp(['方程组的解alpha1', num2str(J(1))])
disp(['方程组的解alpha2', num2str(J(2))])
disp(['方程组的解N1', num2str(J(3))])
disp(['方程组的解N2', num2str(J(4))])
运行结果：
m1=2
m2=3
theta=45
方程组的解alpha1：8.6621
方程组的解alpha2：2.45
方程组的解N1：10.3945
方程组的解N2：36.75
3.一个自然数是素数，且它的数字位置经过任意对换后仍为素数，则称是绝对
素数。

例如12是绝对素数。

试求所有两位数素数。

要求：定义一个判断素数的函数文件。

函数文件：prime.m
function a= prime(b)
%absolute prime
x=fix(b/10);
y=rem(b,10);
a=10*y+x;
signab=0;
for m=2:sqrt(b)
if rem(b,m)==0
signab=signab+1;
end
end
for n=2:sqrt(a)
if rem(a,n)==0
signab=signab+1; end end
if signab==0 a=b; else
a=0; %ÕâÀ︳ֵ0£¬±íʾ²»ÊÇË×Êý end end
执行文件：exp0503.m clear all; close all; clc;
twoprime=[];
for n=10:99
[a]=prime(n); if a>0
twoprime=[twoprime,a]; end end
disp(['所有的两位绝对素数是', num2str(twoprime)]) 运行结果：
所有的两位绝对素数是：11 13 17 31 37 71 73 79 97 4. 设01
.0)3(1
1.0)2(1)(4
2+-++-=
x x x f ，编写一个MATLAB 函数文件fx.m ，使得调用f(x)时，x 可用矩阵代入，得出的f(x)为同阶矩阵。

函数文件：fx.m function f=fx(x)
f=1./((x-2).^2+0.1)+1./((x-3).^4+0.01); end
执行文件:exp0504.m clear all; close all; clc;
x=[1,2;3+4i,5] f=fx(x 运行结果: x =
1.0000
2.0000
3.0000 +
4.0000i
5.0000 f =
0.9716 10.9901 -0.0482 - 0.0280i 0.1724
5. 已知)
20()30()
40(f f f y +=
（1） 当)5ln(19)(2++=n n n f 时，求y 的值。

函数文件：fun5.m function y=fun5(n) y=n+10*log((n^2)+5); end
执行文件：exp0505.m clear all; close all; clc;
y=fun5(40)/(fun5(30)+fun5(20)) 运行结果： y =
0.6390
（2） 当f(n)=1×2+2×3+2×3+…+n×(n+1)时，y 的值。

函数文件：fun.m function y=fun(n) y=0;
for m=1:n
y=y+m*(m+1); end
执行文件：exp0506.m clear all; close all; clc;
y=fun(40)/(fun(30)+fun(20)) 运行结果：y = 1.7662 实验六、 例5.37
[x,camp]=imread('greens.jpg'); image(x);
colormap(cmap); axis imageoff
例5.38
[x,y,z]=sphere(50);
m=moviein(30);
for i=1:30
surf(i*x,i*y,i*z)
m(:,i)=getframe;
end
movie(m,10)
例5.39
n=30;
s=0.02;
x=rand(n,1)-0.5;
y=rand(n,1)-0.5;
h=plot(x,y,'.');
axis([-1 1 -1 1]);
axis square
grid off;
set(h,'EraseMode','Xor','Markersize',20); for i=linspace(1,10,5000)
drawnow
x=x+s*randn(n,1); y=y+s*randn(n,1);
set(h,'XData','x','YData',y); end
1. 设ocsx x x y ⎥⎦⎤
⎢⎣⎡++=21sin 35.0，在x=0～2π区间取101点，绘制函数的曲线。

程序如下：
>> x=0:2*pi/100:2*pi;
>> y=(0.5+3*sin(x)/(1+x.^2)).*cos(x); >> plot(x,y); 结果：
2.已知y1=
,y2=cos(2x)y3=y1*y2，完成下列操作
:
1
2
3
4
5
6
7
（1）在一坐标系下用不同的颜色和线型绘制三条曲线。

程序：>> x=linspace(-2*pi,2*pi,500);
>> y1=x.^2;
>> y2=cos(2*x);
>> y3=y1.*y2;
>> plot(x,y1,'k-',x,y2,'r:',y1,y2,'g--')
-10-50510152025303540
（2）以子图形式绘制条曲线。

程序：x=linspace(-2*pi,2*pi,500);
y1=x.^2;
y2=cos(2*x);
y3=y1.*y2;
subplot(2,2,1);
plot(x,y1);
title('y1=x.^2');
subplot(2,2,2);
plot(x,y2);
title('y2=cos(2*x)');
subplot(2,1,2);
plot(y1,y2);
title(' y3=y1.*y2');
（3） 分别用条形图，阶梯图，杆图和填充图绘制三条曲线。

程序：x=linspace(-2*pi,2*pi,500); y1=x.^2;
y2=cos(2*x); y3=y1.*y2;
subplot(3,4,1);bar(x,y1); title('bar(y1)');
subplot(3,4,2);stairs(x,y1); title('stairs(y1)');
subplot(3,4,3);stem(x,y1); title('stem(y1)');
subplot(3,4,4);fill(x,y1,'k'); title('fill(y1)');
subplot(3,4,5);bar(x,y2); title('bar(y2)');
subplot(3,4,6);stairs(x,y2); title('stairs(y2)');
subplot(3,4,7);stem(x,y2); title('stem(y2)');
subplot(3,4,8);fill(x,y2,'g'); title('fill(y2)');
subplot(3,4,9);bar(x,y3); title('bar(y3)');
subplot(3,4,10);stairs(x,y3); title('stairs(y3)');
subplot(3,4,11);stem(x,y3);
-10-50510
y1=x.2
-10
-50510
-1-0.500.51y2=cos(2*x)
0510152025303540
-0.50.5 y3=y1.*y2
title('stem(y3)');
subplot(3,4,12);fill(x,y3,'r'); title('fill(y3)') >>
3.已知：
，在区间绘制函数图像。

程序： clear all
x=linspace(-5,5,100); y=[]; for x0=x if x0<=0
y1=(x0+sqrt(pi))/exp(2); elseif x0>0
y1=0.5*log(x0+sqrt(1+x0.^2)); end
y=[y,y1]; end
>> plot(x,y)
-10010020
40bar(y1)-10
010020
40stairs(y1)
020
40stem(y1)
-10
010
02040fill(y1)
-10010-10
1bar(y2)
-10
010-101stairs(y2)
-10
1stem(y2)
-10
010
-101fill(y2)
-10
010-500
50bar(y3)
-10
010-50050stairs(y3)
-10
010-50050stem(y3)
-10
010
-50050fill(y3)
4.绘制极坐标曲线)sin(θρn b a +=，分析参数a ，b ，n 对曲线形状的影响。

程序：clear all
a=input('输入a 的值：'); b=input('输入b 的值：'); n=input('输入n 的值：'); theta=0:0.01:2*pi; rho=a*sin(b+n*theta); polar(theta,rho,'r')
参数a 、b 、n 对曲线形状的影响如下：
输入a 的值：34：输入b 的值：15：输入n 的值：23
-5
-4-3-2-1012345
-0.6-0.4-0.2
0.20.40.60.81.2
1）输入a 的值：1：输入b 的值：15：输入n 的值：23
2）输入a 的值：34：输入b 的值：1：输入n 的值：23
90270
180
90
270
180
3）输入a 的值：34：输入b 的值：15：输入n 的值：1
5.绘制函数的曲面图和等高线。

其中x 的21个值均匀
分布在[-5,5]范围，y 的31个值均匀分布在[0,10]，要求使用subplot （2,1,1）subplot(2,1,2)将产生的曲面图和等高线图在同一个窗口上。

程序：clear all
x1=linspace(-5,5,21);
90270
180
90270
180
y1=linspace(0,10,31);
[x,y]=meshgrid(x1,y1);
z=cos(x).*cos(y).*exp(-0.25.*sqrt(x.^2+y.^2));
subplot(2,1,1);
surf(x,y,z);
subplot(2,1,2);
contour3(x,y,z,30);
6.绘制曲面图形，并进行插值着色处理。

，
程序：
clc;
s=0:pi/100:pi/2;
t=0:pi/100:3*pi/2;
[s,t]=meshgrid(s,t);
x=cos(s).*cos(t);
y=cos(s).*sin(t);
z=sin(s);
subplot(2,2,1);
mesh(x,y,z);
title('未着色的图形');
subplot(2,2,2);
surf(x,y,z);
title('shading faceted (缺省)');
subplot(2,2,3);
surf(x,y,z);shading flat;
title('shading flat');
subplot(2,2,4);
surf(x,y,z);shading interp;
title('shading interp');
实验八、
1.利用MATLAB提供的rand函数生成30000个符合均匀分布的随机数，然后检
验随机数的性质：
（1）均值和标准方差。

（2）最大元素和最小元素。

（3）大于0.5的随机数个数占总数的百分比。

程序：x=rand(1,30000);
>> max(x)
ans =
0.9999
>> min(x)
ans =
3.5706e-005
>> mean(x)
ans =
0.5010
>> std(x)
ans =
0.2882
2.将100个学生5门功课的成绩存入矩阵P中，进行如下处理：
（1）分别求每门课的最高分、最低分及相关学生序号。

（2）分别求每门课的平均分和标准方差。

（3）5门课总分的最高分、最低分及相关学生的序号。

（4）将5门课的总分按从大到小的顺序存入zcj中，相应学生序号存入xsxh。

提示：上机调试时，为避免输入学生成绩的麻烦，可用取值范围在[45,95]之间的随机矩阵来表示学生成绩。

clc
t=45+50*rand(100,5);
P=fix(t);
[x,l]=max(P)
[y,k]=min(p)
mu=mean(P)
sig=std(P)
s=sum(P,2)
[X,m]=max(s)
[Y,n]=min(s)
[zcj,xsxh]=sort(s)
运行结果：
x =
94 94 94 94 94
l =
16 14 16 34 27
y =
0.5005
k =
1
mu =
69.5000 68.5900 72.3800 71.6900 71.0400
sig =
14.8157 13.9480 14.5490 14.8376 14.7798
X =
422
m =
42
Y =
266
n =
36
zcj =
266
……
422
xsxh =
36
……
42
3.某气象观测站测得某日6：00~18:00之间每隔2h的室外温度（摄氏度）如实验表1所示。

实验表1 室内外温度观测结果（摄氏度）
时间h 6 8 10 12 14 16 18
室外温度t1 18.0 20.0 22.0 25.0 30.0 28.0 24.0
室内温度t2 15.0 19.0 24.0 28.0 34.0 32.0 30.0
试用三次样条插值分别求出该日室内外6:00~17:30之间每隔2h各点的近似温度（摄氏度）。

clc
x=[6,8,10,12,14,16,18];
y1=[18.0,20.0,22.0,25.0,30.0,28.0,24.0];
y2=[15.0,19.0,24.0,28.0,34.0,32.0,30.0];
x0=[6.30:2.0:17.30];
y11=spline(x,y1,x0)
y12=spline(x,y2,x0)
运行结果：
y11 =
18.3014 20.2991 22.3055 25.7999 30.1996 27.2962
y12 =
15.3558 19.7984 24.5613 28.9490 34.2441 31.3618 4．已知lgx在[1,101]区间10个整数采样点的函数值如实验表2所示。

数曲线。

x=1:10:101;
y=log10(x);
p=polyfit(x,y,5)
p =
0.0000 -0.0000 0.0001 -0.0058 0.1537 -0.1326 y1=polyval(p,x)
plot(x,y,':o',x,y1,'-*') y1 =
Columns 1 through 8
0.0153 0.9868 1.3720 1.5103 1.5887 1.6813 1.7885 1.8771
Columns 9 through 11
1.9189 1.9308
2.0140
5
．
有
3
个
多
项
式
32)(,2)(,542)(2322341++===+++=x x x P x x P x x x x P ，试进行
下列操作：
（1）求)()()()(P 321x P x P x P x +=。

clc
p1=[1,2,4,0,5];p2=[0,1,2];p3=[1,2,3]; p4=conv(p2,p3); p=p1+p4 运行结果： p =
1 3 8 7 11 （2）求P （x ）的根。

x=roots(p) x =
-1.3840 + 1.8317i -1.3840 - 1.8317i -0.1160 + 1.4400i -0.1160 - 1.4400i
（3）当x 取矩阵A 的每一元素时，求P 得值。

其中：
⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡--=5.2505.3275.04.12.11A
clc
p1=[1,2,4,0,5];p2=[0,1,2];p3=[1,2,3]; A=[-1,1.2,-1.4;0.75,2,3.5;0,5,2.5]; y=polyval(p,A) 运行结果： y =
1.0e+003 *
0.0100 0.0382 0.0125 0.0223 0.0970 0.4122 0.0110 1.2460 0.1644 （5） 当以矩阵A 为自变量时，求P(x)的值。

其中A 的值与第（3）题相同。

y1=polyvalm(p,A) y1 =
1.0e+003 *
0.0076 -0.1281 -0.0775 0.1328 1.3900 1.1644 0.1824 1.7364 1.5198 实验九、
1求函数在指定点的数值导数。

3,2,1,620321)(232=⎥⎥⎥
⎦
⎤⎢⎢⎢
⎣⎡=x x x x
x x x x f clc
x=1; i=1;
f=inline('det([x,x^2,x^3;1,2*x,3*x^2;0,2,6*x])'); while x<=3.01 g(i)=f(x); i=i+1; x=x+0.01; end g;
t=1:0.01:3.01; dx=diff(g)/0.01;
f1=dx(1) f2=dx(101)
f3=dx(length(g)-1) 运行结果： f1 =
6.0602 f2 =
24.1202 f3 =
54.1802
2．用符号方法求定积分。

（1）⎰
++=
π
20
2211)2sin(4cos I dt t t 的近似值。

clc
x=sym('x');t=sym('t');
f=sqrt(cos(t.^2)+4*sin(2*t).^2+1); I1=int(f,t,0,2*pi) 运行结果： I1 =
int((cos(t^2) + 4*sin(2*t)^2 + 1)^(1/2), t = 0..2*pi)
（2）dx x
x ⎰
++=
1
2
21)
1ln(I 。

clc
x=sym('x');t=sym('t'); g=log(1+x)./(1+x.^2); I2=int(g,x,0,1) 运行结果： I2 =
(pi*log(2))/8 3．求函数的极限。

（1）4
58
6lim 24+-+-→x x x x x （2）x x x 0lim →（3）x x x 11lim 20-+→（4）x x x x )1(lim +∞→
(1)clear all
close all clc syms x;
f=(x^2-6*x+8)/(x^2-5*x+4); limit(f,4) 运行结果： ans = 2/3
（2）clc
syms x;
f=abs(x)/x;
limit(f)
运行结果：
ans =
NaN
（3）
f=(sqrt(1+x^2)-1)/x; limit(f)
ans =
(4)f=(x+1/x)^x;
limit(f,x,inf,'left') ans =
Inf
f=(x+1/x)^x;
limit(f,x,inf,'right') ans =
Inf。