建立一个吊链模型
找出悬挂在两点间长度为 L 的链条或缆绳在势能最小时的位置.
![](assets.zh/model-a-hanging-chain/O_41.png)
设置链条长度 、 左端高度
和右端高度
的值.
In[1]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_73.png)
L = 4; a = 1; b = 3;
令 为链条高度,它是水平位置的函数,其中
.
In[2]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_74.png)
xf = 1; nh = 201; h := xf/nh;
为链条高度 设置变量.
In[3]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_75.png)
varsy = Array[y, nh + 1, {0, nh}];
用 表示在位置
处的斜率,并为之设置变量.
In[4]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_76.png)
varsm = Array[m, nh + 1, {0, nh}];
用 表示从
到
的势能部分.
In[5]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_77.png)
varsv = Array[v, nh + 1, {0, nh}];
用 表示链条在位置
的长度,并为之设置变量.
In[6]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_78.png)
varss = Array[s, nh + 1, {0, nh}];
将所有变量结合在一起.
In[7]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_79.png)
vars = Join[varsm, varsy, varsv, varss];
我们的目的是将势能 最小化.
In[8]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_80.png)
objfn = v[nh];
下面是由几何形状决定的边界值约束.
In[9]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_81.png)
bndcons = {y[0] == a, y[nh] == b, v[0] == 0, s[0] == 0, s[nh] == L};
将常微分方程 ,
,
离散化.
In[10]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_82.png)
odecons = {Table[
y[j + 1] == y[j] + 0.5*h*(m[j] + m[j + 1]), {j, 0, nh - 1}],
Table[v[j + 1] ==
v[j] + 0.5*
h*(y[j]*Sqrt[1 + m[j]^2] + y[j + 1]*Sqrt[1 + m[j + 1]^2]), {j,
0, nh - 1}],
Table[s[j + 1] ==
s[j] + 0.5*h*(Sqrt[1 + m[j]^2] + Sqrt[1 + m[j + 1]^2]), {j, 0,
nh - 1}]};
选取变量的初始点.
In[11]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_83.png)
tmin = If[b > a, 0.25 , 0.75]; init =
Join[Table[4*Abs[b - a]*((k/nh) - tmin), {k, 0, nh}],
Table[4*Abs[b - a]*(k/nh)*(0.5*(k/nh) - tmin) + a, {k, 0, nh}],
Table[(4*Abs[b - a]*(k/nh)*(0.5*(k/nh) - tmin) + a)*4*
Abs[b - a]*((k/nh) - tmin), {k, 0, nh}],
Table[4*Abs[b - a]*((k/nh) - tmin), {k, 0, nh}]];
在满足约束的条件下将势能最小化.
In[12]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_84.png)
sol = FindMinimum[{objfn, Join[bndcons, odecons]},
Thread[{vars, init}]];
提取解.
In[13]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_85.png)
solpts = Table[{i h, y[i] /. sol[[2]]}, {i, 0, nh}];
绘制势能最小时链条位置的图形.
In[14]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_86.png)
ListPlot[solpts, ImageSize -> Medium, PlotTheme -> "Marketing"]
Out[14]=
![](assets.zh/model-a-hanging-chain/O_42.png)
用 FindFit 将结果与悬链曲线拟合.
In[15]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_87.png)
catenary[t_] = c1 + (1/c2) Cosh[c2 (t - c3)];
In[16]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_88.png)
fitsol = FindFit[solpts, catenary[t], {c1, c2, c3}, {t}]
Out[16]=
![](assets.zh/model-a-hanging-chain/O_43.png)
同时绘制解及悬链曲线.
In[17]:=
![Click for copyable input](assets.zh/model-a-hanging-chain/In_89.png)
Show[Plot[catenary[t] /. fitsol, {t, 0, 1},
PlotStyle -> Directive[Green, Thickness[0.01]],
ImageSize -> Medium],
ListPlot[Take[solpts, 1 ;; nh ;; 5], PlotStyle -> PointSize[.02]]]
Out[17]=
![](assets.zh/model-a-hanging-chain/O_44.png)