MATLAB实验练习题计算机南邮MATLAB数学实验大作业答案.docx

要求：抄题、写出操作命令、运行结果，并根据要求，贴上运行图。
1、求的所有根。（先画图后求解）（要求贴图）
>> solve('exp(x)-3*x^2',0)
ans =
-2*lambertw(-1/6*3^(1/2))
-2*lambertw(-1,-1/6*3^(1/2))
-2*lambertw(1/6*3^(1/2))

2、求下列方程的根。
1)
a=solve('x^5+5*x+1',0);a=vpa(a,6)
a =
1.10447+1.05983*i
-1.00450+1.06095*i
-.199936
-1.00450-1.06095*i
1.10447-1.05983*i
2）至少三个根

>> fzero('x*sin(x)-1/2', 3)
ans =
2.9726
>> fzero('x*sin(x)-1/2',-3)
ans =
-2.9726
>> fzero('x*sin(x)-1/2',0)
ans =
-0.7408
3）所有根

>> fzero('sin(x)*cos(x)-x^2',0)
ans =
0
>> fzero('sin(x)*cos(x)-x^2',0.6)
ans =
0.7022
3、求解下列各题：
1）
>> sym x;
>> limit((x-sin(x))/x^3)
ans =
1/6
2）
>> sym x;
>> diff(exp(x)*cos(x),10)
ans =
(-32)*exp(x)*sin(x)
3）
>> sym x;
>> vpa((int(exp(x^2),x,0,1/2)),17)
ans =
0.54498710418362222
4）
>> sym x;
>> int(x^4/(25+x^2),x)
ans =
125*atan(x/5) - 25*x + x^3/3
5）求由参数方程所确定的函数的一阶导数与二阶导数。
>> sym t;
>> x=log(sqrt(1+t^2));y=atan(t);
>> diff(y,t)/diff(x,t)
ans =
1/t
6）设函数y=f(x)由方程xy +ey=e所确定，求y′(x)。
>> syms x y;
f=x*y+exp(y)-exp(1);
>> -diff(f,x)/diff(f,y)
ans =
-y/(x + exp(y))
7）
>> syms x;
>> y=exp(-x)*sin(2*x);
>> int(y,0,inf)
ans =
2/5
8）
>> syms x
f=sqrt(1+x);
taylor(f,0,9)
ans =
- (429*x^8)/32768 + (33*x^7)/2048 - (21*x^6)/1024 + (7*x^5)/256 - (5*x^4)/128 + x^3/16 - x^2/8 + x/2 + 1
9）
>> syms x y;
>> y=exp(sin(1/x));
>> dy=subs(diff(y,3),x,2)
dy =
-0.5826
10）求变上限函数对变量x的导数。
>> syms a t;
>> diff(int(sqrt(a+t),t,x,x^2))
Warning: Explicit integral could not be found.
ans =
2*x*(x^2 + a)^(1/2) - (a + x)^(1/2)
4、求点(1,1,4)到直线L: 的距离
>> M0=[1,1,4];M1=[3,0,1];M0M1=M1-M0;
v=[-1,0,2];
d=norm(cross(M0M1,v))/norm(v)
d =
1.0954
5、已知分别在下列条件下画出的图形：（要求贴图）
，在同一坐标系里作图
>> syms x;
>> fplot('(1/sqrt(2*pi))*exp(-((x)^2)/2)',[-3,3],'r')
>> hold on
>> fplot('(1/sqrt(2*pi))*exp(-((x-1)^2)/2)',[-3,3],'y')
>> hold on
>> fplot('(1/sqrt(2*pi))*exp(-((x+1)^2)/2)',[-3,3],'g')
>> hold off

，在同一坐标系里作图。
>> syms x;
fplot('(1/sqrt(2*pi))*exp(-((x)^2)/2)',[-3,3],'r')
hold on
fplot('(1/(sqrt(2*pi)*2))*exp(-((x)^2)/(2*2^2))',[-3,3],'y')
hold on
fplot('(1/(sqrt(2*pi)*4))*exp(-((x)^2)/(2*4^2))',[-3,3],'g')
hold off

