求简谐下滑曲线的方程

hujunhua

微信视频号“喵喵喵世界”的一个视频（点击链接）：简谐运动视频
演示了10个小球在一个盘子内各自沿着直径槽道上下滚动，看起来保持圆轮队形沿着盘子的边缘匀速滚动。
所谓盘子，是把小球看作足够小的质点（小球的球 心），无视宏观槽道，当作贴着圆盘内壁沿着直径（盘子内壁的母线）滑动。

熟悉圆内旋轮线者应知其奥妙：直径1:2的圆内旋轮线正是大圆的直径。当小圆沿着大圆匀速滚动时，小圆上的一个定点沿着大圆的直径作简谐振动（图1）。

图3 1:2圆内旋轮线

图3 沿槽道的竖直截面

               图1  1 ：2圆内旋轮线                         图2  沿槽道的竖直截面（盘子内壁的母线）
注：当小圆沿的大圆内侧从(0,-1)滚动到A处时，
      P点从圆心前进到A在 x 轴上的投影处，
      OP = sinφ。
        
如果盘子内壁的母线是一段比较平坦的圆弧（图2），那么小球的运动就等同于一个小角度单摆，在伽利略时代就知道它近似于简谐振动。所以小球的队形看起来就像保持一个旋轮，实际上只是近似于圆形。
到了伯努利时代，已知严格的简谐振动要沿着倒置的旋轮线（又称为等时曲线）。注意这里所说的简谐振动是沿着弧线的，即 F= -kx 中的 x  是母线的弧长。
所以，圆盘内壁的母线如果是旋轮线，即使整体不是小角度，小球的队形也是一个卷曲地粘贴在盘子内壁的“圆形”（某种意义上）。
但是，如果要俯拍视频中小球的队形看起来特别圆，那就要求小球在水平面上的投影是简谐振动，那么问题来了：

问：当圆盘内壁的母线是一个什么曲线时，小球在水平方向的投影是严格的简谐振动？试求该曲线的方程。
不妨设定几个边界条件：
1、-1≤x≤1. 即水平简谐振动的振幅为 1 （盘子的直径是2）。
2、在最高点时水平方向加速度大小等于g/2.

northwolves

假设最高点 $H, A=2H, y=(\sqrt[A^2  (1 + x^2) - x^4] - A)/(A^2 - x^2)$

wayne

以下用 `\dot x, \dot y, \ddot x, \ddot y`都表示对时间 `t`求导的导数, `y', y'' `表示对 `x` 求导。
设滑块质量为 1， 端点高度为`H`，重力加速度为 `g`。
在重力作用下, 根据机械能守恒, 得到 $1/2(dotx^2+dot y^2)+gy=gH$,
又因为$x(t)=\sin(ωt),  dot x=ω cos ωt$ , 代入上式得到微分方程 \[ω^2\cos^2 ωt+\dot y^2+2gy=2gH\tag1
\]由于$t=0,\dot y=0, y=0$, 得到$2gH=w^2$,于是 \[
\dot y^2+2gy=w^2\sin^2wt\tag2
\] 代入$x(t)=sin wt,  \dot x=w\cos wt,\dot y=y' \dot x$ .得到\[
(1-x^2)y'^2+\frac{y}{H}=x^2\tag3
\]
这个好像不容易计算. 只能数值计算了.

这题目可以一般化, 比如$x(t)=a-b t$,, 但就算是这样的表达式,虽然可以得到解析解,但是解析解 仍然非常复杂

282842712474

有点没懂，“在最高点时水平方向加速度大小等于g/2”这个条件是否必要？

282842712474

贴一个Kimi的解（还没仔细审查）：https://www.kimi.com/share/19b5a6c0-c212-8a89-8000-0000d51fa099

Ickiverar

