一个问题解了好久也没解出来：一竖直棍上连接一质量分布均匀的柔软有限长绳子，当棍以某角速度旋转起来时绳子被甩开，求稳定时绳子的形状。

对数函数曲线显然是不对的，可答案到底是什么？

沿绳的自然坐标系

假定绳子总长为 [公式] ，质量线密度为常数 [公式] ，自由端为坐标起点，固定端做匀速圆周运动，绳子张力分布为 [公式] ，则有边界条件如下： [公式] , [公式]

最重要的，绳子的运动方程可以写为：[公式]

以及绳子不可伸长的约束： [公式] 。

首先，我们来解决一个大家还在争论中的问题，稳定状态时，绳子的形状是否是一条平面曲线。。

于是我们假定稳定状态的解可以表达为： [公式] ，[公式] ， [公式] 。若 [公式] 则说明绳子的形状是一条平面曲线。。

于是将上述试探解带入运动方程，算一算就可以得到

实部： [公式] ， [公式]

虚部： [公式]

其中虚部的方程很简单，可以直接积分得到 [公式] ，利用边界条件 [公式] 可知， [公式] 对任意 [公式] 均成立，于是得到 [公式] ，亦即[公式] ，说明绳子的形状的确是一条平面曲线。。

╮(╯_╰)╭

既然如此，我们重新整理一下重要的运动方程： [公式] ， [公式] ，及约束 [公式] 。

其中关于 [公式] 的运动方程显然是可以直接积分的，利用边界条件 [公式] ，得到： [公式] 。

若我们为了简化计算，类似地假定 [公式] ，利用约束条件和 [公式] ，则可以简便地计算出 [公式] ，于是结合另一个边界条件 [公式] 和 [公式] ， [公式] 和 [公式] 都可以形式地算出来：

[公式] ， [公式]

同时结合另一个运动方程， [公式] ，取微分得到： [公式] ， 以及边界条件 [公式] ， [公式] 。

所以，这就是一个二阶微分方程定解问题。。

结合微分方程，我们可以进一步简化 [公式] 和 [公式] 的表达式，比如

[公式]

[公式]

当然，这个非线性的二阶微分方程，我也不会解。。但是做个数值画画图看看绳子是什么形状的，非线性微分方程的数值解，留做习题好了。。

一个数值计算的例子。。算算看，如果其中的参数k取得很大会怎么样？

╮(╯_╰)╭

不过我必须要提醒一下大家，问题到了这里，已经开始变得有意思了。。因为这是二阶微分方程，通常会具有非平凡解，会出现震荡或者共振行为，也就是接下来我所要说的事情。。

换句话说，其他答主，其实仅仅只计算了这个体系中，最简单的一个平凡解。。

╮(╯_╰)╭

首先，我们来个最简单的近似，假定 [公式] ，这样我们可以把这个非线性的微分方程线性化，得到 [公式] 。这个线性化的方程有解，可以用一阶贝塞尔函数表达出来，结合边界条件，可以定解：

[公式]

换个简明一些的符号， 记 [公式] ，则有 [公式] 。。

那么一阶近似下，就有 [公式] ， [公式] 。。

可以看到，绳子的横向偏移量 [公式] ，正比于 [公式] ，也就是振幅。。而事实上，在后面的时候我会重新提到 [公式] ，它事实上具有比振幅更多的意义。。

简单画一下这个贝塞尔函数，大家就可以看出来，绳子的大致形状了。。

平凡解的大致样子，参数取x0=0.1，L=1

╮(╯_╰)╭

对于熟悉贝塞尔函数的童鞋，看到绳子的横向偏移量 [公式] 跟贝塞尔函数有关，到这里就已经可以看出来了，一旦参数合适，即系统的转动角速度 [公式] 足够大的时候，也就是 [公式] 足够大的时候，那么 [公式] 很可能存在零点。。而零点的存在则意味著，绳子将会穿过转动轴，而不是简简单单地一直向外延伸，这便是该问题中的非平凡解。。

