MA AB数学建模 乒乓球的弹跳和罗基斯帝模型

乒乓球的弹跳罗基斯第模型[问题]罗基斯第模型
一个乒乓球离球拍的高度为h0，落在球拍上反弹，设恢复系数为e，不计空气阻力。

(1)如果e为常数，讨论球的高度变化的规律。

如果e2与高度h n成线性关系
e2=μ(1–h n/H0)(2.1)其中H0是最大高度，μ是参数。

对于不同的参数讨论小球高度的变化规律。

(2)当参数连续变化时，分析最后分布的高度。

(3)计算前几个分岔点。

(4)用李雅普洛夫指数判断混沌的发生。

[解析](1)当球从高度h n下落到球拍上之前速度为
v(2.2)
n
球与球拍碰撞后反弹的速度为
v'n=ev n(2.3)球反弹的高度为
h n+1=e2h n(2.4)如果e<1，则球的反弹高度随次数不断减小；如果e=1，则球反弹后始终保持初始高度；如果e>1，例如球拍每次加一个向上的冲击力，则球的高度随次数不断增加。

e2与高度的线性关系说明：如果球的高度较大，则恢复系数较小，反之较大。

设相对高度为x n=h n/H0，则下一次上升的相对高度为
x n+1=μ(1–x n)x n，(n=0,1,2,…)(2.5)这是著名的罗基斯第模型。

由于相对高度0≤x n≤1，而(1–x n)x n的最大值为1/4，所以参数的值在0到4之间。

球的高度强烈依赖参数。

[算法](1)先取一个参数，再取一个相对高度，通过迭代算法计算下一次碰撞后的高度，画出高度点，依此类推。

再取另一高度参数，重新通过迭代算法计算高度，画出高度点，依此类推。

[程序]MATH2_1.m如下。

%乒乓球与球拍的碰撞高度
clear%清除变量
u=input('请输参数(参考值:0.5,2,3.25,3.5,3.56,3.8):');%键盘输入初始相对高度(1)
xn=0.9;%第1个的初始相对高度(2)
figure%开创图形窗口
plot(0,xn,'.')%画高度点
text(0,xn,num2str(xn),'FontSize',16)%标记第1个的初始高度
grid minor%加细网格
title(['乒乓球与球拍的碰撞高度(\it\mu\rm=',num2str(u),')'],'FontSize',16)%标题
n=50;%迭代次数
axis([0,n,0,1])%坐标范围
hold on%保持图像
for j=1:n%按次数循环
xn=u*(1-xn)*xn;%计算下一次的相对高度(3)
plot(j,xn,'.')%画高度点
end%结束循环
xn=0.1;%取初始相对高度(4)
plot(0,xn,'ro')%画高度点
text(0,xn,num2str(xn),'FontSize',16)%初始高度
for j=1:n%按次数循环
xn=u*(1-xn)*xn;%计算下一次的相对高度(5)
plot(j,xn,'ro')%画高度点
end%结束循环
[说明](1)程序执行时要用户用键盘输入参数，提供6个参数选择。

(2)取第1个较大的初始高度。

(3)迭代计算下一个高度。

(4)取第2个较小的初始高度。

在说明混沌时，将此句改写如下，使第2个高度比第1个高度大一点。

e=1e-8;%小量
xn=0.9+e;%取初始相对高度
text(0,0,num2str(e),'FontSize',16)%初始高度
(5)同样迭代计算下一个高度。

M2.1a图M2.1b图
[图示](1)如M2.1a图所示，当参数μ为0.5时，如果初始相对高度取0.9，球与球拍碰撞之后高度不断降低，最终的高度为零。

即使初始相对高度取0.2，球与球拍碰撞之后高度也不断降低，最终的高度为零。

M2.1c图M2.1d图
(2)如M2.1b图所示，当μ为2时，如果初始相对高度取0.9，球第一碰撞之后高度降低，以后碰撞则高度升高，最后碰撞保持一定的高度。

如果初始相对高度取0.2，则碰撞高度不
断增加，最后稳定在一定的高度。