令水平方向运动为简谐运动 $x=A\sin(\omega t)=A\sin(\phi)$, 则竖直方向运动 $y(\phi)$ 满足
$y''(\phi)=B+A^2(\frac{\cos\phi\sin\phi}{y'(\phi)}-1/B), if \sin\phi!=0;$
$y''(\phi)=B, if \sin\phi=0.$
$y(0)=y'(0)=0$
其中 $\omega=\sqrt{\frac{gB}{A^2-B^2}}$，$g$为重力加速度。
注意求导的自变量是$\phi=\omega t$，所以$y'(\phi)=\frac{v_y(t)}{omega}$。

令无量纲量$u=B/A$，无量纲坐标$(\xi,\eta)=(x/A,y/A)$，则有
$\xi=\sin(\phi);$
$\eta''(\phi)=u+(\frac{\cos\phi\sin\phi}{\eta'(\phi)}-1/u), if \sin\phi!=0;$
$\eta''(\phi)=u, if \sin\phi=0.$
其中相位$\phi$可视为为无量纲时间。

Ickiverar

当$u>0.42...$时，解似乎无物理意义，因为在$\phi=\pi/2$时$y'(\phi)$不能抵达$0$。

这个原因可能是满足水平方向简谐运动的轨道在水平方向的两端处斜率不能太陡，否则势能和动能的转换不能满足简谐运动的局域条件？

当$u<0.2...$时，数值求解会比较困难，$y$很容易就掉到负无穷去了，但似乎降低步长就可以解决，$u$似乎可以取到任意小的正值。

Ickiverar

我把方程简化了一下，最终需要求解的方程是$\frac{\text{d}y}{\text{d}t}=\sqrt{\sin^2t-ay}$：

1. Clear[x,equ,Energy,EnergyConst,EnergyPhi];
   
2. x=A Sin[\[Omega] t];
   
3. (*能量，实际上是两倍能量=vx^2+vy^2+2gy*)
   
4. Energy=D[x,t]^2+D[y[t],t]^2+2g y[t]
   
5. (*将能量表达式自变量由t换为\[Phi]=\[Omega]t，相位\[Phi]和t是线性关系，可以看作无量纲时间。*)
   
6. EnergyPhi=Energy/.{y[t]->y[\[Phi]],y'[t]->y'[\[Phi]]\[Omega]}/.t->\[Phi]/\[Omega]
   
7. (*根据能量守恒，上式应为常数。考虑t\[Equal]0时，x\[Equal]0，小球位于水平方向简谐运动平衡点，x方向速度最大，位于y方向最低点y\[Equal]0，且y方向速度为0，即y'[\[Phi]]\[Equal]0。则该常数为*)
   
8. EnergyConst=EnergyPhi/.{y[\[Phi]]->0,y'[\[Phi]]->0}/.\[Phi]->0
   
9. (*现在已经得到微分方程*)
   
10. EnergyPhi-EnergyConst
    
11. (*令\[Eta]为 y/A，a\[Equal](2 g)/(A \[Omega]^2)，可以将上式无量纲化：*)
    
12. equ=FullSimplify[(EnergyPhi-EnergyConst)/(A \[Omega])^2/.{y[\[Phi]]->A \[Eta][\[Phi]],y'[\[Phi]]->A \[Eta]'[\[Phi]],g->a A \[Omega]^2/2}]
    
13. (*考虑 0<\[Phi]<\[Pi]/2的区间，此时小球在上升，取上式中 \[Eta]'[\[Phi]]>0 的解，则最终需解的方程为*)
    
14. DSolve[{\[Eta]'[\[Phi]]==Sqrt[Sin[\[Phi]]^2-a \[Eta][\[Phi]]],\[Eta][0]==0},\[Eta][\[Phi]],\[Phi]]

复制代码

hujunhua

wayne在3#得出的方程是正确的, 应用边界条件能进一步得到 $H=1/4$

设滑行速度为 `v`, 曲线切线倾角为  `φ`, 切线斜率`y'=\tan φ`。
沿用3#的符号， `y',y''`表示对 `x` 求导, $dot x,ddot x, dot y, ddot y$表示对时间 `t` 求导。
则切向加速度为 `\dot v=-g \sin φ`, 水平速度为 $\dot x=vcosφ$.
所以水平加速度为\[
\ddot x=(v\cos φ)'_t=\dot v\cos φ-v\sin φ·\dot φ=-g\sin φ\cos φ-v\sin φ·\dot φ\tag{1}
\]
在端点 `v=0`, 所以水平方向最大加速度为 $-g sin φ cos φ$, 其大小不超 `g/2`。题设端点加速度大小为 $g/2$, 蕴含端点坡度为 45°，原来是极大值. Ickiverar在7#也指出了端点斜率不能大太。
  $x=sin ωt,\dot x=ωcosωt,\ddot x=-ω^2x$
可见水平方向速度最大值平方和加速度最大值均为 $ω^2$, 故有 `ω^2=2gh=g/2`→ `h=1/4, ω=\sqrt{g/2}`
将$H=1/4$代入3#的 (3) 式得
\[
(1-x^2)y'^2+4y=x^2\tag{3}
\]试了一下，Mathematica14都解不出这个非线性方程。
将`y'=\dot y/\dot x`代入可得参数方程\[
\begin{cases}x=\sin ωt\\\dot y^2=ω^2(\sin^2 ωt-4y)\end{cases}\tag{4}
\]取 `g=2→ω=1`, 代入得\[
\begin{cases}x=\sin t\\\dot y^2=\sin^2 t-4y\end{cases}\tag{5}
\]第2个微分方程与楼上是一致的，用Mathematica14也解不出来。

wayne

\[
(1-x^2) y'^2+4y=x^2
\]解得的级数解是

1. FullSimplify[url=home.php?mod=space&uid=6175]@[/url] Normal[AsymptoticDSolveValue[{y'[x]^2 + 1 == (1 - 4 y[x])/(1 - x^2), y[0] == 0}, y[x], {x, 0, 15}]]

复制代码

\[-\frac{5 \left(539096+398507 \sqrt{2}\right) x^{10}}{60082624}-\frac{\left(991+741 \sqrt{2}\right) x^8}{15778}-\frac{1}{392} \left(38+29 \sqrt{2}\right) x^6-\frac{1}{28} \left(5+4 \sqrt{2}\right)
   x^4-\frac{1}{2} \left(1+\sqrt{2}\right) x^2\]

\[
\begin{array}{l}
\{2\}\to \frac{1}{2} \left(\sqrt{2}-1\right) \\
\{4\}\to \frac{1}{28} \left(4 \sqrt{2}-5\right) \\
\{6\}\to \frac{1}{392} \left(29 \sqrt{2}-38\right) \\
\{8\}\to \frac{741 \sqrt{2}-991}{15778} \\
\{10\}\to \frac{5 \left(398507 \sqrt{2}-539096\right)}{60082624} \\
\{12\}\to \frac{720697 \sqrt{2}-982152}{28815136} \\
\{14\}\to \frac{5 \left(1556988441871 \sqrt{2}-2132772637848\right)}{394632047321344} \\
\{16\}\to \frac{3 \left(146153667670419 \sqrt{2}-200965448001404\right)}{27278940271087904} \\
\{18\}\to \frac{5577416513496157005 \sqrt{2}-7691649054995007296}{415512818209210953728} \\
\{20\}\to \frac{382274845963828979775 \sqrt{2}-528417116345422804928}{33448781865841481775104} \\
\{22\}\to \frac{8825496040763685318196437 \sqrt{2}-12222733902541405671775552}{892948680690504197468176384} \\
\{24\}\to \frac{15815453984392449125957646369 \sqrt{2}-21938161381073333601837761344}{1826749763522598961970521837568} \\
\{26\}\to \frac{526864747435226624256851925681610149 \sqrt{2}-731814521130374706391111765711823360}{68728873174872637944010733879574003712} \\
\{28\}\to \frac{2572133743225521495968838539079515036807 \sqrt{2}-3576809893268935045151999051418931038208}{375500198590916657407102644551052569280512} \\
\{30\}\to \frac{1301511075185300333676450611035735994747343571 \sqrt{2}-1811692688679492890043701881732847702383593984}{210974035577909342130865806635719783945117106176} \\
\{32\}\to \frac{11442145613755372645211283277249453402129263652007 \sqrt{2}-15941343912322709126027496733168265258106473208288}{2045439425498113739627335122228504869050648338743296} \\
\{34\}\to \frac{3 \left(423972869450180339192265214202670695402023814502197329 \sqrt{2}-591142647652215400577478832801243270222614081887686656\right)}{249249066633498147856028548654276689323035803965903077376}
   \\
\{36\}\to \frac{3 \left(241497139595400960755646913953543561605064168603825289 \sqrt{2}-336950670131267026730932132073807106065799019401875456\right)}{154810406438453001869055794060355244752542999349064368128}
   \\
\{38\}\to \frac{5 \left(2146090975517365160727709924590836799632418482056181035890600347
   \sqrt{2}-2996206299652137648892611263684715781412476601205754915738005504\right)}{2488437095863444102559947377260233969067374517738407168513427898368} \\
\{40\}\to \frac{155676033619008589548102244675766500616812671859230872185363458155
   \sqrt{2}-217464497704669351857310146858272216863264307836920345226375483392}{39015138753001855750850603522044382586449193331684312392049815977984} \\
\{42\}\to \frac{21699163138967256491129534029822458107820422496591195515953069859115591918041
   \sqrt{2}-30327047276537941519099465298535227467088457922360313459008343183779789553664}{5854629608312031730215436522489073181478384065195519534867591894383111570456576} \\

\end{array}
\]

1. a[1] = 0;
   
2. a[n_] := If[Mod[n, 2] == 1, 0, Subscript[a, n]]; Do[
   
3. fx = Sum[If[i == 2 m, Subscript[a, 2 m], a[i]] x^i, {i, 1, 2 m + 1}];
   
4.   sol = SolveValues[
   
5.    Coefficient[1 - 4 fx - (D[fx, x]^2 + 1) (1 - x^2), x^(2 m)] == 0,
   
6.    Subscript[a, 2 m]];
   
7. If[Length[sol] > 0, a[2 m] = FullSimplify[Last@sol]], {m, 20}];
   
8. fx

复制代码