而每当 [公式] ，其中 [公式] 表示 [公式] 的第 [公式] 个零点时，系统发生共振。。而共振时则说明著，一个非平凡解向著另一个非平凡解的过渡。。

由于 [公式] 表示振幅，为了表达共振解则只要取 [公式] 就好了，此条件意味著边界条件 [公式] ，在每个 [公式] 的零点 [公式] 处存在非平凡解， [公式] ， [公式] ，以及 [公式] 。。

简单画一下这个贝塞尔函数，大家就可以看出来，非平凡解的形状了。。可以看出，系统每跨过一个 [公式] 的零点，那么绳子跟转轴就多相交一次。。如果我们把绳子穿过转轴的次数，称作解的阶，那么平凡解也就是0阶解，只出现低角速度情形。。当然，根据初始条件的不同，低阶解是可以出现在高角速度情形的，但是高阶解是不可能出现在低角速度情形的， [公式] 具有其特有的阈值，也就是 [公式] 。。

前4个非平凡解的大致样子，参数取x0-amp;amp;amp;>0，L=1

于是我们可以对这个系统简单地画个相图：

线性化系统的相图，不同颜色区用来分解的阶数，从平凡解（阶=0）开始，向右递增。。

这个相图中，可能出现的构型数目，随著 [公式] 的不断增加，而向右无限增加，因为 [公式] 具有无穷多个零点。。

╮(╯_╰)╭

至此，做了近似之后的线性化的方程已经研究得差不多了。。我们还是回头看一看原先那么非线性的方程会给出什么样的相图好了。。╮(╯_╰)╭

所以我们回过头看看，当初这个线性化的近似，扔掉了什么。。

假定 [公式] ，这样我们将复杂的 [公式] 项近似为了 [公式] ，而原本 [公式] 这一项的是来自约束 [公式] ，也就是绳子不可伸长。。

正是这个约束，导致了方程的非线性效应，结果就是绳子的振幅并不是独立参数，而是转动角速度的函数。。

换句话说，随著 [公式] 的振幅（正比于 [公式] ）不断增加，绳子本身会在竖直方向上缩短，这样对应于系统的转动角速度也要增加，结果就是相变的临界点也要不断向右偏移，高阶解的存在，需要依靠更高的转动角速度才能实现。。

另一方面，绳子不可伸长，意味著边界条件 [公式] 将限制绳子本身可以穿过转轴的次数，进而解的阶数将会被限制。。也就是说，高阶解不仅有角速度的阈值，也有振幅的阈值。。

所以，我们就可以粗略的描述出非线性效应，也就是绳子不可伸长，对线性系统的影响。。于是一个更加贴近实际情况的相图，就可以画出来了。。

考虑非线性效应之后的相图，不同颜色区用来分解的阶数，从平凡解（阶=0）开始，向右下递增。。

比如， [公式] 较大的时候，只存在平凡解。。

当然，这个相图只是我随手画的，具体这个相边界，满足什么样的方程 [公式] ，留作习题好了。。

╮(╯_╰)╭

前几天在李永乐老师那儿看到了这个题，作为李永乐的前学生，就动手算了一下……不过李永乐老师想要的解析解还是没给出来。

首先，稳定情况下，绳子一定在过棍子轴线的竖直平面内。看到仍有人受生活中彩带形状的直觉影响，以为是螺旋线，所以先明确这一点。证明如下：换到与棍子共同旋转的参考系，考虑绳上任一点M到末端的一段绳子，稳定状况下该绳在旋转参考系中静止，无科氏力、无空气阻力且受力平衡。这段绳子受重力、离心力（的积分）和M 点张力，于是张力和另外两力共面。又因为是软绳，绳子不产生侧向力，于是绳子任一点切向与张力平行，故绳子整体位于竖直平面内。

相关参量：线密度 [公式] 、角速度 [公式] 、绳长 [公式] 、重力加速度 [公式] 、棍的半径 [公式] 。简单的量纲分析可知， 绳子形状与 [公式] 无关，容易构造两个无量纲量 [公式] ， [公式] 。若调整参量使得 [公式] 不变，比如等比放大长度量，则绳子的形状不会变化。

