《理论力学》课程教学资源（文献资料）Matlab求解理论力学问题系列（四）乒乓球的滚动_高云峰

力学与实践第43卷第5期2021年10月Matlab求解理论力学问题系列(四)乒乓球的滚动高云峰1)(清华大学航天航空学院，北京100084)摘要本文讨论了乒乓球运动的理论解和数值解。在考虑摩擦、滚动摩阻、空气阻力等条件的一般情况下，只能获得部分解析解，而采用数值计算的方法可获得系统所有的运动规律。文章介绍了如何利用Matlab处理分段函数问题，包括求解微分代数方程、利用积分中断条件确定分段的时间点，以及利用部分解析解验证数值计算的准确性等，最后通过大量计算获得乒乓球向前运动及向后运动的分界线，以及乒乓球正好回到出发位置的条件。关键词微分代数方程，解析解，数值解，积分中断中图分类号：0313.3文献标识码：Adoi：10.6052/1000-0879-20-240THEROLL OF A TABLE TENNISGAO Yunfengl)(School of Aerospace Engineering, Tsinghua University, Beijing 100084, China)Abstract Theoretical and numerical solutions of the kinematics of table tennis are discussed in this paper.Under the general conditions of friction (sliding, rolling friction and air resistance), analytical solutions carbe obtained only in limited conditions,but numerical solutions are universally available.This paper handlesvariousnumerical problems in solving thekinematicswith Matlab,including using piecewisefunction,solvingdifferential-algebraic equation, obtaining the interruption point of the integral equation, and verifying the accuracy of the numerical results through comparing with available analytical solutions.At last, two criticalconditions are investigated, one for the determination of forward or backward motion of the ball, the other forthe ball jumping back to the starting position.Key words differential-algebraic equation, analytical solution, numerical solution, integral interruption可能很多人尝试过用手指按压乒乓球，在一定所有的解析解，这也是本系列文章的一个目的：绝大条件下可以使它向前运动，然后又会返回（图1）。本部分动力学问题没有解析解，但是可以利用数值计文从理论分析和数值计算的角度，对这一现象进行算的方法获得系统所有参数变化的规律，加深对问分析。题的理解。乒乓球滚动问题的特点是，乒乓球可能打滑也1、乒乓球滚动的动力学方程及部分解析解可能不打滑，摩擦力或滚动摩阻的方向也可能发生变化，因此需要在不同阶段列写不同的方程再进行假设乒乓球质量为m，半径为r，水平面摩擦处理。系数为μ，滚动摩阻为，空气阻力系数为n。乒乓如果同时考虑空气阻力或滚动摩阻，不能获得球的运动和受力如图2所示，其中c为乒乓球质心2020-06-01收到第1稿，2020-07-01收到修改稿。1)E-mail:gaoyunfeng@tsinghua.edu.cn引用格式：高云峰。Matlab求解理论力学问题系列（四）乒乓球的滚动.力学与实践，2021，43(5)：753-759Gao Yunfeng. The roll of a table tennis. Mechanics in Engineering, 2021, 43(5): 753-759(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 43 卷 第 5 期 力 学 与 实 践 2021 年 10 月 Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 高云峰1) (清华大学航天航空学院, 北京 100084) 摘要 本文讨论了乒乓球运动的理论解和数值解。在考虑摩擦、滚动摩阻、空气阻力等条件的一般情况下， 只能获得部分解析解，而采用数值计算的方法可获得系统所有的运动规律。文章介绍了如何利用 Matlab 处理 分段函数问题，包括求解微分代数方程、利用积分中断条件确定分段的时间点，以及利用部分解析解验证数值 计算的准确性等，最后通过大量计算获得乒乓球向前运动及向后运动的分界线，以及乒乓球正好回到出发位置 的条件。 关键词 微分代数方程, 解析解, 数值解, 积分中断 中图分类号: O313.3 文献标识码: A doi: 10.6052/1000-0879-20-240 THE ROLL OF A TABLE TENNIS GAO Yunfeng1) (School of Aerospace Engineering, Tsinghua University, Beijing 100084, China) Abstract Theoretical and numerical solutions of the kinematics of table tennis are discussed in this paper. Under the general conditions of friction (sliding, rolling friction and air resistance), analytical solutions can be obtained only in limited conditions, but numerical solutions are universally available. This paper handles various numerical problems in solving the kinematics with Matlab, including using piecewise function，solving differential-algebraic equation, obtaining the interruption point of the integral equation, and verifying the ac￾curacy of the numerical results through comparing with available analytical solutions. At last, two critical conditions are investigated, one for the determination of forward or backward motion of the ball, the other for the ball jumping back to the starting position. Key words differential-algebraic equation, analytical solution, numerical solution, integral interruption 可能很多人尝试过用手指按压乒乓球，在一定 条件下可以使它向前运动，然后又会返回 (图 1)。本 文从理论分析和数值计算的角度，对这一现象进行 分析。 乒乓球滚动问题的特点是，乒乓球可能打滑也 可能不打滑，摩擦力或滚动摩阻的方向也可能发生 变化，因此需要在不同阶段列写不同的方程再进行 处理。 如果同时考虑空气阻力或滚动摩阻，不能获得 所有的解析解，这也是本系列文章的一个目的：绝大 部分动力学问题没有解析解，但是可以利用数值计 算的方法获得系统所有参数变化的规律，加深对问 题的理解。 1 乒乓球滚动的动力学方程及部分解析解 假设乒乓球质量为 m，半径为 r，水平面摩擦 系数为 µ，滚动摩阻为 δ，空气阻力系数为 n。乒乓 球的运动和受力如图 2 所示，其中 xC 为乒乓球质心 2020–06–01 收到第 1 稿，2020–07–01 收到修改稿。 1) E-mail: gaoyunfeng@tsinghua.edu.cn 引用格式: 高云峰. Matlab 求解理论力学问题系列 (四) 乒乓球的滚动. 力学与实践, 2021, 43(5): 753-759 Gao Yunfeng. The roll of a table tennis. Mechanics in Engineering, 2021, 43(5): 753-759

