第一次练习
教学要求：熟练掌握Matlab软件的基本命令和操作，会作二维、三维几何图形，能够用Matlab软件解决微积分、线性代数与解析几何中的计算问题。
补充命令
vpa(x,n) 显示x的n位有效数字，教材102页
fplot(‘f(x)’,[a,b]) 函数作图命令，画出f(x)在区间[a,b]上的图形
在下面的题目中为你的学号的后3位（1-9班）或4位（10班以上）
1.1 计算与
程序：
syms x
limit((627*x-sin(627*x))/x^3,x,0)
结果：
1003003001/6
程序：
syms x
limit((627*x-sin(627*x))/x^3,x,inf)
结果：
0
1.2 ，求
程序：
syms x
diff(exp(x)*cos(627*x/1000),2)
结果：
-2001/1000000*exp(x)*cos(1001/1000*x)-1001/500*exp(x)*sin(1001/1000*x)
1.3 计算
程序：
dblquad(@(x,y) exp(x.^2+y.^2),0,1,0,1)
结果：
2.13935019514228
1.4 计算
程序：
syms x
int(x^4/(627^2+4*x^2))
结果：
1/12*x^3-1002001/16*x+1003003001/32*atan(2/627*x)
1.5
程序：
syms x
diff(exp(x)*cos(627*x),10)
结果：
- 9389137388146839380380277888*cos(627*x)*exp(x) -149759579095532896918284384*sin(627*x)*exp(x)
1.6 给出在的泰勒展式(最高次幂为4).
程序：
syms x
taylor(sqrt(627/1000+x),4)
结果：
(62500*627^(1/2)*1000^(1/2)*x^3)/246491883- (125*627^(1/2)*1000^(1/2)*x^2)/393129 + (627^(1/2)*1000^(1/2)*x)/1254 +(627^(1/2)*1000^(1/2))/1000
1.7 Fibonacci数列的定义是，
用循环语句编程给出该数列的前20项（要求将结果用向量的形式给出）。
程序：
x=[1,1];
for n=3:20
x(n)=x(n-1)+x(n-2);
end
x
结果：
Columns 1 through 10
1 1 2 3 5 8 13 21 34 55
Columns 11 through 20
89 144 233 377 610 987 1597 2584 4181 6765
1.8 对矩阵，求该矩阵的逆矩阵，特征值，特征向量，行列式，计算，并求矩阵（是对角矩阵），使得。
程序与结果：
a=[-2,1,1;0,2,0;-4,1,627 /1000];
inv(a)
0.2283 0.0679 -0.3642
0 0.5000 0
1.4567 -0.3642 -0.7283
eig(a)
-0.6865 + 1.5082i
-0.6865 - 1.5082i
2.0000
[p,d]=eig(a)
p =
0.2937 - 0.3372i 0.2937 +0.3372i 0.2425
0 0 0.9701
0.8944 0.8944 0.0000
注：p的列向量为特征向量
d =
-0.6865 + 1.5082i 0 0
0 -0.6865 -1.5082i 0
0 0 2.0000
a^6
11.9680 13.0080 -4.9910
0 64.0000 0
19.9640 -4.9910 -3.0100
1.9 作出如下函数的图形（注：先用M文件定义函数，再用fplot进行函数作图）：
函数文件f.m:
function y=f(x)
if 0<=x&x<=1/2
y=2.0*x;
else 1/2<x&x<=1
y=2.0*(1-x);
end
程序：fplot(@f,[0,1])
1.10 在同一坐标系下作出下面两条空间曲线（要求两条曲线用不同的颜色表示）
（1） （2）
程序：
t=-10:0.01:10;
x1=cos(t);
y1=sin(t);
z1=t;
plot3(x1,y1,z1,'k');hold on
x2=cos(;
y2=sin(2*t);
z2=t;
plot3(x2,y2,z2,'r');hold off
1.11 已知，在MATLAB命令窗口中建立A、B矩阵并对其进行以下操作：
(1) 计算矩阵A的行列式的值
(2) 分别计算下列各式：
解：（1）程序：
a=[4,-2,2;-3,0,5;1,5*627,3];b=[1,3,4;-2,0,3;2,-1,1];det(a)
-81538
(2) 2*a-b
7 -7 0
-4 0 7
0 6271 5
a*b 12 10 12
7 -14 -7
-6263 0 9412
a.*b 4 -6 8
6 0 15
2 -3135 3
a*inv(b) -0.0000 0 2.0000
0.0286 1.6000 0.0857
716.8286 -626.6000 -984.5143
inv(a)*b 0.3464 0.5766 0.5382
0.0007 -0.0008 -0.0007
-0.1921 0.3460 0.9229
a^2 24 6262 4
-7 15681 9
-9398 9403 15686
A' 4 -3 1
-2 0 5005
2 5 3
1.12 已知分别在下列条件下画出的图形：
(1),分别为(在同一坐标系上作图);
(2),分别为(在同一坐标系上作图).
(1)程序：
x=-5:0.1:5;
h=inline('1/sqrt(2*pi)/s*exp(-(x-mu).^2/(2*s^2))');
y1=h(0,1001/600,x);y2=h(-1,1001/600,x);y3=h(1,1001/600,x);
plot(x,y1,'r+',x,y2,'k-',x,y3,'b*')
(2)程序：
z1=h(0,1,x);z2=h(0,2,x);z3=h(0,4,x); z4=h(0,1001/100,x);
plot(x,z1,'r+',x,z2,'k-',x,z3,'b*',x,z4, 'y:')
1.13 作出的函数图形。
程序：x=-5:0.1:5;y=-10:0.1:10;
[X Y]=meshgrid(x,y);Z=627*X.^2+Y.^4;
mesh(X,Y,Z);
1.14对于方程，先画出左边的函数在合适的区间上的图形，借助于软件中的方程求根的命令求出所有的实根,找出函数的单调区间,结合高等数学的知识说明函数为什么在这些区间上是单调的,以及该方程确实只有你求出的这些实根。最后写出你做此题的体会。
解：作图程序：（注：x范围的选择是经过试探而得到的）
x=-1.7:0.02:1.7;y=x.^5-627/200*x-0.1;
plot(x,y)；grid on;
由图形观察，在x=-1.5,x=0,x=1.5附近各有一个实根
求根程序：solve('x^5-627/200*x-0.1')
结果：
-1.4906852047544424910680160298802
-0.019980020616193485540810824654811
1.500676329192316320110463906588700421518815060273901630060819255
1.495764171739511484743570420265584278874768154469167692755643546*i+ 0.004994448089159828249181473973153383352756761740138087409772356778
0.004994448089159828249181473973153383352756761740138087409772356778 -1.4957641717395114847435704202656*i
三个实根的近似值分别为：
-1.490685，-0.019980，1.500676
由图形可以看出，函数在区间单调上升，在区间单调下降，在区间单调上升。
diff('x^5-1001/200*x-0.1',x)
结果为5*x^4-1001/200
solve('5*x^4-1001/200.')得到两个实根：-1.0002499与1.0002499
可以验证导函数在内为正，函数单调上升
导函数在内为负，函数单调下降
导函数在内为正，函数单调上升
根据函数的单调性，最多有3个实根。
1.15 求的所有根。（先画图后求解）（要求贴图）
作图命令：（注：x范围的选择是经过试探而得到的）
x=-5:0.001:15;y=exp(x)-3*627*x.^2;
plot(x,y);grid on;
可以看出，在(-5,5)内可能有根，在(10,15)内有1个根
将(-5,5)内图形加细，最终发现在(-0.03,0.03)内有两个根。
用solve('exp(x)-3*627.0*x^2',x)可以求出3个根为：
.18417113274368129311145677478702e-1
13.162041092091149185726742857195
-.18084038990284796648194134222365e-1
即：-0.018417,0.018084,13.16204
第二次练习
教学要求：要求学生掌握迭代、混沌的判断方法，以及利用迭代思想解决实际问题。
2.1 设，数列是否收敛？若收敛，其值为多少？精确到8位有效数字。
解：程序代码如下（m=627）：
>> f=inline('(x+627/x)/2');
x0=3;
for i=1:20;
x0=f(x0);
fprintf('%g,%g\n',i,x0);
end
运行结果：
1,3
2,3
3,3
4,3
5,3
6,3
7,3
8,3
9,3
10,3
11,3
12,3
13,3
14,3
15,3
16,3
17,3
18,3
19,3
20,3
由运行结果可以看出，，数列收敛，其值为3
2.2 求出分式线性函数的不动点，再编程判断它们的迭代序列是否收敛。
解：取m=627.
（1）程序如下：
f=inline('(x-1)/(x+627)');
x0=2;
for i=1:20;
x0=f(x0);
fprintf('%g,%g\n',i,x0);
end
运行结果：
1,0.00158983
2,-0.00159236
3,-0.00159744
4,-0.00159745
5,-0.00159745
6,-0.00159745
7,-0.00159745
8,-0.00159745
9,-0.00159745
10,-0.00159745
11,-0.00159745
12,-0.00159745
13,-0.00159745
14,-0.00159745
15,-0.00159745
16,-0.00159745
17,-0.00159745
18,-0.00159745
19,-0.00159745
20,-0.00159745
由运行结果可以看出，，分式线性函数收敛，其值为-0.00159745。易见函数的不动点为--0.00159745（吸引点）。
（2）程序如下：
f=inline('(x+393129)/(x+627)');
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））
解：程序如下：
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
运行结果：
2.4 函数称为Logistic映射，试从“蜘蛛网”图观察它取初值为产生的迭代序列的收敛性，将观察记录填人下表，若出现循环，请指出它的周期．（56页练习8）
| 3.3 | 3.5 | 3.56 | 3.568 | 3.6 | 3.84 | |
| 序列收敛情况 | T=2 | T=4 | T=8 | T=9 | 混沌 | 混沌 |
解：当=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运行结果：
当=3.5时，上述程序稍加修改，得：
当=3.56时，得：
当=3.568时，得：
当=3.6时，得：
当=3.84时，得：
2.5 对于Martin迭代，取参数为其它的值会得到什么图形？参考下表（取自63页练习13）
| m | m | m |
| -m | -m | m |
| -m | m/1000 | -m |
| m/1000 | m/1000 | 0.5 |
| m/1000 | m | -m |
| m/100 | m/10 | -10 |
| -m/10 | 17 | 4 |
解：取m=627；迭代次数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），得：
执行Martin（-10000，-10000，10000，20000），得：
执行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 能否找到分式函数（其中是整数）,使它产生的迭代序列（迭代的初始值也是整数）收敛到(对于为整数的学号，请改为求)。如果迭代收敛,那么迭代的初值与收敛的速度有什么关系.写出你做此题的体会.
提示：教材54页练习4的一些分析。
若分式线性函数的迭代收敛到指定的数，则为的不动点，因此
化简得：。
若为整数，易见。
取满足这种条件的不同的以及迭代初值进行编。
解：取m=10000；根据上述提示，取：
a=e=1,b=10000,c=1,d=0.
程序如下（初值为1200）：
f=inline('(x+9)/(x^2+1)');
x0=1;
for i=1:100;
x0=f(x0);
fprintf('%g,%g\n',i,x0);
end
运行结果如下：
1,5
2,0.538462
3,7.3945
4,0.294449
5,8.55291
6,0.236714
7,8.74661
8,0.228979
9,8.7692
10,0.228106
11,8.77169
12,0.22801
13,8.77197
14,0.227999
15,8.772
16,0.227998
17,8.772
18,0.227998
19,8.772
20,0.227998
21,8.772
22,0.227998
23,8.772
24,0.227998
25,8.772
26,0.227998
27,8.772
28,0.227998
29,8.772
30,0.227998
31,8.772
32,0.227998
33,8.772
34,0.227998
35,8.772
36,0.227998
37,8.772
38,0.227998
39,8.772
40,0.227998
41,8.772
42,0.227998
43,8.772
44,0.227998
45,8.772
46,0.227998
47,8.772
48,0.227998
49,8.772
50,0.227998
51,8.772
52,0.227998
53,8.772
54,0.227998
55,8.772
56,0.227998
57,8.772
58,0.227998
59,8.772
60,0.227998
61,8.772
62,0.227998
63,8.772
64,0.227998
65,8.772
66,0.227998
67,8.772
68,0.227998
69,8.772
70,0.227998
71,8.772
72,0.227998
73,8.772
74,0.227998
75,8.772
76,0.227998
77,8.772
78,0.227998
79,8.772
80,0.227998
81,8.772
82,0.227998
83,8.772
84,0.227998
85,8.772
86,0.227998
87,8.772
88,0.227998
89,8.772
90,0.227998
91,8.772
92,0.227998
93,8.772
94,0.227998
95,8.772
96,0.227998
97,8.772
98,0.227998
99,8.772
100,0.227998
若初值取为1000，运行结果：
1,0.011
2,9998.8
3,0.000200036
4,10000
5,0.0002
6,10000
7,0.0002
8,10000
9,0.0002
10,10000
11,0.0002
12,10000
13,0.0002
14,10000
15,0.0002
16,10000
17,0.0002
18,10000
19,0.0002
20,10000
21,0.0002
22,10000
23,0.0002
24,10000
25,0.0002
26,10000
27,0.0002
28,10000
29,0.0002
30,10000
31,0.0002
32,10000
33,0.0002
34,10000
35,0.0002
36,10000
37,0.0002
38,10000
39,0.0002
40,10000
41,0.0002
42,10000
43,0.0002
44,10000
45,0.0002
46,10000
47,0.0002
48,10000
49,0.0002
50,10000
51,0.0002
52,10000
53,0.0002
54,10000
55,0.0002
56,10000
57,0.0002
58,10000
59,0.0002
60,10000
61,0.0002
62,10000
63,0.0002
64,10000
65,0.0002
66,10000
67,0.0002
68,10000
69,0.0002
70,10000
71,0.0002
72,10000
73,0.0002
74,10000
75,0.0002
76,10000
77,0.0002
78,10000
79,0.0002
80,10000
81,0.0002
82,10000
83,0.0002
84,10000
85,0.0002
86,10000
87,0.0002
88,10000
89,0.0002
90,10000
91,0.0002
92,10000
93,0.0002
94,10000
95,0.0002
96,10000
97,0.0002
98,10000
99,0.0002
100,10000
若初值取为-1，运行结果：
1,4999.5
2,0.0006001
3,10000
4,0.0002
5,10000
6,0.0002
7,10000
8,0.0002
9,10000
10,0.0002
11,10000
12,0.0002
13,10000
14,0.0002
15,10000
16,0.0002
17,10000
18,0.0002
19,10000
20,0.0002
21,10000
22,0.0002
23,10000
24,0.0002
25,10000
26,0.0002
27,10000
28,0.0002
29,10000
30,0.0002
31,10000
32,0.0002
33,10000
34,0.0002
35,10000
36,0.0002
37,10000
38,0.0002
39,10000
40,0.0002
41,10000
42,0.0002
43,10000
44,0.0002
45,10000
46,0.0002
47,10000
48,0.0002
49,10000
50,0.0002
51,10000
52,0.0002
53,10000
54,0.0002
55,10000
56,0.0002
57,10000
58,0.0002
59,10000
60,0.0002
61,10000
62,0.0002
63,10000
64,0.0002
65,10000
66,0.0002
67,10000
68,0.0002
69,10000
70,0.0002
71,10000
72,0.0002
73,10000
74,0.0002
75,10000
76,0.0002
77,10000
78,0.0002
79,10000
80,0.0002
81,10000
82,0.0002
83,10000
84,0.0002
85,10000
86,0.0002
87,10000
88,0.0002
89,10000
90,0.0002
91,10000
92,0.0002
93,10000
94,0.0002
95,10000
96,0.0002
97,10000
98,0.0002
99,10000
100,0.0002
第三次练习
教学要求：理解线性映射的思想，会用线性映射和特征值的思想方法解决诸如天气等实际问题。
3.1 对，，求出的通项.
程序：
A=sym('[4,2;1,3]');
[P,D]=eig(A)
Q=inv(P)
syms n;
xn=P*(D.^n)*Q*[1;2]
结果：
P =
[ 2,-1]
[ 1, 1]
D =
[ 5, 0]
[ 0, 2]
Q =
[ 1/3, 1/3]
[-1/3, 2/3]
xn =
2*5^n-2^n
5^n+2^n
3.2 对于练习1中的，，求出的通项.
程序：
A=sym('[2/5,1/5;1/10,3/10]'); %没有sym下面的矩阵就会显示为小数
[P,D]=eig(A)
Q=inv(P)
xn=P*(D.^n)*Q*[1;2]
结果：
P =
[ 2, -1]
[ 1, 1]
D =
[ 1/2, 0]
[ 0, 1/5]
Q =
[ 1/3, 1/3]
[ -1/3, 2/3]
xn =
2*(1/2)^n-(1/5)^n
(1/2)^n+(1/5)^n
3.3 对随机给出的，观察数列.该数列有极限吗?
>> A=[4,2;1,3];
a=[];
x=2*rand(2,1)-1;
for i=1:20
a(i,1:2)=x;
x=A*x;
end
for i=1:20
if a(i,1)==0
else t=a(i,2)/a(i,1);
fprintf('%g,%g\n',i,t);
end
end
结论：在迭代18次后，发现数列存在极限为0.5
1,-0.597298
2,-0.282275
3,0.0445866
4,0.277259
5,0.402189
6,0.459283
7,0.483443
8,0.493333
9,0.497326
10,0.498929
11,0.499572
12,0.499829
13,0.499931
14,0.499973
15,0.499989
16,0.499996
17,0.499998
18,0.499999
19,0.5
20,0.5
3.4 对120页中的例子，继续计算.观察及的极限是否存在. （120页练习9）
>>A=[2.1,3.4,-1.2,2.3;0.8,-0.3,4.1,2.8;2.3,7.9,-1.5,1.4;3.5,7.2,1.7,-9.0];
x0=[1;2;3;4];
x=A*x0;
fori=1:1:100
a=max(x);
b=min(x);
m=a*(abs(a)>abs(b))+b*(abs(a)<=abs(b));
y=x/m;
x=A*y;
end
x %也可以用f0，不能把x1,y一起输出
y
m
程序输出：
x1 =
0.9819
3.2889
-1.2890
-11.2213
y =
-0.0875
-0.2931
0.1149
1.0000
m =
-11.2213
结论：及的极限都存在.
3.5 求出的所有特征值与特征向量，并与上一题的结论作对比. （121页练习10）
>> A=[2.1,3.4,-1.2,2.3;0.8,-0.3,4.1,2.8;2.3,7.9,-1.5,1.4;3.5,7.2,1.7,-9.0];
[P,D]=eig(A)
P =
-0.3779 -0.8848 -0.0832 -0.3908
-0.5367 0.3575 -0.2786 0.4777
-0.6473 0.2988 0.1092 -0.7442
-0.3874 -0.0015 0.9505 0.2555
D =
7.2300 0 0 0
0 1.1352 0 0
0 0 -11.2213 0
0 0 0 -5.8439
结论：A的绝对值最大特征值等于上面的的极限相等，为什么呢？
还有，P的第三列也就是-11.2213对应的特征向量和上题求解到的y也有系数关系，两者都是-11.2213的特征向量。
3.6 设，对问题2求出若干天之后的天气状态，并找出其特点（取4位有效数字）. （122页练习12）
>> A2=[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:1:20
P(:,i+1)=A2*P(:,i);
end
P
P =
Columns 1 through 14
0.5000 0.5625 0.5938 0.6035 0.6069 0.6081 0.6085 0.6086 0.6087 0.6087 0.6087 0.6087 0.6087 0.6087
0.2500 0.2500 0.2266 0.2207 0.2185 0.2178 0.2175 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174 0.2174
0.2500 0.1875 0.1797 0.1758 0.1746 0.1741 0.1740 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
结论：9天后，天气状态趋于稳定P*=（0.6087,0.2174,0.1739）T
3.7 对于问题2，求出矩阵的特征值与特征向量，并将特征向量与上一题中的结论作对比. （122页练习14）
>> [P,D]=eig(A2)
P =
-0.9094 -0.8069 0.3437
-0.3248 0.5116 -0.8133
-0.2598 0.2953 0.4695
D =
1.0000 0 0
0 0.3415 0
0 0 -0.0915
分析：事实上，q=k(-0.9094, -0.3248, -0.2598)T均为特征向量，而上题中P*的3个分量之和为1，可令k(-0.9094, -0.3248, -0.2598)T=1，得k=-0.6696.有q=(0.6087, 0.2174, 0.1739),与P*一致。
3.8 对问题1，设为的两个线性无关的特征向量，若
，具体求出上述的，将表示成的线性组合，求的具体表达式，并求时的极限，与已知结论作比较. （123页练习16）
>> A=[3/4,7/18;1/4,11/18];
[P,D]=eig(A);
syms k pk;
a=solve(‘u*P(1,1)+v*P(1,2)-1/2’,’u*P(2,1)+v*P(2,2)-1/2’,’u’,’v’);
pk=a.u*D(1,1).^k*P(:,1)+a.v*D(2,2).^k*P(:,2)
pk =
-5/46*(13/36)^k+14/23
5/46*(13/36)^k+9/23
或者：
p0=[1/2;1/2];
[P,D]=eig(sym(A));
B=inv(sym(P))*p0
B =
5/46
9/23
syms k
pk=B(1,1)*D(1,1).^k*P(:,1)+B(2,1)*D(2,2).^k*P(:,2)
pk =
-5/46*(13/36)^k+14/23
5/46*(13/36)^k+9/23
>>vpa(limit(pk,k,100),10)
ans
=
.6086956522
.3913043478
结论：和用练习12中用迭代的方法求得的结果是一样的。
第四次练习
教学要求：会利用软件求勾股数，并且能够分析勾股数之间的关系。会解简单的近似计算问题。
4.1 求满足，的所有勾股数，能否类似于（11.8），把它们用一个公式表示出来？
程序：for b=1:998
a=sqrt((b+2)^2-b^2);
if(a==floor(a))
fprintf('a=%i,b=%i,c=%i\n',a,b,b+2)
end
end
运行结果：
a=4,b=3,c=5
a=6,b=8,c=10
a=8,b=15,c=17
a=10,b=24,c=26
a=12,b=35,c=37
a=14,b=48,c=50
a=16,b=63,c=65
a=18,b=80,c=82
a=20,b=99,c=101
a=22,b=120,c=122
a=24,b=143,c=145
a=26,b=168,c=170
a=28,b=195,c=197
a=30,b=224,c=226
a=32,b=255,c=257
a=34,b=288,c=290
a=36,b=323,c=325
a=38,b=360,c=362
a=40,b=399,c=401
a=42,b=440,c=442
a=44,b=483,c=485
a=46,b=528,c=530
a=48,b=575,c=577
a=50,b=624,c=626
a=52,b=675,c=677
a=54,b=728,c=730
a=56,b=783,c=785
a=58,b=840,c=842
a=60,b=899,c=901
a=62,b=960,c=962
勾股数，的解是：
以下是推导过程：
由，有
显然，，从而是2的倍数.设，代入上式得到：
因为，从而.
4.2 将上一题中改为,,,,分别找出所有的勾股数.将它们与时的结果进行比较，然后用公式表达其结果。
（1）时通项：
a=8,b=6,c=10
a=12,b=16,c=20
a=16,b=30,c=34
a=20,b=48,c=52
a=24,b=70,c=74
a=28,b=96,c=100
a=32,b=126,c=130
a=36,b=160,c=164
a=40,b=198,c=202
a=44,b=240,c=244
a=48,b=286,c=290
a=52,b=336,c=340
a=56,b=390,c=394
a=60,b=448,c=452
a=64,b=510,c=514
a=68,b=576,c=580
a=72,b=646,c=650
a=76,b=720,c=724
a=80,b=798,c=802
a=84,b=880,c=884
a=88,b=966,c=970
（2）5时通项：
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
（3）6时通项
a=12,b=9,c=15
a=18,b=24,c=30
a=24,b=45,c=51
a=30,b=72,c=78
a=36,b=105,c=111
a=42,b=144,c=150
a=48,b=189,c=195
a=54,b=240,c=246
a=60,b=297,c=303
a=66,b=360,c=366
a=72,b=429,c=435
a=78,b=504,c=510
a=84,b=585,c=591
a=90,b=672,c=678
a=96,b=765,c=771
a=102,b=864,c=870
a=108,b=969,c=975
（4）7时通项
a=21,b=28,c=35
a=35,b=84,c=91
a=49,b=168,c=175
a=63,b=280,c=287
a=77,b=420,c=427
a=91,b=588,c=595
a=105,b=784,c=791
综上：当c-b=k为奇数时，通项
当c-b=k为偶数时，通项
4.3 对，（），对哪些存在本原勾股数？（140页练习12）
程序：for k=1:200
for b=1:999
a=sqrt((b+k)^2-b^2);
if((a==floor(a))&gcd(gcd(a,b),(b+k))==1)
fprintf('%i,',k);
break;
end
end
end
运行结果：1,2,8,9,18,25,32,49,50,72,81,98,121,128,162,169,200,
4.4 设方程(11.15)的解构成数列,观察数列,,
,,.你能得到哪些等式？试根据这些等式推导出关于的递推关系式. （142页练习20）
解：1000以内解构成的数列,, , , 如下：
n 1 2 3 4 5 6
2 7 26 97 362 1351
1 4 15 56 209 780
3 11 41 153 571 2131
4 15 56 209 780 2911
1 3 11 41 153 571
我们发现这些解的关系似乎是：
=
=
因为=，所以。
有以下结论：
（4.1）
可以看成一个线性映射,令
,
（4.1）可写成：
4.5 选取对随机的，根据的概率求出的近似值。（取自130页练习7）
提示：(1)最大公约数的命令:gcd(a,b)
(2)randint(1,1,[u,v])产生一个在[u,v]区间上的随机整数
程序：
m=10000;s=0;
for i=1:m
a=randint(1,2,[1,10^9]);
if gcd(a(1),a(2))==1;
s=s+1;
end
end
pi=sqrt(6*m/s)
运行结果：
pi =
3.1510
4.6 用求定积分的MonteCarlo法近似计算。（102页练习16）
提示：Monte Carlo法近似计算的一个例子。
对于第一象限的正方形,内画出四分之一个圆
向该正方形区域内随即投点,则点落在扇形区域内的概率为.
投次点,落在扇形内的次数为,则,因此.
程序如下
n=100000;nc=0;
for i=1:n
x=rand;y=rand;
if(x^2+y^2<=1)
nc=nc+1;
end
end
pi=4*nc/n
解：程序：
a=0;b=1;m=1000;
H=1;s=0;
for i=1:m
xi=rand();
yi=H*rand();
if yi<sqrt(1-xi^2);
s=s+1;
end
end
pi=4*H*(b-a)*s/m
运行结果：
pi =
3.1480
综合题
一、方程求根探究
设方程
1.用matlab命令求该方程的所有根；
2.用迭代法求它的所有根，设迭代函数为
1)验证取该迭代函数的正确性；
2)分别取初值为-1.1,-1,-0.9,….,0.9,1,1.1,观察迭代结果，是否得到了原方程的根；
3)总结出使得迭代序列收敛到每个根时，初值的范围，比如要使迭代序列收敛到0（方程的一个根）初值应该在什么集合中选取，找出每个根的这样的初值集合。寻找的方法，可以是理论分析方法或数值实验方法。
解答：
1． 用solve命令即可求出所有解；
2． 1）提示：验证原方程与同解，以及验证迭代函数在不动点附近的导数绝对值是否小于1
2）代码省略，结果：初值取-1.1，-1，-0.9,-0.8，0.7时收敛到-1，初值取-0.7,0.8,0.9，1,1.1时收敛到1，初值取-0.6，-0.5，。。。，0.5,0.6时收敛到0；
3）在中分别取初值，最后分别收敛到-1,1,0；在内有无穷多个收敛到-1的初值小开区间，也有无穷多个收敛到0的小开区间，它们相互交替着；这种状态反射到内，即：在内有无穷多个收敛到1的初值小开区间，也有无穷多个收敛到0的小开区间，它们也是相互交替着，这些小区间与内小开区间对应。
二、1.三次曲线
(a)对k=0及其邻近的k的正值和负值，把的图形画在一个公共屏幕上。k的值是怎样影响到图形的形状的？
(b)求，它是一个二次函数。求该二次函数的判别式，对什么样的k值，该判别式为正？为零？为负？对什么k值有两个零点？一个或没有零点？现在请说明k的值对f 图形的形状有什么影响。
(c)对其他的k值做实验。当会发生什么情形？当呢？
解答：
(a)先用m文件定义函数f(x,k)=x^3+k*x
由fplot('[f(x,-0.6),f(x,-0.3),f(x,0),f(x,0.3),f(x,3)]',[-3,3])
得下图
可见k值不影响凹凸性，但单调性、单调区间以及极值随k值发生改变；k在0附近，小于0时，函数在某[-a,a]区间上单调递减，该区间长度随着k值增大而减小，k大于等于0时，函数单调增加。
(b) ；判别式，k为负、零、正时判别式分别为正、零、负；故k<0时，有两个零点，k=0时有一个零点，k>0时没有零点。以上说明原函数f(x)的驻点个数随着k值符号而变化，当k由负变正时，驻点由两个变成一个再到没有驻点，相应的单调区间由三个变成一个，单增单减单增，变为单增。
(c) k值越小单减区间长度越大，当时，f(x)单减区间变为无穷大对称区间，图形近乎垂直直线；当时，单增区间变为无穷大对称区间，图形近乎垂直直线。
2.四次曲线
(a)对k=-4及其邻近的k值，把的图形画在一个公共屏幕上。k的值是怎样影响到图形的形状的？
(b)求，它是一个二次函数。求该二次函数的判别式，对什么样的k值，该判别式为正？为零？为负？对什么k值有两个零点？一个或没有零点？现在请说明k的值对f 图形的形状有什么影响。
解答：
(a)先用m文件定义函数f(x,k)=x^4+k*x^3+6*x^2
fplot('[f(x,-4.2),f(x,-5),f(x,-4.5),f(x,-4),f(x,-3.5),f(x,-2.5)]',[-1,4])
得图
由图可以看出，在x<1时，图形受k值影响不大，x>1时k值对图形的影响比较显著，通过改变k值画图发现：在-4附近，k小于-4时，曲线在某[a,b](a>0)区间内是上凸的，在其他区间内上凹；k大于-4时，上面的凸区间不存在，也就是曲线总是上凹的。
(b) ，判别式，当时，判别式为0，时判别式大于0，时判别式小于0；也就是时有两个零点，时有一个零点时没有零点。由二阶导数与凹凸性的关系可知，在k=-4附近，（a）中关于曲线凹凸的判断基本上是正确的
三、对于级数，通过下面的步骤探索它的行为
1. 对于其部分和数列，当你试图求时，发生了什么？
解答：用命令sk=symsum(1/n^3/(sin(n))^2,1,k)及limit(sk,k,inf)得不到结果，命令symsum(1/n^3/(sin(n))^2,1,inf)也得不到结果。这表明极限可能并不存在。
2. 画出部分和数列的前100个点，它们是否显示出收敛？你估计极限是多少？
解答：前100个点图形如下
上图似乎显示着sk的极限存在，并且极限值约为4.8左右
3. 接着画出部分和数列的前200个点，用你自己的话论述部分和数列的行为。
此图可以更加确定，部分和数列sk的极限是存在的，结论跟2中的一样
4. 画出前400个点，当＝355时发生了什么？计算数355/113，通过你的计算解释当＝355时发生了什么。你猜测对的什么值同一现象可能还会出现，并通过实验加以验证。
解答：
此图否定了2与3的推断，因为部分和数列在＝355时发生了跳跃；355/113=3.141592920353983近似等于圆周率(约为3.141592653589793)，也就是355，而sin113=0，因此sin335的值很小，对应于部分和sk，在＝355时由于分母很小因而得到一个很大的加项，于是图形上的点发生了跳跃。我们可以通过观察或计算的倍数来获得sk的比较大的加项，由于710=355*2，因此sk在＝710时也会发生跳跃；我们也可直接由命令(1:500)*pi观察1500以内的数哪些接近的倍数（此略）。
另外，由的各种分数表示（近似）可知，以上的部分和sk在k=22时也会发生跳跃，因为。同上，当k=44,66,88,110,132等等时，sk也会发生跳跃，但由于误差扩大，跳跃幅度相对应该比较小。
四、通过本课程学习，谈谈你开设对这门课的认识，对教学以及上机实验提出自己的和建议
略
扫一扫
2082
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
