考试说明

① 在线考试（50%）
② 实验报告（50%）

Part1 实验报告

利用高数知识，结果：m*(m+1)/2，√2

syms x
simplify( diff(exp(x) * cos(m*x/1000), 2) )
q = simplify( diff(exp(x) * cos(m*x/1000), 6) );
subs(q, 0)

syms x
q = int( exp((-m)*x^2), x, 0, inf );
vpa(q, 6)

syms x; 
f = (m/1000+x) ^(1/3);
simplify(taylor(f, 5, x))

x = [];
x = randint(1, 2, [1, m])
for n = 3:10
    x(n) = x(n-1) + x(n-2);
end
x

A = [ 4, -2, 2; -3, 0, 5; 1, 5, 3*m ];
B = [ 1, 3, 4; -2, 0, -3; 2, -1, 1 ];
det(A)
inv(A)
[P,D] = eig(A)
A*inv(B)
inv(A)*B

1）M文件
f.m :
function y = f(x)
if 0<=x & x<=1/2
    y = 2.0*x;
elseif 1/2<x & x<=1
    y = 2.0*(1-x);
end
2）画图命令
fplot(@f, [0, 1])

t = -m/10: 0.01: m/10; 
x1 = cos(t)*m/20; y1 = sin(t)*m/20; z1=t;
plot3(x1, y1, z1, 'b');
grid on

figure
t = -m/10: 0.01: m/10; 
x2 = cos(t) + t.*sin(t);
y2 = sin(t) - t.*cos(t);
z2 = (-1)*t;
plot3(x2, y2, z2, 'r');
grid on

function y = f(x,a)
if (x>0)
    y = a*exp(-a*x);
else
    y=0;
end
fplot( 'f(x,1000/m)', [-1,10], 'r' ) 
fplot( 'f(x,500/m)', [-1,10], 'c' ) 
fplot( 'f(x,400/m)', [-1,10], 'k' ) 
fplot( 'f(x,100/m)', [-1,10], 'm' ) 

syms x y 
ezplot( 'sin(x^2 + m/1000*y^2)- cos(x*y)', [-6,6,-6,6]) 

x = -8: 0.1: 8; y = -8: 0.1: 8; 
[X Y] = meshgrid(x, y);
Z = m*X.^2 + Y.^4;
mesh(X, Y, Z);

1）画图
syms x
fplot('exp(x) - 3*m*x^2/(m+100)', [-4,4])
grid on

2）近似求根
fsolve( 'exp(x) - 3*m*x^2/(m+100)', -1)
fsolve( 'exp(x) - 3*m*x^2/(m+100)', 1)
fsolve( 'exp(x) - 3*m*x^2/(m+100)', 3)

3）单调区间
syms x
diff( exp(x) - 3*m*x^2/(m+100) )   
>> ans= exp(x) - (954*x)/209

syms x
fsolve(' exp(x) - (954*x)/209', 0)   
fsolve(' exp(x) - (954*x)/209', 1.5)   
先求出导数，再求出导数为零的点，结合图形可得单调区间 (高数知识)

f = inline( '(x+m/x)/2' );

x0 = 3;
for i = 1:20:
x0 = f(x0);
fprintf( '%g, %2.8f \n', i, x0 );
end
x0 = -3;
for i = 1:20:
x0 = f(x0);
fprintf( '%g, %2.8f \n', i, x0 );
end

q=solve('(x-1)/(x+m)-x')
vpa(q,4)
f=inline('(x-1)/(x+m)');
x0=2;
for i=1:20;
x0=f(x0);
fprintf('%g,%g\n',i,x0);
end
由运行结果可以看出，分式线性函数为吸引点。

q=solve('(x+m^2)/(x+m)-x')
vpa(q,4)
f=inline('(x+m^2)/(x+m)');
x0=2;
for i=1:20;
x0=f(x0);
fprintf('%g,%g\n',i,x0);
end
由运行结果可以看出，分式线性函数吸引点。

function  y=f(x)
xx=eval(x);
if   xx>=0  &  xx<=1/2
y=2*x;
elseif  xx>1/2 & xx<=1
       y=2*(1-x);
end
function p=dd(x0,n)
p=sym(x0);
for i=2:n
    p(i)=f(p(i-1));
end
>> dd(3/10,20)
ans =  
[ 3/10, 3/5, 4/5, 2/5, 4/5, 2/5, 4/5, 2/5, 4/5, 2/5, 4/5, 2/5, 4/5, 2/5, 4/5, 2/5, 4/5, 2/5, 4/5, 2/5]
结论：T=2，不产生混沌