754力学与实践2021年第43卷阶段，乒乓球打滑前进（图2），可以列出动力学方程用力按压及补充的摩擦关系，为mic=-F-ntc0=N-mg(1)图1乒乓球游戏Jco=-oN+FrF=μN+y其中Jc=2mr2/3是薄壁圆球对质心的转动惯量。RCC令β=n/(2m），方程（1）中的未知力消去后可以得mg到系统的运动微分方程，为M.FAt)A1c+2+μg=0NI(2)茶-(μr - 0)g/r2 = 0图2乒乓球运动和受力图方程（2）存在解析解，注意Matlab除了可以数的运动，e为乒乓球的转角，N为压力，F为摩擦值计算，还可以进行符号推导，包括对微分方程的积力，M=aN为滚动摩阻，R=nic为空气阻力。分[1]。利用Matlab对方程(2)进行积分，源代码见根据实际经验，乒乓球在初始阶段会打滑，后来图3。会进入纯滚动状态，因此下面分两阶段来处理。第一%定义符号syms xc(t) theta(t) beta mu g r delta geql=diff(xc,t,2)+2*beta*diff(xc,t,1)+mu*g；%位移的微分方程eq2=diff（theta，t，2)-3/2*（mu*r-delta）*g/r*2：%角度的微分方程xc=dsolve(eq)%位移的积分结果%角度的积分结果theta-dsolve(eq2)图3微分方程的积分度ur=ro，代入式(3)处理后有图3中syms为定义符号：dif(xc.t.2)表示xc对时间的2阶导数；dsolve表示对微分方程进行符-2Bt-0)gt/r20-9μ/B)(ur号积分。运行后屏幕上显示方程（2）的解为29μ/β-r0 = 0(4)xc=C1+C2*exp(-2*beta*t)-(g*mu*t)/(2*beta)theta=C4+C3*t-(t"2*(3*delta*g-3*g*mu*r))/设*时刻触点A速度为零，在退化情况下方程（4)(4*r*2)可以有解析解，但是在一般情况下没有*的解析表达式(这将导致后续分析出现问题)。不过在数值积代入初始条件后有分中t可以利用积分中断获得，这表明了在一般情 (gμ/β2 2io/B) (1 - e-2pt)-况下利用数值计算的潜力。TC:1第二阶段，乒乓球作纯滚动，但乒乓球可能向前-59μt/B(3)运动也可能向后运动，滚动摩阻的方向相应也应改变，定义方向函数6=(ur-8)gt2/r2+0ot{+1, if uc >0这一阶段结束的标志是接触点A速度为零。A点速(5)入=if vc=00.度根据图2可以分解为牵连速度Ue=c和相对速-1, if uc<0(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.Allrightsreserved.http://www.cnki.net

754 力 学 与 实 践 2021 年 第 43 卷 ঻᤹࣋⭘ 图 1 乒乓球游戏 xC xC vA Mf θ θ . . y O C A mg N F C A x R 图 2 乒乓球运动和受力图 的运动，θ 为乒乓球的转角，N 为压力，F 为摩擦 力，Mf = δN 为滚动摩阻，R = nx˙ C 为空气阻力。 根据实际经验，乒乓球在初始阶段会打滑，后来 会进入纯滚动状态，因此下面分两阶段来处理。第一 阶段，乒乓球打滑前进 (图 2)，可以列出动力学方程 及补充的摩擦关系，为 mx¨C = −F − nx˙ C 0 = N − mg JC ¨θ = −δN + F r F = µN    (1) 其中 JC = 2mr2/3 是薄壁圆球对质心的转动惯量。 令 β = n/(2m)，方程 (1) 中的未知力消去后可以得 到系统的运动微分方程，为 x¨C + 2βx˙ + µg = 0 ¨θ − 3 2 (µr − δ)g/r2 = 0    (2) 方程 (2) 存在解析解，注意 Matlab 除了可以数 值计算，还可以进行符号推导，包括对微分方程的积 分 [1]。利用 Matlab 对方程 (2) 进行积分，源代码见 图 3。 图 3 微分方程的积分 图 3 中 syms 为定义符号；diff(xc,t,2) 表示 xc 对时间的 2 阶导数；dsolve 表示对微分方程进行符 号积分。运行后屏幕上显示方程 (2) 的解为 xc=C1+C2*exp(-2*beta*t) - (g*mu*t)/(2*beta) theta=C4+C3*t-(tˆ2*(3*delta*g-3*g*mu*r))/ (4*rˆ2) 代入初始条件后有 xC = 1 4 ( gµ/β2 − 2 ˙x0/β) (1 − e −2βt)− 1 2 gµt/β θ = 3 4 (µr − δ)gt2/r2 + ˙θ0t    (3) 这一阶段结束的标志是接触点 A 速度为零。A 点速 度根据图 2 可以分解为牵连速度 ve = ˙xC 和相对速 度 vr = r ˙θ，代入式 (3) 处理后有 ( x˙ 0 − 1 2 gµ/β) e −2βt − 3 2 (µr − δ)gt/r − 1 2 gµ/β − r ˙θ0 = 0 (4) 设 t ∗ 时刻触点 A 速度为零，在退化情况下方程 (4) 可以有解析解，但是在一般情况下没有 t ∗ 的解析表 达式 (这将导致后续分析出现问题)。不过在数值积 分中 t ∗ 可以利用积分中断获得，这表明了在一般情 况下利用数值计算的潜力。 第二阶段，乒乓球作纯滚动，但乒乓球可能向前 运动也可能向后运动，滚动摩阻的方向相应也应改 变，定义方向函数 λ =    +1, if vC > 0 0, if vC = 0 −1, if vC < 0 (5)

第5期755高云峰：Matlab求解理论力学问题系列（四）乒乓球的滚动X=inv(A)*B求出当前时刻的和i(求解后已注意这一阶段A点速度为零，此时动力学方程及补充运动学关系为经是具体数值），再把动力学方程转换为标准的一阶方程组形式，设mic=-F-nic0=N-mg1=c,92=c,93=,94=(11)(6)Jc=-N+Fric=re从而得到一阶微分方程组和初始条件为y1(0) = 0方程（6)中消去未知的力，可以简化为j1 =y23/2(0) = c(0)92=c63(12)(7)=0J3(0) = 093 = 94094=6*ya(0) = 0(0)类似可以利用Matlab的符号推导得到解析解，为第一阶段和第二阶段的分界点时刻*可以由C = Ci + C2e-6βt/5 _ 入gt/(2βr)(8)Matlab积分时的中断条件获得。在积分之前，先把积分中具体的要求在options中说明，其格式为问题在于，由于第一阶段结束的时间t*以及ac(t*)和c(t*）没有解析表达式，方程（8）中的积分常数options =odeset('RelTol',numberl,'AbsTol,无法确定，因此无法得到更多的理论解。number2,'events',@event.name);(13)由于在一般情况下很难得到动力学问题的全部解析解，下面介绍数值计算。格式中黑正体表示固定格式，白斜体表示变量。注意“events，@表示在积分时会调用名为event_name2微分方程的分段准备工作这个子程序，具体的中断条件可以在该子程序中表首先把第一、二阶段的动力学方程(1)和(6)转达，在本问题中即A点速度为零。在设定积分条件换为AX=B的形式，分别为后，再调用积分函数，其格式为m010ic-ni[t,] = ode45(@rg-kt.name, [start.time :?00Jc-r(9)Fstep-time:end.timel,y0,options);(14)0100mgLN0L001-H格式中黑正体表示固定格式，白斜体表示变量。其0[m1-ni0[ic中[,是积分的时间和结果；动力学方程在子程序苗入o00Jc-rrg-kt-name中描述；y是初始条件：如果没有中断条(10)F0001mg件，积分时按等步长step_time从start.time开始到0-r0ON1end_time结束，有中断条件后，最后一步的积分步这是混合形式的微分代数方程。在数值计算中，长会自动调整，结束时间也会自动提前。当前t时刻的c,,c,是已知量，而c,,F,在程序中先给参数赋值，然后分段计算，这一部是未知量，因此先把c和é视为代数量，调用分的源代码见图4。options=odeset(RelTol',le-8,AbsTol',le-10,events,@event_pingpang1)[ttl,iyl]=ode45(@rg_kt_pingpangl,[tto:hh:tt0+10],y0,options);%计算积分的步数NUM-length(tt1):new_y0=iy1（NUM，：）；%把最后一步的参数设为初始值%最后一步的时间tf=tt1(NUM):options=odeset(RelTol',le-8,"AbsTol,le-10):[tt2,iy2]=ode45(rg_kt_pingpang2",[tf:hh:tf+10],new_y0,options);图4带中断条件的分段积分(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 5 期 高云峰：Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 755 注意这一阶段 A 点速度为零，此时动力学方程及补 充运动学关系为 mx¨C = −F − nx˙ C 0 = N − mg JC ¨θ = −λδN + F r x¨C = r ¨θ    (6) 方程 (6) 中消去未知的力，可以简化为 x¨C + 6 5 βx˙ C + 3 5 gλδ/r = 0 (7) 类似可以利用 Matlab 的符号推导得到解析解，为 xC = C1 + C2e −6βt/5 − λgt/(2βr) (8) 问题在于，由于第一阶段结束的时间 t ∗ 以及 xC (t ∗ ) 和 x˙ C (t ∗ ) 没有解析表达式，方程 (8) 中的积分常数 无法确定，因此无法得到更多的理论解。 由于在一般情况下很难得到动力学问题的全部 解析解，下面介绍数值计算。 2 微分方程的分段准备工作 首先把第一、二阶段的动力学方程 (1) 和 (6) 转 换为 AX = B 的形式，分别为       m 0 1 0 0 JC −r δ 0 0 0 1 0 0 1 −µ             x¨C ¨θ F N       =       −nx˙ 0 mg 0       (9)       m 0 1 0 0 JC −r λδ 0 0 0 1 1 −r 0 0             x¨C ¨θ F N       =       −nx˙ 0 mg 0       (10) 这是混合形式的微分代数方程。在数值计算中， 当前 t 时刻的 xC , θ, x˙ C , ˙θ 是已知量，而 x¨C , ¨θ, F, N 是未知量，因此先把 x¨C 和 ¨θ 视为代数量，调用 X =inv(A)*B 求出当前时刻的 x¨ ∗ C 和 ¨θ ∗ (求解后已 经是具体数值)，再把动力学方程转换为标准的一阶 方程组形式，设 y1 = xC , y2 = ˙xC , y3 = θ, y4 = ˙θ (11) 从而得到一阶微分方程组和初始条件为 y˙1 = y2 y˙2 = ¨x ∗ C y˙3 = y4 y˙4 = ¨θ ∗ , y1(0) = 0 y2(0) = ˙xC (0) y3(0) = 0 y4(0) = ˙θ(0)    (12) 第一阶段和第二阶段的分界点时刻 t ∗ 可以由 Matlab 积分时的中断条件获得。在积分之前，先把 积分中具体的要求在 options 中说明，其格式为 options = odeset( ′RelTol′ ,number1, ′AbsTol′ , number2, ′ events′ , @event name); (13) 格式中黑正体表示固定格式，白斜体表示变量。注意 ‘events’，@ 表示在积分时会调用名为 event name 这个子程序，具体的中断条件可以在该子程序中表 达，在本问题中即 A 点速度为零。在设定积分条件 后，再调用积分函数，其格式为 [t, y] = ode45(@rg kt name, [start time : step time : end time], y0, options); (14) 格式中黑正体表示固定格式，白斜体表示变量。其 中 [t, y] 是积分的时间和结果；动力学方程在子程序 rg kt name 中描述；y 是初始条件；如果没有中断条 件，积分时按等步长 step time 从 start time 开始到 end time 结束，有中断条件后，最后一步的积分步 长会自动调整，结束时间也会自动提前。 在程序中先给参数赋值，然后分段计算，这一部 分的源代码见图 4。 图 4 带中断条件的分段积分

756力学与实践2021年第43卷图4中rg-kt-pingpangl，rg-kt-pingpang2分别c(0) =1 m/s, (0)= 0根据方程（9）和方程（10)解出的结果进行积分：在而9o.β=n/(2m），根据情况改变。event-pingpang1中满足A点速度为零就中断积分；取o=-50rad/s.β=0.1,=r/100,第—阶第一阶段结束后把最后时刻的参数赋给tf和new_yo段的位移和速度曲线见图6和图7。数值解（圆点）作为第二阶段的初值，再进行第二阶段积分。第一阶与解析解（位移用圆圈，半径与角度乘积用方框）对段具体中断的条件见图5。比表明，两者很吻合，只有最后一个圆圈或方框与圆3数值计算与解析解的比较点不重合。原因是数值计算中运用了中断，最后一步不是等步长的；而解析解的最后一步的时间没有解在正式进行大量计算前，首先对程序进行一些析表达式，按等步长来处理。测试。由于第一阶段有解析解，可以用于验证数值积图6和图7的结果表明数值计算的精度很好（在分。下面参数在所有计算中保持不变前面文章中已介绍精度~10-10），另一方面，把解m=3g, r=2 cm, g=9.8m/s?析解与数值解结合起来，可以取长补短，更好地处理22mr2,rc(0)=0μ=0.3,Jc=第二阶段的情况。2function [value,isterminal,direction]=event_pingpangl(tt,iy)global radius%%which_mode：0结束：1打滑：2不打滑：xc=iy(1);dxc=iy(2):theta-iy(3);dtheta_dt-iy(4):%A点速度表达式vA-dxc-radius*dtheta_dt：%A点速度判断A点速度36if(vA0不会退出end一退出Svalue-flag+vA;%value=o就退出积分%结束标志isterminal-l:%从什么方向结束direction=o;图5中断子程序1.00.1520000.50.10c积分值。rc理论值(t-s-)/e迪r积分值or理论值drc积分值0.050.5drc理论值-O0.10o积分值田自由家电由rde理论值-0.15-1.000.050.100.150.200.250.3000.050.100.150.200.250.30时间/s时间/s图6位移的数值解与解析解图7速度的数值解与解析解(C)1994-2022ChinaAcademic Journal ElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

756 力 学 与 实 践 2021 年 第 43 卷 图 4 中 rg kt pingpang1，rg kt pingpang2 分别 根据方程 (9) 和方程 (10) 解出的结果进行积分；在 event pingpang1 中满足 A 点速度为零就中断积分； 第一阶段结束后把最后时刻的参数赋给 tf 和 new y0 作为第二阶段的初值，再进行第二阶段积分。第一阶 段具体中断的条件见图 5。 3 数值计算与解析解的比较 在正式进行大量计算前，首先对程序进行一些 测试。由于第一阶段有解析解，可以用于验证数值积 分。下面参数在所有计算中保持不变 m = 3 g, r = 2 cm, g = 9.8 m/s 2 µ = 0.3, JC = 2 3 mr2 , xC (0) = 0 x˙ C (0) = 1 m/s, θ(0) = 0 而 ˙θ0, β = n/(2m), δ 根据情况改变。 取 ˙θ0 = −50 rad/s, β = 0.1, δ = r/100，第一阶 段的位移和速度曲线见图 6 和图 7。数值解 (圆点) 与解析解 (位移用圆圈，半径与角度乘积用方框) 对 比表明，两者很吻合，只有最后一个圆圈或方框与圆 点不重合。原因是数值计算中运用了中断，最后一步 不是等步长的；而解析解的最后一步的时间没有解 析表达式，按等步长来处理。 图 6 和图 7 的结果表明数值计算的精度很好 (在 前面文章中已介绍精度 ∼ 10−10)，另一方面，把解 析解与数值解结合起来，可以取长补短，更好地处理 第二阶段的情况。 图 5 中断子程序 0 0.05 0.10 0.15 0.20 0.25 0.30 ᰦ䰤/s -0.15 -0.10 -0.05 0 0.05 0.10 0.15 ս〫/m xC〟࠶٬ xC⨶䇪٬ r *θ〟࠶٬ r *θ⨶䇪٬ 图 6 位移的数值解与解析解 0 0.05 0.10 0.15 0.20 0.25 0.30 ᰦ䰤/s -1.0 -0.5 0 0.5 1.0 䙏ᓖ m.s-1  dxC〟࠶٬ dxC⨶䇪٬ r *dθ〟࠶٬ r *dθ⨶䇪٬ 图 7 速度的数值解与解析解

第5期757高云峰：Matlab求解理论力学问题系列（四）乒乓球的滚动第二阶段也有解析解，但把数值计算与解析解数值解)。进行对比，发现第一阶段数值解与理论解吻合，第二如果用第一阶段数值计算中断时的最后结果来阶段差别很大（图8和图9)。原因是方程（8）中的参确定第二阶段的解析解中的初始条件，可以发现两数与第一阶段的中止时间有关，由于没有解析解，于阶段的数值解（圆点）全部落在代表解析解的圆圈或是按等步长来计算时间，导致解出的常数Ci和C菱形中了（图10和图11)。这一结论既证明了数值有偏差（表1给出了第一阶段终止时刻的理论解和解的精确性，也表明了解析解与数值解结合的潜力。0.302....·0积分值0.2500理论值100理论值2-0080三0.20eEe0D1015心积分值oooooooo0000000.100理论值1理论值20.051002300.51.01.52.02.53.001时间/s时间/s图8位移的数值解与解析解图9转角的数值解与解析解表1积分常数及其他参数的变化e/radic/(m-s-1)/(rad-s-1)t*/src/mCi/mC2/m解析解0.16150.08575.408213.9450.3004.95904.7975数值解0.15810.16505.70488.24780.27335.61615.45800.30ee积分值20.25理论值10积分值理论值20%理论值1()/季0-理论值24转-2G8扫0.100.050d02310.51.01.52.02.53.00时间/s时间/s图10位移的数值解与修正的解析解图11转角的数值解与修正的解析解有了前面的验证，计算的可靠性有了保证。乒根据式(3)，得到速度表达式为乓球游戏有两种结果：模式1初始角速度小，乒乓1(gμ/β-2z0)e-2Btic-29/8球只会前进不能后退；模式2初始角速度大，乒乓(15)(ur-)gt/r2+009-球先前进再后退。(C)1994-2022ChinaAcademic JournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 5 期 高云峰：Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 757 第二阶段也有解析解，但把数值计算与解析解 进行对比，发现第一阶段数值解与理论解吻合，第二 阶段差别很大 (图 8 和图 9)。原因是方程 (8) 中的参 数与第一阶段的中止时间有关，由于没有解析解，于 是按等步长来计算时间，导致解出的常数 C1 和 C2 有偏差 (表 1 给出了第一阶段终止时刻的理论解和 数值解)。 如果用第一阶段数值计算中断时的最后结果来 确定第二阶段的解析解中的初始条件，可以发现两 阶段的数值解 (圆点) 全部落在代表解析解的圆圈或 菱形中了 (图 10 和图 11)。这一结论既证明了数值 解的精确性，也表明了解析解与数值解结合的潜力。 0 1 2 3 0 0.05 0.10 0.15 0.20 0.25 0.30 ս〫/m ᰦ䰤/s xC〟࠶٬ ⨶䇪٬ ⨶䇪٬ 图 8 位移的数值解与解析解 0 0.5 1.0 1.5 2.0 2.5 3.0 -6 -4 -2 0 2 4 䖜䀂/(°) ᰦ䰤/s θ 〟࠶٬ θ ⨶䇪٬ θ ⨶䇪٬ 图 9 转角的数值解与解析解 表 1 积分常数及其他参数的变化 xC /m θ/rad x˙ C /(m·s−1 ) ˙θ/(rad·s−1 ) t ∗/s C1/m C2/m 解析解 0.161 5 0.085 7 −5.408 2 13.945 0.300 4.959 0 4.797 5 数值解 0.158 1 0.165 0 −5.704 8 8.247 8 0.273 3 5.616 1 5.458 0 0 1 2 3 0 0.05 0.10 0.15 0.20 0.25 0.30 ս〫/m ᰦ䰤/s xC〟࠶٬ ⨶䇪٬ ⨶䇪٬ 图 10 位移的数值解与修正的解析解 0 0 0.5 1.0 1.5 2.0 2.5 3.0 -6 -4 -2 2 4 䖜䀂/(°) ᰦ䰤/s θ 〟࠶٬ θ ⨶䇪٬ θ ⨶䇪٬ 图 11 转角的数值解与修正的解析解 有了前面的验证，计算的可靠性有了保证。 乒 乓球游戏有两种结果：模式 1 初始角速度小，乒乓 球只会前进不能后退；模式 2 初始角速度大，乒乓 球先前进再后退。 根据式 (3)，得到速度表达式为 x˙ C = 1 2 (gµ/β − 2 ˙x0) e−2βt − 1 2 gµ/β ˙θ = 3 2 (µr − δ)gt/r2 + ˙θ0    (15)

758力学与实践2021年第43卷在第一阶段，质心速度在小阻尼时，利用停在出发点前方），分界线右侧（角速度较大）乒乓球e-2Bt～1-2Bt，其斜率k1～2βto-gp；牵连速最终会停在出发点的后方。如果希望乒乓球向前运度re的斜率是kz=(3/2)(μ-8/r)g，在第二阶段，动后返回，正好返回原点，初始角速度是多少（注意质心速度的斜率k3=(3/5)g入s/r。方程（8）没有解析解)？Matlab可以把图形放大，只但是当A点速度为0时（牵连速度与相对速度需要把图14中分界线与曲线的交点区域不断放大，曲线相交），乒乓球就开始作纯滚动，模式1中乒乓就能得到分界线处的速度比值是1.7706653870405，球继续前进，模式2中乒乓球会返回（图12）：而在然后以这样的速度比值单独再算一次，就可以让乒速度图中，牵连速度与相对速度会一直相等且逐渐乓球先前进然后正好退回到原点（图15)。变为0（静止状态)，模式1速度从正趋于0（一直前进），模式2速度从负趋于0（先前进再后退），见21分界线图13。-I00.3E-1模式1：单纯前进0.2o0模式2：先前进再后退-2SAIE0B只能前进可以后退-30的区域的区域吧0.1-420130.2BBRHO速度比0.3图14不同速度比对最终位移的影响12304时间/s0.20图12两种模式下质心位移曲线正好返回原点n0.151.0000ka0.5目0.100(r_s-u)/0.050.5影1.00质心速度100e牵连速度11231.5质心速度2时间/s牵连速度2-2.0图15正好退回原点的情况12340时间/s4小结图13两种模式下速度曲线本文讨论了乒乓球运动的理论解和数值解。在如果关心角速度对最终位移的影响，可以把初退化情况下，可以获得全部解析解，而在一般情况下始速度比（定义为[ro/col，即牵连速度与质心（如同时考虑摩擦、滚动摩阻、空气阻力），只能得到速度之比）作为横坐标，进行大量计算得出曲线如很少的解析解，而采用数值计算的方法可获得一般图14所示，存在一条分界线，分界线左侧（角速度情况下所关心的所有参数。较小）乒乓球最终停在前方（包括先前进再后退，但(1）乒乓球的运动分不同的阶段，本文介绍了在(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.All rights reserved.http://www.cnki.net

758 力 学 与 实 践 2021 年 第 43 卷 在第一阶段， 质心速度在小阻尼时， 利用 e −2βt ≈ 1 − 2βt，其斜率 k1 ≈ 2βx˙ 0 − gµ；牵连速 度 r ˙θ 的斜率是 k2 = (3/2)(µ − δ/r)g，在第二阶段， 质心速度的斜率 k3 = (3/5)gλδ/r。 但是当 A 点速度为 0 时 (牵连速度与相对速度 曲线相交)，乒乓球就开始作纯滚动，模式 1 中乒乓 球继续前进，模式 2 中乒乓球会返回 (图 12)；而在 速度图中，牵连速度与相对速度会一直相等且逐渐 变为 0 (静止状态)，模式 1 速度从正趋于 0 (一直 前进)，模式 2 速度从负趋于 0 (先前进再后退)，见 图 13。 0 1 2 3 4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 ⁑ᔿ1:অ㓟ࡽ䘋 ⁑ᔿ2:ࡽݸ䘋޽ਾ䘰 1:অ㓟 ࡽݸ:2 䘋ࡽ ⁑ᔿ2:ࡽݸ䘋޽ਾ䘰 ᰦ䰤/s 䐍⿫/m 图 12 两种模式下质心位移曲线 0 1 2 1.0 0.5 0 -0.5 -1.0 -1.5 -2.0 3 4 ᰦ䰤/s 䍘ᗳ䙏ᓖ1 牵连速度1 䍘ᗳ䙏ᓖ2 牵连速度2 䙏ᓖ/(m.s-1 ) k k k 图 13 两种模式下速度曲线 如果关心角速度对最终位移的影响，可以把初 始速度比 (定义为 |r ˙θ0/x˙ C0|， 即牵连速度与质心 速度之比) 作为横坐标，进行大量计算得出曲线如 图 14 所示，存在一条分界线，分界线左侧 (角速度 较小) 乒乓球最终停在前方 (包括先前进再后退，但 停在出发点前方)，分界线右侧 (角速度较大) 乒乓球 最终会停在出发点的后方。如果希望乒乓球向前运 动后返回，正好返回原点，初始角速度是多少 (注意 方程 (8) 没有解析解)？Matlab 可以把图形放大，只 需要把图 14 中分界线与曲线的交点区域不断放大， 就能得到分界线处的速度比值是 1.770 665 387 040 5， 然后以这样的速度比值单独再算一次，就可以让乒 乓球先前进然后正好退回到原点 (图 15)。 0 1 2 3 4 ս〫/m 2 1 0 -1 -2 -3 -4 䙏ᓖ∄ 㓯⭼࠶ ਚ㜭ࡽ䘋 Ⲵ४ฏ ਟԕਾ䘰 Ⲵ४ฏ 图 14 不同速度比对最终位移的影响 0 1 2 3 0 0.05 0.10 0.15 0.20 ↓ྭ䘄എ৏⛩ ᰦ䰤/s 䐍⿫/m 图 15 正好退回原点的情况 4 小结 本文讨论了乒乓球运动的理论解和数值解。在 退化情况下，可以获得全部解析解，而在一般情况下 (如同时考虑摩擦、滚动摩阻、空气阻力)，只能得到 很少的解析解，而采用数值计算的方法可获得一般 情况下所关心的所有参数。 (1) 乒乓球的运动分不同的阶段，本文介绍了在

第5期759高云峰：Matlab求解理论力学问题系列（四）乒乓球的滚动Matlab（R2016b版本）中如何处理分段计算的方法，质心最终位移的关系，得到了乒乓球质心最终是否前进的分界线的解，也可以让乒乓球正好返回原处。以及微分代数方程求解的一般过程。这些结论在一般情况下都没有解析解。(2)把解析解与数值解相结合，证明数值精度很高。参考文献(3）适当运用Matlab中的符号推导，可以获得部分解析解表达式（包括部分微分方程的解析解）。1甘勤涛，胡仁喜，程政田等，Matlab数学计算与工程分析从入(4)通过计算，得到了一般情况下乒乓球运动的门到精通.北京：机械工业出版社，2019规律，例如在初始速度一定时，找到了初始角速度和(责任编辑：胡漫)(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 5 期 高云峰：Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 759 Matlab (R2016b 版本) 中如何处理分段计算的方法， 以及微分代数方程求解的一般过程。 (2) 把解析解与数值解相结合，证明数值精度很 高。 (3) 适当运用 Matlab 中的符号推导，可以获得 部分解析解表达式 (包括部分微分方程的解析解)。 (4) 通过计算，得到了一般情况下乒乓球运动的 规律，例如在初始速度一定时，找到了初始角速度和 质心最终位移的关系，得到了乒乓球质心最终是否 前进的分界线的解，也可以让乒乓球正好返回原处。 这些结论在一般情况下都没有解析解。 参 考 文 献 1 甘勤涛, 胡仁喜, 程政田等. Matlab 数学计算与工程分析从入 门到精通. 北京: 机械工业出版社, 2019 (责任编辑: 胡 漫)