南邮2013MATLAB数学实验答案(全)

第一次练习

教学要求：熟练掌握Matlab 软件的基本命令和操作，会作二维、三维几何图形，能够用Matlab 软件解决微积分、线性代数与解析几何中的计算问题。 补充命令

vpa(x,n) 显示x 的n 位有效数字，教材102页

fplot(‘f(x)’,[a,b]) 函数作图命令，画出f(x)在区间[a,b]上的图形 在下面的题目中m 为你的学号的后3位（1-9班）或4位（10班以上） 1.1 计算30sin lim

x mx mx x →-与3

sin lim

x mx mx

x →∞- syms x

limit((902*x-sin(902*x))/x^3) ans =

366935404/3

limit((902*x-sin(902*x))/x^3,inf) ans = 0 1.2 cos

1000

x

mx

y e =，求''y syms x

diff(exp(x)*cos(902*x/1000),2) ans =

(46599*cos((451*x)/500)*exp(x))/250000 - (451*sin((451*x)/500)*exp(x))/250 1.3 计算

22

11

00

x y e

dxdy +??

dblquad(@(x,y) exp(x.^2+y.^2),0,1,0,1) ans = 2.1394

1.4 计算4

2

2

4x dx m x +? syms x

int(x^4/(902^2+4*x^2)) ans =

(91733851*atan(x/451))/4 - (203401*x)/4 + x^3/12 1.5 (10)cos ,x y e mx y =求

syms x

diff(exp(x)*cos(902*x),10) ans =

-356485076957717053044344387763*cos(902*x)*exp(x)-3952323024277642494822005884*sin(902*x)*exp(x)

1.6 0x =的泰勒展式(最高次幂为4).

syms x

taylor(sqrt(902/1000+x),5,x) ans =

-(9765625*451^(1/2)*500^(1/2)*x^4)/82743933602 +(15625*451^(1/2)*500^(1/2)*x^3)/91733851

-(125*451^(1/2)*500^(1/2)*x^2)/406802 + (451^(1/2)*500^(1/2)*x)/902 +(451^(1/2)*500^(1/2))/500 1.7 Fibonacci 数列{}n x 的定义是121,1x x ==12,(3,4,)n n n x x x n --=+= 用循环语句编

程给出该数列的前20项（要求将结果用向量的形式给出）。 x=[1,1]; for n=3:20

x(n)=x(n-1)+x(n-2); end x x=

Columns 1 through 10

1 1

2

3 5 8 13 21 3

4 5

5 Columns 11 through 20

89 144 233 377 610 987 1597 2584 4181 6765

1.8 对矩阵211020411000A m ?? ?- ?= ?

?- ?

?

?，求该矩阵的逆矩阵，特征值，特征向量，行列式，

计算6

A ，并求矩阵,P D （D 是对角矩阵），使得1

A PDP -=。

A=[-2,1,1;0,2,0;-4,1,902/1000];inv(A)

ans =

0.4107 0.0223 -0.4554 0 0.5000 0 1.8215 -0.4554 -0.9107 eig(A) ans =

-0.5490 + 1.3764i -0.5490 - 1.3764i 2.0000 det(A) ans =

4.3920

[P,D]=eig(A)

P = %特征向量

0.3245 - 0.3078i 0.3245 + 0.3078i 0.2425 0 0 0.9701 0.8944 0.8944 0.0000 D =

-0.5490 + 1.3764i 0 0 0 -0.5490 - 1.3764i 0 0 0 2.0000 P*D^6*inv(P) %A^6的值 ans =

15.3661 12.1585 + 0.0000i -5.8531 0 64.0000 0 23.4124 -5.8531 + 0.0000i -1.6196

1.9 作出如下函数的图形（注：先用M 文件定义函数，再用fplot 进行函数作图）：

1202

()12(1)1

2

x x f x x x ?

≤≤??=??-<≤??

m 文件:

function y=fenduan(x)

if x<=1/2 y=2*x

else x<=1 y=2-2*x

end end

执行函数：fplot('fenduan',[0,1]); grid on

title('第1.9题图') 得下图：

0.0.0.0.0.0.0.0.0.第1.9题图

1.10 在同一坐标系下作出下面两条空间曲线（要求两条曲线用不同的颜色表示）

（1）cos sin x t y t z t =??=??=? （2）2cos 2sin x t y t z t =??

=??=?

t=-10:0.01:10; x1=cos(t); y1=sin(t); z1=t;

plot3(x1,y1,z1); hold on

x2=cos(2*t); y2=sin(2*t); z2=t;

plot3(x2,y2,z2,'m'); grid on

title('第1.10题图') 得下图：

第1.10题图