当 α=2.8 时：
f=inline('2.8*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
y
plot (x,y,'r');
hold on;
syms x;
ezplot(x,[0,1]);
ezplot(f(x),[0,1]);
axis([0,1,0,1]);
hold off
运行结果：T=2

当 α=3.4 时，T=2
当 α=3.6 时，产生混沌
当 α=3.84 时，T=3

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

m=10000；N=20000；
>> Martin（m，m，m，N）
>> Martin（-m，-m，m，N）
>> Martin（-m，m/1000，-m，N）
>> Martin（m/1000，m/1000，0.5，N）
>> Martin（m/1000，m，-m，N）
>> Martin（m/100，m/10，-10，N）
>> Martin（-m/10，17，4，N）

syms x
>> f=inline('(100*x+503)/(x^2+100)');
>> x0=10;
>> for i=1:10;
x0=f(x0);
fprintf('%g,%g\n',i,x0);
end

(1) m=412; n=200;
>> s=[];s(1)=(sin(1))^m;
>> for k=2:n
s(k)=s(k-1)+1/(k*k)*(sin(k))^m;
plot(k,s(k))
hold on
end
s
s=[……-0.0052   -0.0052   -0.0052]
该级数收敛，其和为-0.0052.

(2) m=412; n=200; a=m/1000;
>> s=[];s(2)=0;
>> for k=3:n
s(k)=s(k-1)+1/(log(k))^a;
plot(k,s(k))
hold on
end
s
s=[0  0  0.9813  1.9177  2.8265  3.7159  ...  146.8743  147.5897  148.3049]
该级数发散

A=sym('[4,2;1,3]'); [P,D]=eig(A)
Q=inv(P) syms n;
xn=P*(D.^n)*Q*[1;2]

A=sym('[2/5,1/5;1/10,3/10]');   
[P,D]=eig(A)
Q=inv(P) 
xn=P*(D.^n)*Q*[1;2] 

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

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;
for i=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,y 一起输出)
x   
y m
结论：{xn},{yn} 及 m(xn) 的极限都存在

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)
结论： A 的绝对值最大特征值与上面的 m(xn) 的极限相等

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
结论： 9 天后，天气状态趋于稳定，P* =（ 0.6087,0.2174,0.1739 )T

>> [P,D]=eig(A2)
分析：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* 一致。

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

c-b=4 时通项：{a, b, c} {4(u+1), 2(u^2+u), 2(u^2+2u+2)}
...
c-b=7 时通项：{a, b, c} {5(2u+1), 5(2u^2+2u), 5(2u^2+2u+1)}
综上：
c-b=奇数 时通项：{a, b, c} {k(2u+1), k(2u^2+2u), k(2u^2+2u+1)}
c-b=偶数 时通项：{a, b, c} {k(u+1), k(u^2+u), k(u^2+2u+2)}

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

(1) gcd(a,b)
(2) randint(1,1,[u,v])  
m=m; s=0;
for i=1:100*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)

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

Part2 命令总结

• 求两个正整数的最大公约数

gcd(a,b)

• 求a的近似整数

fix(a)
floor(a)
ceil(a)
round(a)
fix 朝零方向取整
floor 朝负无穷方向取整
ceil 朝正无穷方向取整
round 四舍五入取整

• 清理窗口命令

clc

• 产生100×100全是1的矩阵

ones(100)

• 绘制二维曲线

ezplot(f)
f是关于x的函数  （-2*pi < x < 2*pi）

y=[y1 y2 y3 y4 ...];
plot(y)

fplot(@(x) (f))
f是关于x的函数  （-2*pi < x < 2*pi）

• 绘制三维曲线

x=[1,2,3]；
y=[-1,0,1,2]；
z=[1,2,3
2,3,4
-1,2,3
0,2,9]；
surf(x,y,z)
z为4行3列的矩阵

• 求符号积分

计算定积分
syms x
f = f(x);
int (f, x, a, b)

计算不定积分
syms x
f = f(x);
int (f, x)

• 求数值积分

(1) 梯形法
format long
ac = @(x)f;
x = a : pi/100 : b;
y = ac(x);
s = trapz(x, y)

(2) 变步长辛普森法
ac = @(x)f;
s = quad(ac, a, b)

(3) 二重积分
f = @(x,y)f(x);
I = dblquad(f, a, b, c, d)

• 产生0到1的随机数

rand

• 获取矩阵中的元素

取出矩阵A 第2行所有元素
A(2,:)
取出矩阵A 第1列所有元素
A(:,1)
取出矩阵A 第2行第2列元素
A(2,2)

