自行车轮饰物的运动轨迹问题
理工大学暑期数学建模强化训练专题四
自行车轮饰物 的运动轨迹问题
学员： 曹 阳
倪迪杭
学院： 通信工程学院
时间：2010.08.19
自行车轮饰物的运动轨迹问题
摘 要
本文就自行车轮饰物的运动轨迹问题，采用解析几何的方法建立数学模型，求出了自行车在各种不同形状的道路上行驶时饰物和椭圆板中心的运动轨迹方程，并且利用Matlab 软件模拟仿真出了两者的运动轨迹。
对于问题1和问题2，先运用解析几何方法求出自行车轮轴心的轨迹方程，而后利用饰物始终绕车轮轴心作圆周运动建立参数方程，求出饰物的轨迹方程。求出的曲线轨迹分别见图2、图4和图5。 对于问题4，将“圆板”换为“椭圆板”，通过设定参数，结合坐标转换的知识，将转动过程中椭圆板中心的坐标用该参数表示，求出了其运动轨迹的参数方程。其轨迹图像见图8。
关键词:运动轨迹，解析几何，抛物线，椭圆，坐标转换
一、问题的提出
为了改变平淡的自行车外表，给自行车添加一分美妙的动感，同时，也为了增加骑车人的“安全系数”，一些骑车人及自行车厂家在自行车的辐条上安装一块亮丽的饰物。当有这种饰物的自行车在马路上驶过时，这种饰物就如游龙一样，对街边的行人闪过一道波浪形的轨迹。这一波一闪的光亮游龙，也默默地维护着骑车人的安全。建立数学模型解决以下问题：
1、这轨迹是什么曲线？试画出它的图形。
2、当这自行车又在一个抛物线形的拱桥上通过时，或是在一拱一拱的正弦曲线（例如山地摩托车赛场）上通过时，这饰物又画出一条曲中有曲的轨迹，这轨迹是什么曲线？试画出它的图形。
3、这种滚动中圆盘中心的运动轨迹是什么？ 4、将问题中“圆板”换为“凸形板”（例如椭圆板）时，其滚动轨迹会有什么结果？
二、问题的分析
对于问题1、2，装有饰物的自行车在马路上行驶过程中，饰物会形成一道道曲线轨迹。而随着路况的不同，如平坦的公路、抛物线形的拱桥、正弦曲线形的山地摩托车赛场等，饰物会形成不同的曲线轨迹。不管形状的道路怎么样，不管曲线轨迹有多复杂，都可以先运用解析几何方法求出自行车轮轴心的轨迹方程，而后利用饰物始终绕车轮轴心作圆周运动建立参数方程，求出饰物的轨迹方程。饰物绕车轮轴心的参数方程是较易得到的，故问题的关键在于求出不同形状的道路上自行车轮轴心的轨迹方程。求出饰物的轨迹方程后，设定合理的参数，利用Matlab 软件就可以模拟仿真出饰物的曲线轨迹。
对于问题4，研究凸形板（椭圆板）中心的运动轨迹时，可以通过设定参数，将转动过程中椭圆板中心的坐标通过坐标转换用参数表示，求出其运动轨迹方程，设定合理的参数，利用Matlab 软件就可以模拟仿真出椭圆板中心的曲线轨迹。
三、模型假设
1、自行车在行驶过程中车轮不打滑； 2、自行车在行驶过程中速度保持不变；
3、自行车在行驶过程中车轮始终与地面接触。
四、模型的建立及求解
1、符号说明
R 自行车轮的半径
r 饰物距车轮轴心的距离
B 自行车车轮上的饰物
自行车轮的转动速度
a , b 定常量
2、模型建立
2.1 模型1（问题1） 2.1.1模型的分析
自行车在马路上行驶时，可以认为马路是平坦的，因此自行车轮轴心的运动轨迹就是一条直线。求出这条直线方程后联立饰物绕车轮轴心的参数方程，就可以容易的得到饰物的运动轨迹。在检验模型的正确性时，可以利用图像求出几个特殊点的坐标，与利用轨迹方程求出的坐标一一比较，若两者求出的坐标都是相同的，则验证了模型的正确性。
2.1.2模型的建立与求解
设自行车轮的半径为R ，车轮上的饰物B 距离车轮轴心为r 。以自行车轮上某一点为坐标原点建立坐标系，假设饰物B 的初始位置在y 轴上，其示意图如图1所示，则饰物B 的初始坐标B 0为(0, R -r )。
图1 自行车轮示意图
设自行车轮的转动速度为ω，则自行车的速度为v =ωR 。经过时间t 后，饰物B 的绕自行车轮轴心旋转需满足参数方程：
⎧x =-r sin ωt ，
⎨
y =-r cos ωt 。⎩
而自行车轮轴心的运动轨迹需满足参数方程：
⎧x =ωRt ，
⎨
y =R 。⎩
所以，饰物B 的坐标应满足参数方程：
⎧x =-r sin ωt +ωRt ，
⎨
⎩y =-r cos ωt +R 。
根据一般自行车的规格，取R =0.5m ，r =0.4m ，ω=20rad /s ，使用Matlab 软件画出其图像如图2所示：
2.1.3模型的验证
用饰物B 在特殊点的坐标来验证模型。假设自行车沿正x 轴方向行驶，则自
π3
行车轮顺时针转动。饰物B 转动、π和π时的坐标分别为B 1、B 2和B 3。
22
π
饰物B 转动时的图像如图3所示。
2
图2 饰物在平坦道路上的轨迹
π
图3 B 转动时的图像
2
π
时，车轮中心经过的2
πππ
距离为R ，故B 1的横坐标为R -r 。因此饰物B 转动时的坐标为
222
由图像我们可以得到B 1的纵坐标为R 。饰物B 转动
T π⎛π⎫⎛π⎫
带入方程式求得的B 的坐标B 1' R , R -r ⎪，这B 1 R , R -r ⎪。而把t ==
42ω⎝2⎭⎝2⎭
与B 1的坐标是相同的。
3
用同样的方法，我们可以把饰物B 转动π和π时的图像作出，然后利用图
2
33
像求得饰物B 转动π和π时的坐标，这与由方程式计算出饰物B 转动π和π
22
时的坐标都是相同的。从而验证了我们的模型的正确性。
2.2 模型2（问题2）
2.2.1模型的分析
当这自行车在一个抛物线形的拱桥上通过时，自行车轮轴心的轨迹是一个类抛物线。在求饰物的轨迹时，先利用拱桥的抛物线方程求出自行车轮轴心的轨迹，然后联立饰物的坐标的参数方程，就可以得到饰物的轨迹。
当这自行车在一拱一拱的正弦曲线上通过时，使用类似的方法，先求出自行车轮轴心的轨迹，然后确定饰物的轨迹。
2.2.2模型的建立与求解 （1）当这自行车在一个抛物线形的拱桥上通过时，设拱桥的抛物线方程为：
⎧⎪x 1=x 0,
⎨2
⎪⎩y 1=ax 0+bx 0(a 0) 。
抛物线上任意一点的切线斜率k 1为：
k 1=y 1' |x =x 0=2ax 0+b ，
则该点的法向量的斜率k 2为：
k 2=-
1
。 k 1
假设自行车轮的轴心坐标为(x 2, y 2) ，则自行车轮的轴心的轨迹方程为：
⎧(y 2-y 1) 2+(x 2-x 1) 2=R 2, ⎪
⎨y 2-y 1
⎪x -x =k 2。⎩21
解之得：
⎧
x =x 1⎪2
⎪⎪
⎨
⎪y =y +21
⎪⎪⎩
这就是自行车轮轴心的轨迹方程。
抛物线上从原点到任意一点的距离为s ，则有：
s =⎰=
x 00
1⎡
k ln k 1+⎤
⎥⎣⎦4a ⎢
1⎡-ln b ⎤。
⎥⎣⎦4a ⎢
(
(
又根据行驶路程与行驶时间之间的关系，有
s =vt =ωRt ，
令u =ωt ，则：
s u =。
R
而经过时间t 后，饰物B 的绕自行车轮轴心旋转需满足参数方程：
⎧x =-r sin u ，
⎨
⎩y =-r cos u 。
而自行车轮轴心的运动轨迹需满足参数方程：
⎧x =x 2，
⎨y =y 。⎩2
所以，饰物B 的坐标应满足参数方程：
⎧x =-r sin u +x 2，
⎨
⎩y =-r cos u +y 2。
这就是当这自行车在一个抛物线形的拱桥上通过时饰物B 的运动轨迹方程。 取a =-0.05，b =2，R =0.5m ，r =0.4m ，使用Matlab 软件画出其图像如图4所示：
图4 饰物在抛物线拱桥上的轨迹
（2）当这自行车在一拱一拱的正弦曲线上通过时，设正弦曲线的方程为：
⎧x 1=x 0,
⎨
y =a sin bx (a >0, b >0) 。0⎩1
正弦曲线上任意一点的切线斜率k 1为：
k 1=y 1' |x =x 0=ab cos x 0，
则该点的法向量的斜率k 2为：
k 2=-
1
。 k 1
假设自行车轮的轴心坐标为(x 2, y 2) ，则自行车轮的轴心的轨迹方程为：
⎧(y 2-y 1) 2+(x 2-x 1) 2=R 2, ⎪
⎨y 2-y 1
=k 。2⎪x -x
⎩21
解之得：
⎧
x =x 1⎪2
⎪⎪
⎨
⎪y =y +1
⎪2⎪
⎩
这就是自行车轮轴心的轨迹方程。
设自行车在正弦曲线上行驶的距离为s ，则有：
s =⎰
x 0
x 0'
。
这个定积分较为复杂，故在求解s 的时候，运用微积分的思想，将图形等步长微分成n 份，对每一部分求出面积后累加即可得到s 的解。n 越大，求出的s 就却精确。故有：
s =⎰
x 0
'
=i =1
n
x 0
同样令u =ωt ，根据s =vt =ωRt ，得到饰物B 的坐标应满足参数方程：
⎧x =-r sin u +x 2，
⎨
y =-r cos u +y 。⎩2
这就是当这自行车在一拱一拱的正弦曲线上通过时饰物B 的运动轨迹方程。
取a =1，b =0.6，n =100，R =0.5m ，r =0.4m ，使用Matlab 软件画出其图像如图5所示：
图5 饰物在正弦曲线上的轨迹
2.3 模型3（问题4） 2.3.1模型的建立与求解
将问题中的“圆板”换为“椭圆板”之后，考虑在水平道路上运动。 设初始状态时椭圆的方程为：
(y -a )
a 2
其参数形式为：
2
x 2
+2=1(a >b >0)， b
⎧x =b sin θ,
⎨
⎩y =a -a cos θ。
则椭圆中心的坐标为O (0, a )，坐标原点为O 1(0,0)，如图6所示。
设椭圆转动任意角度μ后，椭圆中心为O ' (x 2, y 2)，椭圆与x 轴的切点为A (x 0, y 0)。尝试通过坐标转换，找出椭圆转动角度μ后x 2和y 2关于μ的参数表达式，即：
⎧⎪x 2=x 2(μ),
⎨
⎪⎩y 2=y 2(μ)。
图6
则过切点A (x 0, y 0)的切线方程为：
(y 0-a )(y -a )+x 0 x =1。
a 2
b 2
所以，该切线的斜率k 2为：
a 2x 0
k 2=-2=tan μ。
b y 0-a 从而可以推得：
a 2
y 0-a =-2x 0。 （1）
b k 2
点A (x 0, y 0)在椭圆上，故又满足：
(y 0-a )
a 2
联立（1）、（2）可以得到：
2
x 02
+2=1。 （2）
b
⎧2⎪x 0=⎪
⎨2
⎪y =a ⎪0
⎩
这样就得到了x 0和y 0关于μ的参数表达式。
由于点A (x 0, y 0)在椭圆上，亦满足椭圆的参数方程：
⎧x 0=b sin θ0,
⎨
y =a -a cos θ。0⎩0
故得到：
θ0=arc cos 1-
⎛
⎝y 0⎫⎪。 a ⎭
令直线OA 与椭圆长轴的夹角为α，则有：
tan α=
x 0
。 a -y 0
故，
α=arc tan
x 0
。 a -y 0
A 为s ，则有：
令坐标原点O 1和切点A 之间的弧长O 1
s =⎰
θ0
θ。
运用微积分的思想求解s ，将图形等步长微分成n 份，对每一部分求出面积
后累加即可得到s 的解。n 越大，求出的s 就却精确。故有：
s =⎰
θ0
θ
=i =1
n
切点A (x 0, y 0)到椭圆中心O (0, a )的距离
d =。 观察图7，容易得到椭圆中心O ' 在坐标轴中满足：
⎧x 2=O 1A ' +O ' A ' sin (μ-α),
⎪
⎨' '
⎪⎩y 2=O A cos (μ-α)。
所以，椭圆中心O ' 在椭圆转动过程当中的轨迹方程为：
⎧⎪x 2=s +d sin (μ-α),
⎨
⎪⎩y 2=d cos (μ-α)。
取a =4.5，b =2.5，使用Matlab 软件画出其图像如图8所示：
图7
图8 椭圆中心在水平道路上运动轨迹
五、模型的评价
本文通过对自行车运动过程中饰物以及圆板（凸形板）的运动进行了详细的分析，并通过一定的物理和数学思想及方法，借助于MATLAB 软件模拟了问题中所需的各点的轨迹。在解决问题二时，我们巧妙地运用了数值积分的方法，把自行车在抛物线轨道上通过时饰物的轨迹方程求解出来，从而模拟出了其轨迹。对于问题四，我们设定了椭圆参数方程，通过对其运动过程物理状态的分析，完成了对其中心轨迹的求解。
本文的缺点在于由于时间有限，对于问题4的求解，只考虑椭圆板在水平直线上的运动。
七、模型的推广
可以将问题4拓展为较为简单的凸形板在常见平滑曲线上的运动，通过设定参数，结合坐标转换的知识，将转动过程中简单的凸形板中心的坐标用该参数表示，从而求出其运动轨迹的参数方程。
[1].同济大学数学系年4月出版
参考文献
第六版 上册》高等教育出版社2007 《高等数学
附录一
图2的程序：
t=0:0.01:3; r=0.4;R=0.5; a=3; v=a*t w=v./R
x =-r*sin(w.*t)+w.*R.*t ; y =R-r*cos(w.*t); plot(x,y)
axis([0 7 0 1.2])
附录二
图4的程序：
x0=-1:0.01:40; a=-0.05; b=2; R=0.5; r=0.4; x1=x0;
y1=a*x0.^2+b*x0 k1=2*a*x0+b; plot(x1,y1) hold on
k2=-1./k1;%法线的斜率;
x2=R./(k2.*sqrt(1+(1./k2).^2))+x1; y2=R./sqrt(1+(1./k2).^2)+y1 %plot(x2,y2) hold on
u=1/R*(1/(4*a)*(k1.*sqrt(1+k1.^2)+log(k1+sqrt(k1.^2+1)))-1/(4*a)*(b*sqrt(1+b^2)+log(b+sqrt(1+b^2))));
x=-r*sin(u)+x2; y=y2-r*cos(u); plot(x,y);
axis([-2 45 0 30])
图5程序：
x0=-2:0.01:18; a=1; b=0.6; R=0.5;
r=0.4; x1=x0;
y1=a*sin(b*x0); k1=a*b*cos(b*x0); plot(x1,y1) hold on
k2=-1./k1;%法线的斜率;
x2=R./(k2.*sqrt(1+(1./k2).^2))+x1; y2=R./sqrt(1+(1./k2).^2)+y1; %plot(x2,y2) hold on
for x3=-2:0.01:18 s=0; for i=1:100
s=s+sqrt(1+(a*b*cos(b*i*(x3+2)/100-2))^2)*(x3+2)/100 end u=s/R
y1=a*sin(b*x3); k1=a*b*cos(b*x3); k2=-1./k1;
x2=R./(k2.*sqrt(1+(1./k2).^2))+x3; y2=R./sqrt(1+(1./k2).^2)+y1; x=-r*sin(u)+x2; y=y2-r*cos(u); plot(x,y); hold on end
axis([-3 14 -2.5 2.5])
附录三
图8的程序：
clc
for u=0:0.001:1/2*pi a=4.5; b=2.5; k2=tan(u);
x0=(b^2*k2)/(sqrt(a^2+k2^2*b^2)) y0=a-a^2/sqrt(a^2+k2^2*b^2); d=sqrt(x0^2+(y0-a)^2); w=asin(x0./b) n=100; s=0; for i=1:n
s=s+w/n*sqrt(a^2*(sin(w*i/n))^2+b^2*(cos(w*i/n))^2); end
x2=s+d.*sin(u-w);
y2=d.*cos(u-w); plot(x2,y2) hold on end hold on
for u=0.5*pi:0.001:1*pi k2=tan(u);
x0=(b^2*k2)/(sqrt(a^2+k2^2*b^2)) y0=a-a^2/sqrt(a^2+k2^2*b^2); d=sqrt(x0^2+(y0-a)^2); w=pi+asin(x0./b) n=100;s=0; for i=1:n
s=s+w/n*sqrt(a^2*(sin(w*i/n))^2+b^2*(cos(w*i/n))^2); end
x2=s+d.*sin(u-w); y2=d.*cos(u-w); plot(x2,y2) hold on end
axis([0 15 2 5])