6、画下列函数的图形：（要求贴图）
（1）
>> ezmesh('u*sin(t)','u*cos(t)','t/4',[0,20,0,2])

(2)
>> x=0:0.1:3;y=x;
[X Y]=meshgrid(x,y);
Z=sin(X*Y);
>> mesh(X,Y,Z)

（3）
ezmesh('sin(t)*(3+cos(u))','cos(t)*(3+cos(u))','sin(u)',[0,2*pi,0,2*pi])

7、 已知，在MATLAB命令窗口中建立A、B矩阵并对其进行以下操作：
(1) 计算矩阵A的行列式的值
>> A=[4,-2,2;-3,0,5;1,5,3];
>> det(A)
ans =
-158
(2) 分别计算下列各式：
>> A=[4,-2,2;-3,0,5;1,5,3];B=[1,3,4;-2,0,-3;2,-1,1];
>> 2*A-B
ans =
7 -7 0
-4 0 13
0 11 5
>> A*B
ans =
12 10 24
7 -14 -7
-3 0 -8
>> A.*B
ans =
4 -6 8
6 0 -15
2 -5 3
>> A*inv(B)
ans =
-0.0000 -0.0000 2.0000
-2.7143 -8.0000 -8.1429
2.4286 3.0000 2.2857
>> inv(A)*B
ans =
0.4873 0.4114 1.0000
0.3671 -0.4304 0.0000
-0.1076 0.2468 0.0000
>> A*A
ans =
24 2 4
-7 31 9
-8 13 36
>> A'
ans =
4 -3 1
-2 0 5
2 5 3
>>
8、 在MATLAB中分别利用矩阵的初等变换及函数rank、函数inv求下列矩阵的秩：
(1) 求 rank(A)=?
>> A=[1,-6,3,2;3,-5,4,0;-1,-11,2,4];
>> rank(A)
ans =
3
(2) 求。
>> B=[3,5,0,1;1,2,0,0;1,0,2,0;1,2,0,2]
>> inv(B)
ans =
2.0000 -4.0000 -0.0000 -1.0000
-1.0000 2.5000 0.0000 0.5000
-1.0000 2.0000 0.5000 0.5000
0 -0.5000 0 0.5000
9、在MATLAB中判断下列向量组是否线性相关，并找出向量组
中的一个最大线性无关组。
>> a1=[1 1 3 2]'
a2=[-1 1 -1 3]'
a3=[5 -2 8 9]'
a4=[-1 3 1 7]'
A= [a1, a2 ,a3 ,a4] ;[R jb]=rref(A)
a1 =
1
1
3
2
a2 =
-1
1
-1
3
a3 =
5
-2
8
9
a4 =
-1
3
1
7
R =
1.0000 0 0 1.0909
0 1.0000 0 1.7879
0 0 1.0000 -0.0606
0 0 0 0
jb =
1 2 3
>> A(:,jb)
ans =
1 -1 5
1 1 -2
3 -1 8
2 3 9
10、在MATLAB中判断下列方程组解的情况，若有多个解，写出通解。
（1）
一：
>> A=[1,-1,4,2;1,-1,-1,2;3,1,7,-2;1,-3,-12,6];
>> rank(A)
ans =
3
>> rref(A)
ans =
1 0 0 0
0 1 0 -2
0 0 1 0
0 0 0 0
二：
>> A=[1,-1,4,2;1,-1,-1,2;3,1,7,-2;1,-3,-12,6];
>> format rat
n=4;
RA=rank(A)
RA =
3
>> if(RA==n)
fprintf('%方程只有零解')
else
b=null(A,'r')
end
b =
0
2
0
1
>> syms k
X=k*b
X =
0
2*k
0
k
(2)
>> A=[2 3 1;1 -2 4;3 8 -2;4 -1 9];
b=[4 -5 13 -6]';
B=[A b];
>> n=3;
>> RA=rank(A)
RA =
2
>> RB=rank(B)
RB =
2
rref(B)
ans =
1 0 2 -1
0 1 -1 2
0 0 0 0
0 0 0 0
>> format rat
if RA==RB&RA==n %判断有唯一解
X=A\b
elseif RA==RB&RA<n %判断有无穷解
X=A\b %求特解
C=null(A,'r') %求AX=0的基础解系
else X='equition no solve' %判断无解
end
Warning: Rank deficient, rank = 2, tol = 8.9702e-015.
X =
0
3/2
-1/2
C =
-2
1
1
11、求矩阵 的逆矩阵 及特征值和特征向量。
A=[-2 1 1;0 2 0;-4 1 3];
>> a1=inv(A)
a1 =
-3/2 1/2 1/2
0 1/2 0
-2 1/2 1
>> [P,R]=eig(A)
P =
-985/1393 -528/2177 379/1257
0 0 379/419
-985/1393 -2112/2177 379/1257
R =
-1 0 0
0 2 0
0 0 2
A的三个特征值是：
r1=-1，r2=2，r3=2。
三个特征值分别对应的特征向量是
P1=[1 0 1];p2=[1 0 4];p3=[1 3 1]
12、化方阵为对角阵。
>> A=[2 2 -2;2 5 -4;-2 -4 5];
[P,D]=eig(A)
P =
-0.2981 0.8944 0.3333
-0.5963 -0.4472 0.6667
-0.7454 0 -0.6667
D =
1.0000 0 0
0 1.0000 0
0 0 10.0000
>> B=inv(P)*A*P
B =
1.0000 -0.0000 0.0000
0.0000 1.0000 0.0000
-0.0000 0 10.0000
程序说明：
所求得的特征值矩阵D即为矩阵A对角化后的对角矩阵，D和A相似。
13、求一个正交变换，将二次型化为标准型。
>> A=[5 -1 3;-1 5 -3;3 -3 3];
>> syms y1 y2 y3
y=[y1;y2;y3];
[P,D]=eig(A)
P =
881/2158 985/1393 -780/1351
-881/2158 985/1393 780/1351
-881/1079 0 -780/1351
D =
* 0 0
0 4 0
0 0 9
>> x=P*y
x =
(6^(1/2)*y1)/6 + (2^(1/2)*y2)/2 - (3^(1/2)*y3)/3
(2^(1/2)*y2)/2 - (6^(1/2)*y1)/6 + (3^(1/2)*y3)/3
- (3^(1/2)*y3)/3 - (2^(1/2)*3^(1/2)*y1)/3
>> f=[y1 y2 y3]*D*y
f =
- y1^2/2251799813685248 + 4*y2^2 + 9*y3^2
14、设，数列是否收敛？若收敛，其值为多少？精确到6位有效数字。
f=inline('(x+7/x)/2');
>> x0=3;
>> for i=1:20
x0=f(x0);
fprintf('%g,%g\n',i,x0);
end
1,2.66667
2,2.64583
3,2.64575
4,2.64575
5,2.64575
6,2.64575
7,2.64575
8,2.64575
9,2.64575
10,2.64575
11,2.64575
12,2.64575
13,2.64575
14,2.64575
15,2.64575
16,2.64575
17,2.64575
18,2.64575
19,2.64575
20,2.64575
该数列收敛于三，它的值是
15、设 是否收敛？若收敛，其值为多少？精确到17位有效数字。
（注：学号为单号的取，学号为双号的取）
>> f=inline('1/(x^8)');
x0=0;
for i=1:20
x0=(x0+f(i));
fprintf('%g , %.16f\n',i,x0);
end
1 , 1.0000000000000000
2 , 1.0039062500000000
3 , 1.0040586657902759
4 , 1.0040739245793384
5 , 1.0040764845793384
6 , 1.0040770799535192
7 , 1.0040772534200448
8 , 1.0040773130246896
9 , 1.0040773362552626
10 , 1.0040773462552626
11 , 1.0040773509203365
12 , 1.0040773532460168
13 , 1.0040773544719115
14 , 1.0040773551495150
15 , 1.0040773555396993
16 , 1.0040773557725300
17 , 1.0040773559158835
18 , 1.0040773560066281
19 , 1.0040773560655085
20 , 1.0040773561045711
>>
16、求二重极限
>> clear
>> syms x y;
>> f=(log(x+exp(y))/sqrt(x^2+y^2));
>> fx=limit(f,'x',1);
>> fxy=limit(fx,'y',0)
fxy =
log(2)
17、已知。
>> clear
syms x y z;
>> F=exp(x)-x*y*z;
>> Fx= diff(F, 'x')
Fx =
exp(x) - y*z
>> Fz= diff(F, 'z')
Fz =
-x*y
>> G=-Fx/Fz
G =
(exp(x) - y*z)/(x*y)
18、已知函数，求梯度。
一：
>> clear
syms x y z;
>> f=x^2+2*y^2+3*z^2+x*y+3*x-3*y-6*z;
>> dxyz=jacobian(f)
dxyz =
[ 2*x + y + 3, x + 4*y - 3, 6*z - 6]
二：
>> clear
>> syms x y z;
>> f=x^2+2*y^2+3*z^2+x*y+3*x-3*y-6*z;
>> gr=jacobian(f)
gr =
[ 2*x + y + 3, x + 4*y - 3, 6*z - 6]
19、计算积分，其中由直线围成。
>> A=int(int ((2-x-y),'y',x^2,x),'x',0,1)/2
A =
11/120
20、计算曲线积分，其中曲线。
clear
syms x y z t
x=cos(t);
y=sin(t);
z=t;
dx=diff(x,t);
dy=diff(y,t);
dz=diff(z,t);
ds=sqrt(dx^2+dy^2+dz^2);
f=z^2/(x^2+y^2);
I=int(f*ds,t,0,2*pi)
I =
(8*2^(1/2)*pi^3)/3
21、计算曲面积分，其中。
>> clear
>> syms x y z a;
>> z=sqrt(a^2-x^2-y^2);
>> f=x+y+z;
>> I=int(int(f,'y',0,sqrt(a^2-x^2)),'x',0,a)
I=
1/2*a^3+1/4*a^3*pi+1/3*a^2*(a^2)^(1/2)+1/3*(-1/2-1/4*pi)*a^3
22、求解二阶微分方程：。
>> clear
>> syms x y;
>> d_equa='D2y-10*Dy+9*y=exp(2*x)'
d_equa =
D2y-10*Dy+9*y=exp(2*x)
>> Condit= 'y(0)=6/7,Dy(0)=33/7'
Condit =
y(0)=6/7,Dy(0)=33/7
>> y1=dsolve( d_equa , Condit , 'x')
y1 =
exp(9*x)/2 - exp(2*x)/7 + exp(x)/2
23、求数项级数的和。
>> clear
>> syms n;
>> f=1/(n*(n+1));
>> I=symsum(f,n,1,inf)
I =
1
24、将函数展开为的幂级数。
>> clear
>> syms x;
>> f=1/x;
>> taylor(f,10,x,3)
ans =
(x - 3)^2/27 - x/9 - (x - 3)^3/81 + (x - 3)^4/243 - (x - 3)^5/729 + (x - 3)^6/2187 - (x - 3)^7/6561 + (x - 3)^8/19683 - (x - 3)^9/59049 + 2/3
25、能否找到一个分式线性函数，使它产生的迭代序列收敛到给定的数？用这种办法近似计算。
>> f=inline('(2+x^2)/(2*x)');
x1=2;
for i=1:20
x1=f(x1);
fprintf('%g,%g\n',i,x1);
end;
1,1.5
2,1.41667
3,1.41422
4,1.41421
5,1.41421
6,1.41421
7,1.41421
8,1.41421
9,1.41421
10,1.41421
11,1.41421
12,1.41421
13,1.41421
14,1.41421
15,1.41421
16,1.41421
17,1.41421
18,1.41421
19,1.41421
20,1.41421
26、函数的迭代是否会产生混沌？
>> x1=0:0.05:0.5;
y1=2*x1;
x2=0.5:0.05:1;
y2=2*(1-x2);
figure
plot(x1,y1,x2,y2)
gtext('2*x')
gtext('2*(1-x)')