接下来列曲线的微分方程，同样用上面的受力分析方法， @Super Mario 的方法会得到含两个积分的微分方程，很难消掉。不如用弧长做参数写参数方程：设曲线形状为 [公式] ,切线与受力同向： [公式] ,结合弧长做参数的条件 [公式] ,初值条件 [公式] ,可以做 Mathematica 数值解：

不妨取 [公式] ，让 [公式] 等间距取值，这是 [公式] 时的情况：

拿 Manipulate 写个变动轴半径的演示：

变动轴半径 R 的情况

本来还做了好多演示，不过……检查的时候发现这个数值解带回原方程积分似乎有点偏差，而且和 @CDK6182CHR 画的图不大一样，不知是我方程写错了还是数值解解错了。之前看这题有一位答主也是这个思路解的方程，但不知为什么删了，或许也是发现了哪里不对。最近时间略紧，有空再换个思路做下计算，比如等周约束下做变分什么的。

附上MMA代码：

Module[{g = 1, L = 1, [CapitalPi]1 = 0.1, sol},
sol = ParametricNDSolveValue[{y[l]/
x[l] == -((g (L - l))/([Omega]^2 !(
*SubsuperscriptBox[([Integral]), (l), (L)](x[
l] [DifferentialD]l)))), x[l]^2 + y[l]^2 == 1,
x[0] == [CapitalPi]1 L, y[0] == 0}, {x[l], y[l]}, {l, 0,
L}, [Omega]];
ParametricPlot[
Evaluate@Table[sol[[Omega]], {[Omega], 0, 5, .5}], {l, 0, L},
PlotRange -&> {{0, (1 + [CapitalPi]1) L}, {-1 L, 0}},
PlotLegends -&> Range[0, 5, .5]]]

我们假设，不计空气阻力，转动的加速度极为缓慢，使得绳子保持在同一平面内而不是螺旋线。棍子与绳子所在的平面中，棍子的半径是1，绳子的长度是l，单位长度质量为 [公式] ，连接点的高度为0，然后建立如下的平面直角坐标系：

然后我们对(x,x+dx)点作受力分析。这一点绳子的长度为 [公式] ，乘以 [公式] 得到绳子的质量。这个点的绳子一共受到四个力，左上方的张力 [公式] ，右下方的张力 [公式] ，重力 [公式] ，三个力的合力提供了绳子的向心力 [公式] 。注意这里 x点处张力与平面的夹角为 [公式] ，x+dx处张力与平面的夹角为 [公式] ，我们可以列出如下的两个方程

[公式]

[公式]

这个时候，我们令 [公式] ，即可得

[公式]

[公式]

通过积分，即可得

[公式]

[公式]

我们假设，绳子的末端位置为t, t&>1，这个时候张力为0，所以积分也可以变换为

[公式]

[公式]

我们知道，因为斜率是负数， [公式] ，所以微分方程可以转化为

[公式]

要想得到解析解，就需要解这个微分方程，找到一个合适的 [公式] 来满足上面的关系，但是这个很难找到。

不过如果求数值解，还是可行的。具体的思路，就是把f(x)看成另一个函数，列出对应的方程，mathematica软体求解（matlab根本算不了）。或者用迭代法，令 [公式] 开始，从t开始向前步步迭代，积分化为求和，解出f(x)每个点对应的导数，最后画出f(x)图像（这里可以强行让t为1），求出长度为多少t，再利用长度的限制，反向解出t的值。

－－－－－7.1 更新－－－－－－－－

（这里特别感谢 @CDK6182CHR 同学，我之前一直忙著课题，只是提了个数值分析的思路，没有亲自实践，而这位同学用代码实现了，思路跟我上面提到的基本一致，不得不说，数值演算法在解微分方程方面太万能了，大家可以看一下）

@qfzklm 的回答

7.6更新: 建议大家不要关注我了，我平时很水的，基本上没什么有价值的回答……

以下为原回答。