高度稳定前的过程称为过渡过程或暂态过程，过渡过程与初始高度有关，但是最后高度的稳定与初始高度无关。

高度值称为不动点，即重复自身轨迹的点。

(3)如M2.1c图所示，当μ为3.25时，不论初始高度如何，经过过程期后，球最后在2个高度之中交替变化。

不动点的个数随参数的增加而增加。

(4)如M2.1d图所示，当μ为3.5时，不论初始高度如何，经过过程期后，球最后在4个高度之中交替变化，只是过渡过程稍微长一点。

不动点的个数随参数的增加而进一步增加。

(5)如M2.1e图所示，当μ为3.57时，不论初始高度如何，经过过渡期后，球最后在8个高度之中交替变化。

(6)如M2.1f图所示，当μ为3.8时，第1个初始相对高度取0.9，第2个初始相对高度取0.2，球的高度杂乱无章地变化。

(7)如M2.1g图所示，μ不变，第1个初始相对高度不变，第2个初始相对高度只增加10-5，以后的高度变化也迥然不同。

(8)如M2.1h图所示，μ不变，第1个初始相对高度不变，第2个初始相对高度稍大一点(10-8)，以后的高度变化也迥然不同。

这种对于初始条件十分敏感的运动称为混沌运动。

M2.1e图M2.1f图
M2.1g图M2.1h图
[解析](2)当参数连续变化时，同样利用(2.5)式计算高度。

当迭代次数n足够多的时候，对于周期性的不动点，x n就代表稳定值x∞。

取μ为自变量，取x=x∞为函数，可作μ-x曲线。

[算法](2)μ从0到4连续取值，先通过迭代算法筛去过渡值，继而用迭代算法获取迭代的结果，画出迭代图。

[程序]M2_2.m如下。

%罗基斯第模型的倍周期分岔和混沌图
clear%清除变量
x=0.2;%初始值(可任取)
u=0:0.0001:4;%参数向量(1)
n=1000;%迭代次数(2)
for i=0:n%按迭代次数循环
x=u.*(x-x.^2);%迭代计算消除暂态过程(3)
end%结束循环
figure%开创图形窗口
grid on%加网格
xlabel('\it\mu','FontSize',16)%横坐标
ylabel('\itx','FontSize',16)%纵坐标
title('罗基斯第模型的倍周期分岔和混沌图','FontSize',16)%标题
hold on%保持图像
cc='bgrk';%颜色代码(4)
for i=1:2^6%再按迭代次数循环
x=u.*(x-x.^2);%迭代一次
plot(u,x,['.'cc(mod(i,4)+1)],'MarkerSize',1)%画点图(5)
end%结束循环
[说明](1)参数向量的间隔很小，可当作连续分布的。

(2)进行1000次迭代，对于周期性运动，可筛去不稳定的点(相对高度)。

(3)变量x的初值是一个数值，第一次迭代之后就变成与u同样大小的向量。

x的每一个元素都代表u的对应元素的点(相对高度)，这就是用向量的好处。

(4)取4种颜色符号。

M2.2图
(5)在循环中画点时，对于周期运动，画出的同一高度；对于混沌运动，则画不同高度。

每循环4轮用同一颜色画点。

[图示]如M2.2图所示，当参数μ从0到1时，高度为零；当μ从0到3时，高度有一个不为零的值；当μ>3时，高度首先有两个值，然后分岔为4个值，再分岔为8个值，…，这种情况称为倍周期分岔；当μ达到某一值时，系统进入混沌状态。

混沌图还有复杂的结构。

[解析](3)在倍周期分岔中，分岔点划分了周期的范围。

设二元函数
f (μ,x )=μ(1–x )x (2.6)
对于周期1不动点，当n →∞时，有x n +1→x ∞，x n →x ∞，x ∞是不动点，用x 表示x ∞，可得
x =f (μ,x )=μ(1–x )x (2.7)
由此解得
x (1)=0，x (2)=1–1/μ(2.8)
不动点x (1)与参数无关，称为平凡不动点。

不动点x (2)与参数有关，称为本征不动点。

由于0≤x n ≤1，所以x (2)≥1。