27、函数称为Logistic映射，试从“蜘蛛网”图观察它取初值为产生的迭代序列的收敛性，将观察记录填人下表，作出图形。若出现循环，请指出它的周期。（要求贴图）

f=inline('3.3*x*(1-x)');
x=linspace(1,202,202);y=linspace(1,202,202);
x(1)=0.5;
y(1)=0;x(2)=x(1);y(2)=x(1);
for i=1:100
x(1+2*i)=x(2*i);
x(2+2*i)=f(x(1+2*i));
y(1+2*i)=x(2+2*i);
y(2+2*i)=y(1+2*i);
end
plot(x,y,'r');
hold on;
syms x y;
y=x;
ezplot(x,[0,1]);
ezplot(f(x),[0,1]);
axis([0,1,0,3.3/4]);
hold off T=0.35
hold on
f=inline('3.5*x*(1-x)');
x=linspace(1,202,202);y=linspace(1,202,202);
x(1)=0.5;
y(1)=0;x(2)=x(1);y(2)=x(1);
for i=1:100
x(1+2*i)=x(2*i);
x(2+2*i)=f(x(1+2*i));
y(1+2*i)=x(2+2*i);
y(2+2*i)=y(1+2*i);
end
plot(x,y,'r');
hold on;
syms x y;
y=x;
ezplot(x,[0,1]);
ezplot(f(x),[0,1]);
axis([0,1,0,3.5/4]);
hold off T=0.4
hold on
f=inline('3.56*x*(1-x)');
x=linspace(1,202,202);y=linspace(1,202,202);
x(1)=0.5;
y(1)=0;x(2)=x(1);y(2)=x(1);
for i=1:100
x(1+2*i)=x(2*i);
x(2+2*i)=f(x(1+2*i));
y(1+2*i)=x(2+2*i);
y(2+2*i)=y(1+2*i);
end
plot(x,y,'r');
hold on;
syms x y;
y=x;
ezplot(x,[0,1]);
ezplot(f(x),[0,1]);
axis([0,1,0,3.56/4]);
hold off
hold on
f=inline('3.568*x*(1-x)');
x=linspace(1,202,202);y=linspace(1,202,202);
x(1)=0.5;
y(1)=0;x(2)=x(1);y(2)=x(1);
for i=1:100
x(1+2*i)=x(2*i);
x(2+2*i)=f(x(1+2*i));
y(1+2*i)=x(2+2*i);
y(2+2*i)=y(1+2*i);
end
plot(x,y,'r');
hold on;
syms x y;
y=x;
ezplot(x,[0,1]);
ezplot(f(x),[0,1]);
axis([0,1,0,3.568/4]);
hold on
f=inline('3.6*x*(1-x)');
x=linspace(1,202,202);y=linspace(1,202,202);
x(1)=0.5;
y(1)=0;x(2)=x(1);y(2)=x(1);
for i=1:100
x(1+2*i)=x(2*i);
x(2+2*i)=f(x(1+2*i));
y(1+2*i)=x(2+2*i);
y(2+2*i)=y(1+2*i);
end
plot(x,y,'r');
hold on;
syms x y;
y=x;
ezplot(x,[0,1]);
ezplot(f(x),[0,1]);
axis([0,1,0,3.6/4]);
hold off
hold on
f=inline('3.84*x*(1-x)');
x=linspace(1,202,202);y=linspace(1,202,202);
x(1)=0.5;
y(1)=0;x(2)=x(1);y(2)=x(1);
for i=1:100
x(1+2*i)=x(2*i);
x(2+2*i)=f(x(1+2*i));
y(1+2*i)=x(2+2*i);
y(2+2*i)=y(1+2*i);
end
plot(x,y,'r');
hold on;
syms x y;
y=x;
ezplot(x,[0,1]);
ezplot(f(x),[0,1]);
axis([0,1,0,3.84/4]);
hold off
表 Logistic迭代的收敛性
a
3.3
3.5
3.56
3.568
3.6
3.84
序列收敛情况
不收敛
不收敛
不收敛
不收敛
不收敛
不收敛
28、由函数与构成的二维迭代Martin迭代。现观察其当时取初值为所得到的二维迭代散点图有什么变化。（要求贴图）
function Martin (a,b,c N)
f=@(x,y)(y-sign(x)*sqrt(abs(a*x-c)));
g=@(x)(a-x);
m=[0;0];
for n=1:N
m(:,n+1)=[f(m(1,n)),m(2,n),g(m(1,n))];
end
plot(m(1,:),m(2,:),'kx');
axis equal
Martin(4.52555120,2,-300,500)
书上62页
29、对，，求出平面映射的通项，并画出这些点的散点图。
A=[4,2;1,3];
t=[];
for i=1:20
x=2*rand(2,1)-1;
t(length(t)+1,1:2)=x;
for j=1:40
x=A*x;
t(length(t)+1,1:2)=x;
end
end
plot(t(:,1),t(:,2),'*')
grid('on')