• 求导数，导数值

求 f=sin(x) 的导数
syms x;
f = sin(x);
diff(f)

求 f=sin(x) 的导数在x=2的导数值
syms x;
f = sin(x);
f1 = diff(f);
subs(f1, x, 2)

• 求矩阵特征值，特征向量

[P,D] = eig(矩阵)

Part3 知识点

1）迭代的可视化也称为蜘蛛网图
2）迭代次数足够大之后，产生序列：0.4，0.8，0.4，0.8...
则迭代序列周期性重复，且周期为2
3）logistic混沌映射：https://blog.csdn.net/qiluofei/article/details/1837562
4）产生线性映射迭代序列的要素是：初始向量
5）本原勾股数：https://blog.csdn.net/weixin_42282037/article/details/97373925
6）选择语句有：if ... else，switch ... case ... end
7）符号 ; 的作用是区分行，取消运行显示
8）循环语句有：while, for
9）任意两个正整数互质的概率为：6/π^2 ≈ 0.61
10）蒙特·卡罗方法：http://blog.sina.com.cn/s/blog_13dd6d82a0102vcgz.html
12）若勾股数 (a,b,c) 中，c-b=1，则它一定是本原勾股数
若勾股数(a,b,c)中，c-b = 3或15，则它一定不是本原勾股数
13）不动点迭代法：https://blog.csdn.net/qq_39521554/article/details/79835632
14）混沌：序列没有规律，杂乱无章
15）画图时如果需要网格线，命令：grid on

Part4 慕课练习

1）单元测试一

m文件：

function y=f(x)

if 0<=x&x<=1/2

    y=2*x;

elseif 1/2<x&x<=1

    y=2*(1-x)

end

end

命令运行窗口：

>> fplot('f(x)',[0,1])

单元测试一1.png

>>fplot('exp(x)-3*733/833*x^2',[-1,5])

单元测试一2_1.png

>> fplot('exp(x)-3*733/833*x^2',-1)
ans =

  -0.4833

>> fplot('exp(x)-3*733/833*x^2',3)

ans =

  3.4440

>> fplot('exp(x)-3*733/833*x^2',2)
ans =

  1.0302

>> fplot('exp(x)-(4398*x)/833',[-1,4])

单元测试一2_2.png

>>fsolve('exp(x)-(4398*x)/833',2)

>>ans =

  2.6314

>>fsolve('exp(x)-(4398*x)/833',0)

>>ans =

  0.2410

>>

单调增区间【负无穷，0.2410】，【2.6314，正无穷】

单调减区间【0.2410，2.6314】

收敛   27.074

>>f=inline('(x+733/x)/2')

f =
     内联函数:
     f(x) = (x+733/x)/2

>> x1=3;
>> for i=1:100
x1=f(x1);
fprintf('%g %g\n',i,x1)
;end
1 123.667
2 64.7969
3 38.0546
4 28.6582
5 27.1178
6 27.074
7 27.074
8 27.074
9 27.074
10 27.074
11 27.074
12 27.074
13 27.074
14 27.074
15 27.074
16 27.074
17 27.074
18 27.074
19 27.074
20 27.074
21 27.074
22 27.074
23 27.074
24 27.074
25 27.074
26 27.074
27 27.074
28 27.074
29 27.074
30 27.074
31 27.074
32 27.074
33 27.074
34 27.074
35 27.074
36 27.074
37 27.074
38 27.074
39 27.074
40 27.074
41 27.074
42 27.074
43 27.074
44 27.074
45 27.074
46 27.074
47 27.074
48 27.074
49 27.074
50 27.074
51 27.074
52 27.074
53 27.074
54 27.074
55 27.074
56 27.074
57 27.074
58 27.074
59 27.074
60 27.074
61 27.074
62 27.074
63 27.074
64 27.074
65 27.074
66 27.074
67 27.074
68 27.074
69 27.074
70 27.074
71 27.074
72 27.074
73 27.074
74 27.074
75 27.074
76 27.074
77 27.074
78 27.074
79 27.074
80 27.074
81 27.074
82 27.074
83 27.074
84 27.074
85 27.074
86 27.074
87 27.074
88 27.074
89 27.074
90 27.074
91 27.074
92 27.074
93 27.074
94 27.074
95 27.074
96 27.074
97 27.074
98 27.074
99 27.074
100 27.074

m文件：

function g(f,x0,n)
close all
for i=1:n
    axis([50,n,0,1]);
    if i>50
        plot(i,f(x0),'.');
        hold on;pause(0.05);
    end
    x0=f(x0);