看到了@Super Mario ?的回答，这学期刚刚学过一点粗浅的数值分析，想尝试利用数值方法求解这个方程。几乎是第一次利用数值方法解决问题，如果不妥之处，望不吝赐教。

@Super Mario 给出的方程是：

[公式] （1）

方便起见，直接令 [公式] ，取步长h=0.01, 在0至1间计算。

对于固定的 [公式] ，（1）式可以看做关于 [公式] 的不动点迭代方程。可以对一系列的 [公式] 分别解出以上方程，得到 [公式] 的函数表；再根据 [公式] 可求解出待求的 [公式] 。

下面考虑 [公式] 的求解过程。对于积分[公式] 可用复化梯形公式计算

[公式]

上式中 [公式] 。上式可以递推化地表示为：

[公式]

现在再来考虑（1）式。我们从1开始，每次减去一个步长倒著计算 [公式] 的函数值。也就是说，在计算 [公式] 时， [公式] 都已经全部计算完毕，相应的积分值也是已知量。继而将（1）式离散化，得到关于 [公式] 的不动点迭代方程

[公式] （2）

上式中 [公式] 应视为整体，是未知量； [公式] 和后面的积分值也应该视为整体，是已知量。

现在只剩下最后一个问题，即 [公式] 的初值问题。考虑极限

[公式]

这就是所要求的初值。有了初值，从1开始，每次减去一个步长，用（2）式计算导数值得到倒数的函数值表，再用数值积分方法就可以得到要求的 [公式] 函数值表。我用Python实现了上述思路，作出的图像为

附上Python源代码：

from math import sqrt
import matplotlib.pyplot as plt
import numpy as np

int1 = 0 # 递推复化梯形公式的初值
int2 = 0
h = 0.01 # 步长
t = 1 # x的范围为(0,t)
ys = [-1,] # 从前到后是从1开始倒著的函数值。t=1时初值为1. 前面的函数值往后插入。

def yk(xk:float):
"""
给出yk=y(xk)的函数值。y=f(x)。
"""
global h,ys
k = round((t-xk)/h,0)
return ys[k]

def int1_next(xk):
"""
递推一次积分1的值。xk是要求解的积分下限，保证此时全局变数int1保存的是xk+1开始的积分值；
保证ys[-1],y[-2]有效且是y(x_k),y(x_k+1)的值。
调用此函数时，x_k已经解算出来。
"""
global int1,h,ys
int1+=h/2*(sqrt(1+ys[-1]**2)+sqrt(1+ys[-2]**2))

def int2_next(xk):
global int2,h,ys
int2+=h/2*(xk*sqrt(1+ys[-1]**2)+(xk+h)*(1+ys[-2]**2))

def fixed_fac(xk):
"""
不动点迭代法的工厂函数。
"""
def phi(yk):
global h,int1,int2,ys
yk1 = ys[-1] # y_{k+1}
xk1 = xk + h
return -(h/2*(sqrt(1+yk**2)+(sqrt(1+yk1**2)))+int1)/
(h/2*(xk*sqrt(1+yk**2)+xk1*sqrt(1+yk1**2))+int2)
return phi

def fixed_point(phi)-&>float:
"""
不动点迭代法求解方程x=phi(x)。以0为初值，返回结果。
"""
x=0
x1=phi(x)
while abs(x1-x)&>1e-8:
x=x1
x1=phi(x)
return x

if __name__ == __main__:
x=1
while x&>0:
x-=h
phi=fixed_fac(x)
y=fixed_point(phi)
ys.append(y)
int1_next(x)
int2_next(x)
ys.reverse()
fs = [0,] # y=f(x)，最终数值解
x=0
i=0
while x&<1: x+=h i+=1 # 每一次取的导数值是ys[i] 和 ys[i-1] f1=h/2*(ys[i]+ys[i-1])+fs[-1] fs.append(f1) plt.plot(np.arange(0,1.01,0.01),fs) plt.xlabel(x) plt.ylabel(f(x)) plt.show()

推荐阅读：
相关文章