一個問題解了好久也沒解出來：一豎直棍上連接一質量分布均勻的柔軟有限長繩子，當棍以某角速度旋轉起來時繩子被甩開，求穩定時繩子的形狀。

對數函數曲線顯然是不對的，可答案到底是什麼？

沿繩的自然坐標系

假定繩子總長為 [公式] ，質量線密度為常數 [公式] ，自由端為坐標起點，固定端做勻速圓周運動，繩子張力分布為 [公式] ，則有邊界條件如下： [公式] , [公式]

最重要的，繩子的運動方程可以寫為：[公式]

以及繩子不可伸長的約束： [公式] 。

首先，我們來解決一個大家還在爭論中的問題，穩定狀態時，繩子的形狀是否是一條平面曲線。。

於是我們假定穩定狀態的解可以表達為： [公式] ，[公式] ， [公式] 。若 [公式] 則說明繩子的形狀是一條平面曲線。。

於是將上述試探解帶入運動方程，算一算就可以得到

實部： [公式] ， [公式]

虛部： [公式]

其中虛部的方程很簡單，可以直接積分得到 [公式] ，利用邊界條件 [公式] 可知， [公式] 對任意 [公式] 均成立，於是得到 [公式] ，亦即[公式] ，說明繩子的形狀的確是一條平面曲線。。

╮(╯_╰)╭

既然如此，我們重新整理一下重要的運動方程： [公式] ， [公式] ，及約束 [公式] 。

其中關於 [公式] 的運動方程顯然是可以直接積分的，利用邊界條件 [公式] ，得到： [公式] 。

若我們為了簡化計算，類似地假定 [公式] ，利用約束條件和 [公式] ，則可以簡便地計算出 [公式] ，於是結合另一個邊界條件 [公式] 和 [公式] ， [公式] 和 [公式] 都可以形式地算出來：

[公式] ， [公式]

同時結合另一個運動方程， [公式] ，取微分得到： [公式] ， 以及邊界條件 [公式] ， [公式] 。

所以，這就是一個二階微分方程定解問題。。

結合微分方程，我們可以進一步簡化 [公式] 和 [公式] 的表達式，比如

[公式]

[公式]

當然，這個非線性的二階微分方程，我也不會解。。但是做個數值畫畫圖看看繩子是什麼形狀的，非線性微分方程的數值解，留做習題好了。。

一個數值計算的例子。。算算看，如果其中的參數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()

推薦閱讀：
相关文章