end
hold off
clear all;

运行命令窗口：

g(@(x)2.8*x*(1-x),0.5,200)

2.8收敛

g(@(x)3.4*x*(1-x),0.5,200)

3.4循环周期2

g(@(x)3.84*x*(1-x),0.5,200)

3.84循环周期3

g(@(x)3.6*x*(1-x),0.5,200)

3.6混沌

2）单元测试二

>> A=[4,2;1,3];
>> [P,D]=eig(A)

P =



    0.8944   -0.7071
    0.4472    0.7071



D =



     5     0
     0     2



>> E=sym('[5^n,0;0,2^n]')


E =


[ 5^n,0]

[ 0,2^n]


>>F=P*E*P^-1;

>>X0=[1;2];

>>K=F*X0


K =


2*5^n-2^n

2^n+5^n

>> A=[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)=A*p(:,i);
end
>> p

p =

  1 至 9 列

    0.5000    0.5625    0.5938    0.6035    0.6069    0.6081    0.6085    0.6086    0.6087
    0.2500    0.2500    0.2266    0.2207    0.2185    0.2178    0.2175    0.2174    0.2174
    0.2500    0.1875    0.1797    0.1758    0.1746    0.1741    0.1740    0.1739    0.1739

  10 至 18 列

    0.6087    0.6087    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.2174    0.2174
    0.1739    0.1739    0.1739    0.1739    0.1739    0.1739    0.1739    0.1739    0.1739

  19 至 21 列

    0.6087    0.6087    0.6087
    0.2174    0.2174    0.2174
    0.1739    0.1739    0.1739

>>
特点：各种天气的概率逐渐趋于一个定值

>> I=0;
>> m=[];
>> n=1000;
>> for a=1:n
for c=a+1:n
b=sqrt(c^2-a^2);
if(b==floor(b))&(b>a)&(c==b+2)
I=I+1;m(:,I)=[a,b,c];
end
end
end
>> m

m =

  1 至 17 列

     6     8    10    12    14    16    18    20    22    24    26    28    30    32    34    36    38
     8    15    24    35    48    63    80    99   120   143   168   195   224   255   288   323   360
    10    17    26    37    50    65    82   101   122   145   170   197   226   257   290   325   362

  18 至 29 列

    40    42    44    46    48    50    52    54    56    58    60    62
   399   440   483   528   575   624   675   728   783   840   899   960
   401   442   485   530   577   626   677   730   785   842   901   962

>>

公式：2m b=m^2-1 c=m^+1(m>2,m为整数）

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

>>sqrt(6*m/s)


ans =


    3.03.2

3）期末考试

1）[d,p,q] = gcd(56,126) 的输出结果：d=14, p=-2, q=1
2）a=12时，在200内的本原勾股数有 2 组
3）逗号的作用是：区分列
4）迭代序列：要么收敛，要么发散 or 混沌 or 周期性重复
5）“分形” 一词是由美国IBM公司数学家 Benoit B.Mandelbrot 提出来的
分形：组成部分以某种方式与整体相似的形体，具有比例性、置换不变性
6）a==floor(a) 可以用来判断 a 是否为整数
7）一元函数的迭代过程，可以用 蜘蛛网图 来直观的显示
8）Feigenbaum 图可以分析 f(x) = αx(1-x) 的迭代行为
9）如果迭代函数存在不动点，则一定收敛

10）吸引点就是一个所有附近点，在迭代过程中都趋向于它的不动点
11）近似计算三次根号2的方法是迭代法

① 若今天晴，则明天晴的概率为 5/9；若今天阴雨，则明天晴的概率为 4/13
则该地区天气模型的转移矩阵 A=

② 若 P(0) 表示天晴和阴雨情况，A是转移矩阵
则 m 天之后，该地区天气情况可以表示为：P(m) = P(0) * A^m
③ 已知初始向量和迭代矩阵，求迭代序列的的命令是：[P,D] = eig(A)
④ 对线性迭代过程进行归一化，是为了使迭代序列分量满足绝对值最大不超过1
⑤ 如果归一化后产生序列的极限存在，则最大分量绝对值所产生序列的极限也存在
⑥ 归一化可增加迭代序列的收敛性
⑦ 在 “天气问题” 中，天气状态趋于稳定，那么转移矩阵存在一个特征值
⑧ 已知：f(x) = e^x * cos(1001 * x/1000)
求 f(x) 的两阶导数的命令是：diff (exp(x) * cos(1001*x/1000), 2)
⑨ 求方阵A的行列式的命令是：det(A)