函数对自变量的导数为
(,)(12)x f f x x x μμ∂=
=-∂(2.9)
不动点的稳定条件是|f x |<1
(2.10)对于平凡不动点，由于
f x (μ,x (1))=f x (μ,0)=μ(2.11)
可知：当μ<1时，x (1)是稳定不动点；当μ>1时，x (1)是不稳定的不动点，或者说x (1)失稳。

μ1=1是一个分岔点。

对于本征不动点，由于
(2)1(,)(,1)2x x f x f μμμμ=-=-(2.12)
只有满足-1<f x <1条件的点才是稳定的，所以当1<μ<3时，x (2)是稳定的不动点。

当μ>3时，f x <-1，x (2)失稳，因此μ2=3是一个分岔点。

分岔值为
x (2)=2/3(2.13)
对于周期2不动点，则有x n +2→x n →x ∞，用函数表示为f [μ,f (μ,x )]=x ，简记为
f 2(μ,x )=x (2.14)
指数表示函数嵌套。

根据(2.6)式可得
μ[1–μ(1–x )x ]μ(1–x )x =x
分解因式得
x [μx –(μ-1)][μ2x 2–(μ+1)μx +μ+1]=0(2.15)
方程除了x (1)和x (2)两个解之外还有两个解
(3)(4)1[12x μμ=+(2.16)
可见：当μ>3时，x 才有不相等的实数解，就是产生周期2的不动点。

不动点稳定条件是
|f x (μ,x (3))f x (μ,x (4))|<1(2.17)
设g (μ)=f x (μ,x (3))f x (μ,x (4))，利用(2.9)式可得
g (μ)=μ2(1–2x (1))(1–2x (2))=μ2[1–2(x (1)+x (2))+4x (1)x (2)]
再利用(2.16)式可得
g (μ)=-μ2+2μ+4
当g (μ)=-1时，解得1μ=g (μ)>-1可得
[(1(10
μμ--<
由于10μ>，必有
31μμ<=(2.18)
当g (μ)=1时，解得μ=3和-1。

由g (μ)<1可得μ>3。

因此，在μ2<μ<μ3范围内有周期2的稳定的不动点。

当μ>μ3时，x (3)和x (4)失稳，因此μ3是分岔点，分岔值为
x (3)=0.440，x (4)=0.8499(2.19)
[算法](3)当迭代次数n 足够多的时候，对于周期性的不动点，x n 就代表稳定值x ∞。

μ从
0到μ3连续取值，通过迭代算法最后获得稳定点x ，可与解析解进行比较。

μ从μ3到4连续取值，先通过迭代算法筛去过渡值，继续用迭代算法获取迭代的结果。

[程序]M2_3.m 如下。

%罗基斯第模型的倍周期分岔图
clear %清除变量
x=0.2;%初始值(可任取)
n=1000;%迭代次数
fs=16;%字体大小
u3=1+sqrt(6);%分岔参数(1)
u=linspace(0,u3);%参数向量
for i=0:n %按迭代次数循环
x=u.*(x-x.^2);%迭代计算形成不动点
end %结束循环
figure %开创图形窗口
plot([u,u],[x,u.*(x-x.^2)],'.')%画两支迭代曲线(点)
hold on %保持图像
plot([0,1],[0,0],'r')%画平凡不动点
u=1:0.1:3;%分岔前的参数向量
x1=1-1./u;%分岔前的不动点
plot(u,x1,'r')%画曲线
u=linspace(3,u3);%分岔后的参数向量
d=(u+1).*(u-3);%根的判别式
x21=(1+u+sqrt(d))./(2*u);%分岔后的上支不动点
x22=(1+u-sqrt(d))./(2*u);%分岔后的下支不动点
plot(u,x21,'r',u,x22,'r')%画2倍分岔后的曲线
grid on %加网格
xlabel('\it\mu','FontSize',fs)%横坐标
ylabel('\itx','FontSize',fs)%纵坐标
title('罗基斯第模型2倍周期分岔前后的不动点','FontSize',fs)%标题
legend('迭代点','解析曲线',2)%图例
axis([0,u3,0,1])%坐标范围
x=0.2;%初始值(可任取)
u=3.401:0.0001:4;
%参数向量(2)for i=0:n
%按迭代次数循环x=u.*(x-x.^2);
%迭代计算消除暂态过程end
%结束循环figure
%开创图形窗口grid on
%加网格xlabel('\it\mu','FontSize',fs)
%横坐标ylabel('\itx','FontSize',fs)
%纵坐标title('罗基斯第模型的倍周期分岔和混沌图','FontSize',fs)%标题hold on
%保持图像cc='bgrk';
%颜色代码for i=1:2^6
%再按迭代次数循环x=u.*(x-x.^2);
%迭代一次plot(u,x,['.'cc(mod(i,4)+1)],'markersize',1)%画点图
end %结束循环[说明](1)平凡解和二倍周期分岔的范围比较宽，单独画出来比较理想。

(2)多倍周期分岔和混沌的范围比较小，单独画出来比较理想。

[图示](1)如M2.3a 图所示，用迭代法和用解析法计算的结果完全相同。

在μ1=1处发生一次分岔。

在μ2=3处发生2倍周期分岔，这是树枝分岔。

(2)如M2.3b 图所示，随着参数的增加，在2倍周期分岔后又发生4倍周期分岔，8倍周期分岔，…，这种连续的分岔又叫做费根包姆分岔。

当周期无限增加时，就失去周期性，倍周期分岔就走向混沌。

混沌之中又有分岔，最明显的是周期5和周期3分岔，每个分支又通过倍周期分岔重新走向混沌。

在μ≈3.678处，两个大分支开始相遇。

相遇之前的上支与下支相似，两支又与总体相似。

不仅如此，更小的局部分支都与整体相似。

这种相似称为自相似。

M2.3a 图
M2.3b 图
[解析](4)对于一维映射x n +1=f (x n )
(2.20)可用初值x 0和附近值x 0+δx 0来计算分离。

作一次迭代后，距离为
010000d ()δ(δ)()δd f x x f x x f x x x =+-=
(2.21)
再作一次迭代后，距离为
121111100d ()δ(δ)()δ()()δd f x x f x x f x x f x f x x x
''=+-=
=(2.22)经过m 次迭代后，距离为1
00δδ()
m m n n x x f x -='=∏(2.23)
令0δδexp()
m m x x m λ=(2.24)则得1=1
0δ11ln(ln|()|δm m m n n x f x m x m λ-'==∑(2.25)李雅普诺夫指数为1=11lim ln|()|m n m n f x m λ-→∞'=∑(2.26)
当λ<0时，表示轨道在局部是稳定的，代表周期运动；当λ>0时，表示轨道分离，即对值具有敏感性，代表混沌。

当λ由负值变为正值时，表示周期性运动转向混沌。

对于罗基斯蒂模型，由于
x n +1=μ(1–x n )x n ，(n =0,1,2,…)
即f (μ,x n )=μ(1–x n )x n (2.27)所以(,)(12)n n n f x x x μμ∂=-∂(2.28)
[程序]M2_4.m 如下。

%罗基斯第模型的李雅普洛夫指数
clear
%清除变量x=0.1;
%初始值(0到1之间可任取)n=1000;
%迭代次数u=linspace(0,4,1000);
%参数向量(1)s=log(abs(u*(1-2*x)));
%第一个对数值for i=0:n
%按迭代次数循环x=u.*(x-x.^2);
%迭代计算s=s+log(abs(u.*(1-2*x)));
%累加导数的对数的绝对值(2)end
%结束循环lambda=s/n;
%李雅普洛夫指数figure
%开创图形窗口plot(u,lambda)
%画曲线grid on
%加网格fs=16;
%字体大小xlabel('\it\mu','FontSize',fs)
%横坐标ylabel('\it\lambda','FontSize',fs)%纵坐标
title('罗基斯第模型的李雅普洛夫指数','FontSize',fs)%标题[说明](1)修改向量范围，可画最后一段的指数。

u=linspace(3.5,4,10000);%参数向量
(2)在循环中求对数的和。

M2.4a图M2.4b图。