30、对及随机给出的，观察数列.该数列有极限吗?
31、若该地区的天气分为三种状态：晴、阴、雨。对应的转移矩阵为：

且，试根据这些数据来求出若干天之后的天气状态，并找出其特点（取4位有效数字）。
>> A1=[3/4,1/2,1/4;1/8,1/4,1/2;1/8,1/4,1/4];
p=[0.5;0.25;0.25];
for i=1:20
p(:,i+1)=A1*p(:,i);
end
p
p =
Columns 1 through 7
0.5000 0.5625 0.5938 0.6035 0.6069 0.6081 0.6085
0.2500 0.2500 0.2266 0.2207 0.2185 0.2178 0.2175
0.2500 0.1875 0.1797 0.1758 0.1746 0.1741 0.1740
Columns 8 through 14
0.6086 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087
0.2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174
0.1739 0.1739 0.1739 0.1739 0.1739 0.1739 0.1739
Columns 15 through 21
0.6087 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087
0.2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174
0.1739 0.1739 0.1739 0.1739 0.1739 0.1739 0.1739
32、对于上例中的，求出矩阵的特征值与特征向量，并将特征向量与上例中的结论作对比。
>> A=[3/4 1/2 1/4;1/8 1/4 1/2;1/8 1/4 1/4];
>> [P,R]=eig(A)
P =
-0.9094 -0.8069 0.3437
-0.3248 0.5116 -0.8133
-0.2598 0.2953 0.4695
R =
1.0000 0 0
0 0.3415 0
0 0 -0.0915
特征值是r1=1,r2=0,3415,r3=-0.0915;
特征向量是R1=[],R2=[],R3=[]
对应于特征值1的特征向量P1=[-0.9094,-0.3248,-0.2598]
因为： P=[0.6087, 0.2174, 0.1739]=-1.494P1
结论：
属于同一特征值的特征向量可以相差k倍
33、编程找出的所有勾股数，并问：能否利用通项表示?
>> for b=1:995
a=sqrt((b+5)^2-b^2);
if(a==floor(a))
fprintf('a=%i,b=%i,c=%i\n',a,b,b+5)
end
end
a=15,b=20,c=25
a=25,b=60,c=65
a=35,b=120,c=125
a=45,b=200,c=205
a=55,b=300,c=305
a=65,b=420,c=425
a=75,b=560,c=565
a=85,b=720,c=725
a=95,b=900,c=905
34、用Monte Carlo方法计算圆周率。
clear
a=10000;
for i=1:a
x=rand();
f=inline('1/(1+(x^2))');
F(1,i)=f(x);
end
jifen=mean(F)
jifen