第一章等切面曲线和相似曲线 Soll:=p(t)=2 arctan/sinh(1 utV/R2=@.So11); Sol1 subs (abs(a"2-R"2)=R"2-a"2 2 (1+sinh(2))(a-P) Toy_Plot subs(Rx, Ry, Child_Subst, Soll, omega =1 a=4.95,R=5,[Toy,t=0.40]) plot(Toy_Plot, scaling constrained, color black) 图111给出此轨迹 图111a<R时玩具的轨道 注意,如果小孩在直线上行走,则有一个解析解:等切面曲线,如本章第一节所示.令(X0,YO) 是小孩的初始位置,(Vx,V)是他的恒定速度向量: Child_Subst : x(t)=Xo Vx*t, Y(t)=Yo Vy*t: Das combine(subs(Child_Subst, Daxial)); Das cos(φ(t))vy+sin(φ(t)) 1 dsolve(das, phi(t)) 2 Vy tan( p(t))+2 Vr Vy"+ v 重复前面的小孩在圆上行走的情况,我们又得到φ(t)的一个隐函数 总之,对于玩具的轨迹,不太可能得到一个解析解.如果考虑小孩在正弦曲线上行走,则只能 数值地计算玩具的轨迹 >X:=t→>t:Y:=t->5*sin(t):a:=10: axial p(t)=- cos((t))cos(t)+io sin(o(t)) F:= dsolve(Daxial, phi(0)=Pi/2, phi(t), numeric); c(rkf45-c)
Gander, s. Barton and. Hrebicek 因为此图形与图1.4相同,所以此处没有给出.然而我们将增加一些命令,以便在屏幕上动态地显 示这些运动.注意,这里还无法用 animate命令 我们想看棒的运动,以及它两端的运动轨迹.一端必须沿正弦曲线移动,另一端将描述所计 算的轨迹我们希望看到,在T秒内绘出N个点序列过程的片段 >T;=6*Pi:N:=200:L:a; SF Levalf(seq(rhs((Ti/N)[2]),i=0.N))] ST Levalf(seq (T*i/N, i=0.N)) plot([seq([X(ST [i])+L*cos(SF[i]), Y(sT[i])+L*sin(SF [i]),i =1.N)]) 上面最后的 MAPLE命令将显示此轨迹.为了此片段,我们必须准备和存储绘图点列,该点列是由 递增数给出的 TS: [seq(op(1, plot( [seq( [X(ST[i])+L+cos(SF[i]) Y(sT[i])+L*8in(SF[i])],i=1..j)])),j=1.N)]: BS :E [seq (op(1, plot([[X(sT [i]),Y(sT[i])] [X(ST[i])+L+cos(SF[i]), Y(ST[i])+L*sin(SF[i]]])) i=1.N)] PS =[seq(op(1, plot( [seq([X(ST[i]), Y(ST[i]) i=1..j)])),j=1.N)] 序列TS包含了轨迹的信息,它存储从1到N的N个部分图形点.序列BS描述棒的位置 PS描述正弦轨迹,它的结构类似于TS 现在用命令 PLOT (ANIMATE(seq([TS[i], BS[i], PS[i]], i=1.N)) AXESLABELS(x, y), VIEW(-2.X(ST [N]), DEFAULT)); 可看到整个过程 参考文献 E. HAIRER, S.P. NORSETT AND G. WANNER, Solving Ordinary Diferential Equations I, Springer-verlag Berlin Heidelberg, 1987 [2] H. HEUSER, Gewohnliche Diferentialgleichungen, B. G. Teubner, Stuttgart, 1989
第二章旋转网球的轨迹 F. Kwania 1引言 考虑一个质量为m,直径为d的网球,在靠近地球表面的空中运动.此球以角速度d旋转(向 量“有旋转轴的方向和大小=dφp(t)/dt=φ(t),其中φ(t)是旋转的角度)我们将笛卡儿坐标系 (xyz)放在地球的表面,其中z轴为垂直方向 图21旋转网球在空气中的移动 我们把网球看成一个质点,它在下列一些力的作用下运动(参见图21) 重力d=mg,其中=(0,0,-9)是重力加速度向量 拉力D=-D1()/,它与速度的方向相反 Magnus力M=ML古/x/v;此力垂直于和 拉力的大小D(v)和 Magnus力的大小ML(v),通常假设由理想的流体理论[】给出: 1丌d2 DE(u)=Cp 24 其中ρ是空气密度,对于实际的流体(空气),系数CD和CM依赖于速度υ,球的旋转和它的表面 材料.通常我们由实验得到这些系数 2]报告了旋转网球的实验结果.实验结果表明,对于一个邀度v∈[136,28]m8-和球的转 数n∈{800,32507m的网球,系数CD和CM只依赖于v/w,其中v=d/2.在某种意义上
d×/v|是转球赤道的速度ud/2在速度向量d上的投影.下面是系数的表达式 0508+ 22053+4196() 022+0.981( 对于网球,我们可以忽略球旋转的减速度,因此w是常数 于是,由位置向量r(t)的 Newton方程,可定义球的轨道 d'r Dr -+ML-x 其中初始条件为 d r 和一x(0) 10 方程(23)是三个微分方程的非线性系统,不存在解析解,因此我们必须数值地求解 实际上,大多数情况是高速旋转的高球,角速度向量位于在水平面内,并且垂直于向量v0,由 方程(2.3),对于t≥0,它也垂直于U(t),因此轨道在垂直平面内.选取x轴在此平面内,于是方 程(23)的最后形式是 CDaυi+mC g-CDaυz-nCMαUx 其中υ=√2+2和a=(pmd)/(8m)参数η=±1描述旋转的方向(对于高速旋转,n=1) 对于t=0,我们取初始条件 r(0)=0,z(0)=h,(0)=0cos(0),(0)=osin() 其中v是初始速度向量vo的大小,是v与x轴之间的夹角 22 MAPLE解法 为了说明拉力和 Magnus力对球的轨迹的影响,我们将考虑三种模型:在真空中的球,在空气 中没有旋转的球,和在空气中旋转的球 在真空中的球 在真空中,只有重力作用在球上,方程(24)有一个非常简单的形式 MAPLE用函数 dsolve求解(2.6)的通解 motion of ball in vacuum eqnid diff(x(t), t$2)=0, diff(z(t),t$2)=-g varid =(x(t), z(t)1 initcid =x(o)=0, z(0)=h, D(x)(O)=vO*cos(theta), D(z)(O)= vo*sin(theta):
第二章旋转网球的轨迹 #solution for ball in vacuum resid: dsolvedfeqnid, initcid, varid resid x()=tv0co(,2()=h+to0n()-29 注意,函数 dsolve以方程组的形式 varname)=erpression 返回解( result),因此没有定义这些方程的顺序,为了访问对应( varname)的解,我们将用函数 subs((result), (varname)) 在空气中的球 对于在空气中的球,我们不能解析地求解方程(24),必须用数值解法因为后面要用到 Newton 法计算飞行时间,我们需要z(t)的导数,因此将二阶微分方程组(24)转变成一阶微分方程组 (用变量vx和vz分别表示速度北和, CLau·vx+ nEma·U (27) g-CLaU·vz- nCMaU 其中v=√+v2.初始条件(2.5)为 a(0)=0, z(0)=h, v=(0)=vo cos(1), vz(0)= vo sin(v) 我们将用SI单位系统,参数的数值是g=981{ms-l],网球的直径d=0063m,其质量m 005kg,空气密度p=1.29kgm-3],我们选取初始条件 h=1[m,o=25[m8-1],=15° 和球旋转的参数u=20{ms-1和η=1(高速旋转),用 MAPLE语句求解此问题1如下所示 numeric solution of ball in air Basic constants in si units >g:=9.81: gravitational acceleration d::=0.063: m =0.05: diameter and mass of ball rov = 1. 29: density of air > alpha: evalf(Pi*d"2/(8*m)rov): >v: =(dx"2+ dz"2)"(1/2): definition of velocity cd:=.508+1/(22.503+4.196/(w/v(t))(5/2))-(2/5) >Cm:=eta/(2.202+.981/(w/v(t))) eta =+-1 defines direction of rotation var =x(t), z(t),dx(t),dz(t) initial conditions 1此例子已用 MAPLE V Release3测试过