1.11 已知422134305,203153211A B m -???? ? ?

=-=-- ? ? ? ?-????

，在MATLAB 命令窗口中建立A 、B 矩阵并

对其进行以下操作：

(1) 计算矩阵A 的行列式的值det()A

(2) 分别计算下列各式：1122,*,.*,,,,T A B A B A B AB A B A A --- 解：A=[4,-2,2;-3,0,5;1,5*902,3]; B=[1,3,4;-2,0,3;2,-1,1];

det(A) ans =

-117288 2*A-B ans =

7 -7 0

-4 0 7 0 9021 5 A*B ans =

12 10 12

7 -14 -7 -9013 0 13537 A.*B ans =

4 -6 8

6 0 15 2 -4510 3 A*inv(B) ans =

1.0e+003 *

-0.0000 0 0.0020 0.0000 0.0016 0.0001 1.0311 -0.9016 -1.4167 inv(A)*B ans =

0.3463 0.5767 0.5383 0.0005 -0.0006 -0.0005 -0.1922 0.3460 0.9230 A*A ans =

24 9012 4

-7 22556 9 -13523 13528 22561 A' ans =

4 -3 1

-2 0 4510 2 5 3

1.12

已知22

()2()x f x μσ--=

分别在下列条件下画出)(x f 的图形：

(1)/600m σ=,μ分别为0,1,1-(在同一坐标系上作图); (2)0μ=,σ分别为1,2,4,/100m (在同一坐标系上作图). (1)x=-5:0.1:5;

h=inline('1/sqrt(2*pi)/s*exp(-(x-mu).^2/(2*s^2))'); y1=h(0,902/600,x);y2=h(-1,902/600,x);y3=h(1,902/600,x); plot(x,y1,'b',x,y2,'m',x,y3,'y')

grid on

title('第1.12题')

第1.12题

(2) z1=h(0,1,x);z2=h(0,2,x);z3=h(0,4,x); z4=h(0,902/100,x);

plot(x,z1,x,z2,'y',x,z3,'m',x,z4, 'g') grid on

title('第1.12题')

z1=h(0,1,x);z2=h(0,2,x);z3=h(0,4,x); z4=h(0,902/100,x);

第1.12题

1.13 作出2

4

z mx y =+的函数图形。

x=-10:0.2:10;y=x;

[X Y]=meshgrid(x,y);Z=902*X.^2+Y.^4; mesh(X,Y,Z); title('第1.13题')

第1.13题

1.14对于方程5

0.10200

m

x x -

-=，

先画出左边的函数在合适的区间上的图形，借助于软件中的方程求根的命令求出所有的实根,找出函数的单调区间,结合高等数学的知识说明函数为什么在这些区间上是单调的,以及该方程确实只有你求出的这些实根。最后写出你做此题的体会。 解：作图程序：（注：x 范围的选择是经过试探而得到的）

x=-1.7:0.02:1.7;y=x.^5-902/200*x-0.1;

plot(x,y);grid on; title('第1.14题')

第1.14题

由图形观察，在x=-1.5,x=0,x=1.5附近各有一个实根 solve('x^5-902/200*x-0.1') ans =

-1.4516870267499636199995749888894 -0.022172950190557703188753959027919 1.4627751059480654637229232196174

1.4573364935933870280941533926624*i + 0.0055424354962279297327028641499658 0.0055424354962279297327028641499658 - 1.4573364935933870280941533926624*i 三个实根的近似值分别为：-1.4517,-0.0222,1.4628

由图形可以看出，函数在区间(,1)-∞-单调上升，在区间(1,1)-单调下降，在区间(1,)∞单调上升。 syms x

diff('x^5-902/200*x-0.1',x)

结果为5*x^4-4.51

solve('5*x^4-902/200') ans =

-(451^(1/4)*500^(3/4))/500 (451^(1/4)*500^(3/4))/500 -(451^(1/4)*500^(3/4)*i)/500 (451^(1/4)*500^(3/4)*i)/500 vpa(ans) ans =

-0.97454440927373918149075795211629 0.97454440927373918149075795211629 -0.97454440927373918149075795211629*i 0.97454440927373918149075795211629*i 得到两个实根：-0.9745与0.9745

可以验证导函数在)9745.0,(--∞内为正，函数单调上升 导函数在)9745.0,9745.0(-内为负，函数单调下降 导函数在),9745.0(∞内为正，函数单调上升 根据函数的单调性，最多有3个实根。

1.15 求23m 0x e x -=的所有根。（先画图后求解）（要求贴图）

作图命令：（注：x 范围的选择是经过试探而得到的） x=-5:0.001:15;y=exp(x)-3*902*x.^2; plot(x,y);grid on; title('第1.15题图') 得到下图

进一步细化

6

第1.15题图

x=-0.05:0.0001:0.05;y=exp(x)-3*902*x.^2; plot(x,y);grid on; title('第1.15题图')

x=10:0.001:15;y=exp(x)-3*902*x.^2; plot(x,y);grid on; title('第1.15题图')

第1.15题

图

6

第1.15题图

可看出解在-0.02，0.02，13附近，进一步求得 fzero('exp(x)-3*902*x^2',0.02) ans =

0.0194

fzero('exp(x)-3*902*x^2',-0.02) ans =

-0.0190

fzero('exp(x)-3*902*x^2',13) ans =

13.0391

第二次练习

教学要求：要求学生掌握迭代、混沌的判断方法，以及利用迭代思想解决实际问题。

2.1 设11

()/23n n

n

m x x x x +?

=+???=?，数列{}n x 是否收敛？若收敛，其值为多少？精确到8位有效数字。

解：程序代码如下（m=902）：

f=inline('(x+902/x)/2'); x0=3; for i=1:20; x0=f(x0);

fprintf('%g %8f\n',i,x0);

end

1 151.833333

2 78.887029

3 45.160551

4 32.566867

5 30.131864

6 30.033476

7 30.033315

8 30.033315

9 30.033315

……

19 30.033315

20 30.033315

由运行结果可以看出，，数列{}n x收敛，其值为30.03315。

2.2求出分式线性函数

2

12

1

(),()

x x m

f x f x

x m x m

-+

==

++

的不动点，再编程判断它们的迭代序

列是否收敛。

解：取m=1000.

（1）程序如下：

f=inline('(x-1)/(x+1000)');

x0=2;

for i=1:20;

x0=f(x0);

fprintf('%g,%g\n',i,x0);

end

运行结果：

1,0.000998004 11,-0.001001

2,-0.000999001 12,-0.001001

3,-0.001001 13,-0.001001

4,-0.001001 14,-0.001001

5,-0.001001 15,-0.001001

6,-0.001001 16,-0.001001

7,-0.001001 17,-0.001001

8,-0.001001 18,-0.001001 9,-0.001001 19,-0.001001 10,-0.001001 20,-0.001001

由运行结果可以看出，，分式线性函数收敛，其值为-0.001001。易见函数的不动点为-0.001001（吸引点）。 （2）程序如下：

f=inline('(x+1000000)/(x+1000)'); x0=2; for i=1:20; x0=f(x0);

fprintf('%g,%g\n',i,x0); end 运行结果：

1,998.006 11,618.332 2,500.999 12,618.302 3,666.557 13,618.314 4,600.439 14,618.309 5,625.204 15,618.311 6,615.692 16,618.31 7,619.311 17,618.311 8,617.929 18,618.31 9,618.456 19,618.31 10,618.255 20,618.31

由运行结果可以看出，，分式线性函数收敛，其值为618.31。易见函数的不动点为618.31（吸引点）。

2.3 下面函数的迭代是否会产生混沌？（56页练习7（1））

1202

()12(1)12

x x f x x x ?

≤≤??=??-<≤??

解：程序如下：

f=inline('1-2*abs(x-1/2)'); x=[];

y=[]; x(1)=rand();

y(1)=0;x(2)=x(1);y(2)=f(x(1)); for i=1:100; x(1+2*i)=y(2*i); x(2+2*i)=x(1+2*i); y(2+2*i)=f(x(2+2*i)); end

plot(x,y,'r'); hold on; syms x;

ezplot(x,[0,1/2]); ezplot(f(x),[0,1]); axis([0,1/2,0,1]); hold off 运行结果：

0.05

0.1

0.15

0.2

0.250.3

0.35

0.4

0.45

0.5

x

1 -

2 abs(x - 1/2)

2.4 函数()(1)(01)f x x x x α=-≤≤称为Logistic 映射，试从“蜘蛛网”图观察它取初值

为00.5x =产生的迭代序列的收敛性，将观察记录填人下表，若出现循环，请指出它的周期．（56页练习8）

解：当α=3.3时，程序代码如下：f=inline('3.3*x*(1-x)');

x=[];

y=[];

x(1)=0.5;

y(1)=0;x(2)=x(1);y(2)=f(x(1)); for i=1:1000;

x(1+2*i)=y(2*i);

x(2+2*i)=x(1+2*i);

y(1+2*i)=x(1+2*i);

y(2+2*i)=f(x(2+2*i));

end

plot (x,y,'r');

hold on;

syms x;

ezplot(x,[0,1]);

ezplot(f(x),[0,1]);

axis([0,1,0,1]);

hold off运行结果：

-(33 x (x - 1))/10

x

当α=3.5时，上述程序稍加修改，得：

x

当α=3.56时，得：

-(89 x (x - 1))/25

x

当α=3.568时，得：

-(446 x (x - 1))/125

x

当α=3.6时，得：

x

当 =3.84时，得：

x

a b c为其它的值会得到什么图形？参考下表（取自63页练

2.5对于Martin迭代，取参数,,

习13

）

解：取m=10000；迭代次数N=20000；

在M-文件里面输入代码：

function Martin(a,b,c,N)

f=@(x,y)(y-sign(x)*sqrt(abs(b*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（10000，10000，10000，20000），得：

x 104

执行Martin（-10000，-10000，10000，20000），得：

x 104

执行Martin（-10000，10，-10000，20000），得：

执行Martin（10，10，0.5，20000），得：

执行Martin（10，10000，-10000，20000），得：

执行Martin （100，1000，-10，20000），得：

执行Martin （-1000，17，4，20000），得：

2.6 能否找到分式函数

2

ax b

cx dx e

+++（其中,,,,a b c d e 是整数）,使它产生的迭代序列（迭代

。如果迭代收敛,那么迭代的初值与收敛的速度有什么关系.写出你做此题的体会. 提示：教材54页练习4的一些分析。

若分式线性函数()ax b

f x cx d

+=+为()f x 的不动点，因此

=

化简得：(2)(0c b d a -+-=。

若,,,a b c d 为整数，易见2,b c d a ==。

取满足这种条件的不同的,,,a b c d 以及迭代初值进行编。 解：

2ax b cx dx e +++迭代收敛到指定的数3902，则3

902为2

ax b cx dx e

+++的不动点，所以

3

902=

e

d c b a +++3

2

3

3902)902(902，解得a=e,d=0,b=902c

取m=902；根据上述提示，取：a=e=1,b=902,c=1,d=0,程序如下：

f=inline('(x+902)/(x^2+1)');

x0=1;

for i=1:100; x0=f(x0);

fprintf('%g %g\n',i,x0); end

结果如下

1 451.5

2 0.00663958

3 901.967

4 0.002217413

5 901.998

6 0.0022173

7 901.99

8 8 0.0022173

9 901.998 10 0.0022173 11 901.998 12 0.0022173 13 901.998 14 0.0022173 15 901.998 16 0.0022173 17 901.998 ……

95 901.998 96 0.0022173 97 901.998

98 0.0022173

99 901.998

100 0.002217

初值为-1，结果为

1 450.5

2 0.00666416

3 901.967

4 0.00221742

5 901.998

6 0.0022173

7 901.998

8 0.0022173

9 901.998

10 0.0022173

11 901.998

12 0.0022173

13 901.998

14 0.0022173

15 901.998

16 0.0022173

17 901.998

18 0.0022173

19 901.998

……

95 901.998

96 0.0022173

97 901.998

98 0.0022173

99 901.998

100 0.0022173

初值为1000，结果为

1 0.001902

2 901.999

3 0.0022173

4 901.998

5 0.0022173

6 901.998

7 0.0022173

8 901.998

9 0.0022173

10 901.998

11 0.0022173

12 901.998

13 0.0022173

14 901.998 15 0.0022173 16 901.998 17 0.0022173 18 901.998 19 0.0022173 20 901.998 ……

93 0.0022173 94 901.998 95 0.0022173 96 901.998 97 0.0022173 98 901.998 99 0.0022173 100 901.998

第三次练习

教学要求：理解线性映射的思想，会用线性映射和特征值的思想方法解决诸如天气等实际问题。

3.1 对???

? ??=3124A ，T

T x x )2,1(),()0(2)0(1=，求出}{n x 的通项.

>> syms n

>> A=sym('[4,2;1,3]');x=[1;2];[P,D]=eig(A) %没有sym 下面的矩阵就会显示为小数 P = [ -1, 2] [ 1, 1] D = [ 2, 0] [ 0, 5]

>> An=P*D^n*inv(P) An =

[ 2^n/3 + (2*5^n)/3, (2*5^n)/3 - (2*2^n)/3] [ 5^n/3 - 2^n/3, (2*2^n)/3 + 5^n/3] >> xn=An*x xn =

2*5^n - 2^n 2^n + 5^n 3.2 =

B ???

? ??=3.01.02.04.0101A 对于练习1中的B ，T T x x )2,1(),()0(2)0(1=，求出}{n x 的通项. >> syms n

>> A=sym('[2/5,1/5;1/10,3/10]'); x=[1;2];[P,D]=eig(A) P =