2
$\begingroup$

Background

Good evening, I have a big problem with num. solution NDSolve of differential equation. To start with, the model:

There is a fast-moving rope in closed path. Where

$T$ is tension,

$a$ is angle beetwen rope and horizontal plane,

$s$ is curvilinear coordinates $[0,1]$ ($0$ - the beginning of the rope, $1$ - the ending of the rope).

$Dr$ drag coefficient

$W$ weight coefficient

And it gives the differential equations for $T$, $a$ and $s$.

$$\frac{d}{ds}(T(s)\sin\alpha(s))=W+Dr\sin\alpha(s) $$

$$\frac{d}{ds}(T(s)\cos\alpha(s))=Dr\cos\alpha(s) $$

D[T[s] Sin[a[s]], s] == W + Dr Sin[a[s]], 
D[T[s] Cos[a[s]], s] == Dr Cos[a[s]],

Moreover, boundary conditions comes from equation for beginning and ending of the rope. $x[0] = x[1] = 0, y[0] = y[1] = 0$, where $x[s]$, $y[s]$ are coordinates of point on the rope with curve distance to beginning $s$. It means, that ending and beginning of the rope are in the same place. Differential equations for coodinates $x[s],y[s]$ are quite easy.

$$\frac{dy(s)}{ds}=\frac{\tan \alpha(s)}{\sqrt{1+\tan^2 \alpha(s)}} $$

$$\frac{dx(s)}{ds}=\frac{1}{\sqrt{1+\tan^2 \alpha(s)}} $$

(Tan[a[s]])/Sqrt[1 + Tan[a[s]]^2] == D[y[s], s],
1/Sqrt[1 + Tan[a[s]]^2] == D[x[s], s],
x[0] == 0,
y[0] == 0,
x[1] == ϵ,
y[1] == ϵ

I've solved it, but solution is unreal and contradicts boundary conditions. But Mathematica's solution in ParametricPlot looks like this:

**Fig.1** Obtained soltion

Fig.1 Obtained soltion

The rope should be closed, but it's not. And it should look like that:

**Fig.2** Shape of rope in dependence of Dr/W

Fig.2 Shape of rope in dependence of $\frac{Dr}{W}$

Please, help. The final code:

x[s] =.
y[s] =.
NumSol = Block[{\[Epsilon] = $MachineEpsilon},
  With[{Dr = 9.9, W = 8},
   NDSolve[{
     D[T[s] Sin[a[s]], s] == W + Dr Sin[a[s]],
     D[T[s] Cos[a[s]], s] == Dr Cos[a[s]],
     (Tan[a[s]])/Sqrt[1 + Tan[a[s]]^2] == D[y[s], s],
     1/Sqrt[1 + Tan[a[s]]^2] == D[x[s], s],
     x[\[Epsilon]] == 0,
     y[\[Epsilon]] == 0,
     x[1 - \[Epsilon]] == \[Epsilon],
     y[1 - \[Epsilon]] == \[Epsilon]
     },
    {T, a, x, y}, {s, \[Epsilon], 1 - \[Epsilon]},
    Method -> {"StiffnessSwitching", "NonstiffTest" -> False}]]]

ParametricPlot[{x[s], y[s]} /. NumSol // Evaluate, {s, 0, 1},
 PlotRange -> Automatic,
 AspectRatio -> 1,
 AxesLabel -> {"x", "y"}
 ]
Improve this question
$\endgroup$
9
  • $\begingroup$ Where did you get this closed loop picture? $\endgroup$
    – Alex Trounev
    Dec 12, 2019 at 20:54
  • $\begingroup$ "$a$ is angle beetwen rope and gorizontal plane. " So, the range of $a$ is not $(-\frac\pi 2,\frac\pi 2)$? $\endgroup$
    – xzczd
    Dec 13, 2019 at 5:14
  • $\begingroup$ @AlexTrounev I get this picture from this article: [ Propelled Strings: Rising from Friction](doi.org/10.1103/PhysRevLett.123.144501). Author's of the aricle solved the same problem and got this picture, I'm trying to get it too. $\endgroup$
    – danyesss
    Dec 13, 2019 at 15:39
  • $\begingroup$ @xzczd Yes, $a$ must be in this range. $\endgroup$
    – danyesss
    Dec 13, 2019 at 15:44
  • 1
    $\begingroup$ @danyesss The article discusses the analytical solution for $T, \alpha$. Help you find this solution? $\endgroup$
    – Alex Trounev
    Dec 13, 2019 at 18:31

1 Answer 1

Reset to default
3
$\begingroup$

Answer has been significantly revised

Begin by obtaining symbolic solutions for T[s] and a[s].

sat = DSolveValue[{D[T[s] Sin[a[s]], s] == W + Dr Sin[a[s]], 
    D[T[s] Cos[a[s]], s] == Dr Cos[a[s]]}, {a[s], T[s]}, s]

but the result is a bit lengthy to reproduce here. However, simpler expressions can be extracted from sat for a[s]

eqa = Simplify[sat[[1, 1]] /. C[1] -> W*C[1]] == 
    Simplify[sat[[1, 0, 1]][a[s]] /. Dr -> r*W]
(* s/C[1] + C[2] == ((Cos[a[s]/2] - Sin[a[s]/2])^(-1 - r) 
   (Cos[a[s]/2] + Sin[a[s]/2])^(-1 + r) (r - Sin[a[s]]))/(-1 + r^2) *)

where r = Dr/W has been introduced for compactness. T[s]also can be obtained in terms of a[s], although it is not needed for the computation below.

eqT = FullSimplify[sat[[2]] /. {sat[[1]] -> a[s], Dr -> r*W}]
(* C[1] Sec[a[s]] (Cos[a[s]/2] - Sin[a[s]/2])^(-r) (Cos[a[s]/2] + Sin[a[s]/2])^r *)

Analyzing eqa, we see that it is valid for RealC[_] only when - Pi/2 < a[s] < Pi/2. But, based on the reference provided by the OP in a comment above, - Pi/2 < a[s] < - 3 Pi/2 also is needed to solve the problem posed in the question. The latter solution can be obtained by the substitution a[s] -> Pi - a[s]and renormalization of C[1]. Combine the two.

eqaext = eqa[[1]] == Piecewise[{{eqa[[2]], a[s] > -Pi/2}, 
    {-Simplify[eqa[[2]] /. a[s] -> Pi - a[s]], a[s] < -Pi/2}}, 0]
(* s/C[1] + C[2] == Piecewise[{{((Cos[a[s]/2] - Sin[a[s]/2])^(-1 - r)*
   (Cos[a[s]/2] + Sin[a[s]/2])^(-1 + r)*(r - Sin[a[s]]))/(-1 + r^2), a[s] > -Pi/2}, 
   {-(((-Cos[a[s]/2] + Sin[a[s]/2])^(-1 - r)*(Cos[a[s]/2] + Sin[a[s]/2])^(-1 + r)*
   (r - Sin[a[s]]))/(-1 + r^2)), a[s] < -Pi/2}}, 0] *)

Note that the constants, C[_] may not be the same above and below a[s] = - Pi/2. In fact, determining C[_] is the essence of the computation. Here is a plot of eqaext[[2]] for r = 9.9/8, the value employed in the question.

ParametricPlot[{Last[eqaext /. {a[s] -> b, r -> 9.9/8}], b}, 
    {b, -3 Pi/2 + .01, Pi/2 - 0.01}, AxesLabel -> {"s/c1+c2", "a[s]"}, 
    ImageSize -> Large, LabelStyle -> {15, Black, Bold}, AspectRatio -> 1]

enter image description here

A numerical expression fora[s] as a function of Last[eqaext /. r -> 9.9/8} can be obtain readily by

int = Interpolation[Table[{Re[Last[eqaext /. {a[s] -> b, r -> 9.9/8}]], b}, 
   {b, -3 Pi/2 + .0001, Pi/2 - .0001, .0001}]];

Surprizingly, this simple result is rather more robust than InverseFunction for the integrations that follow.

The specific problem that the code in the question attempts to address is equivalent to the respective C[_] being the same throughout - 3 Pi/2 < a[s] < Pi/2. To determine these two constants, require that {x[s], y[s]} both are equal to zero at s = 0 and s = 1, in other words that fi[c1, c2] = {0,0}, where

fi[c1_, c2_] := {NIntegrate[Cos[int[s/c1 + c2]], {s, 0, 1}], 
    NIntegrate[Sin[int[s/c1 + c2]], {s, 0, 1}]}

which is solved by

param = FindRoot[Quiet@fi[c10, c20], {{c10, -.7}, {c20, .6}}, 
    Evaluated -> False] // Values
Quiet[fi @@ %]
(* {-0.0909828, 5.49556} *)
(* {-1.13858*10^-16, 6.41848*10^-17} *)

Finally, a plot of x[t] and y[t] is obtained from

ps = ParametricNDSolveValue[{x'[s] == Cos[int[s/c1 + c2]], 
    y'[s] == Sin[int[s/c1 + c2]], x[0] == 0, y[0] == 0}, 
    {x[s], y[s]}, {s, 0, 1}, {c1, c2}];
ps @@ param;
ParametricPlot[%, {s, 0, 1}, AxesLabel -> {"x", "y"}, 
   ImageSize -> Large, LabelStyle -> {15, Black, Bold}]

enter image description here

For completeness, a[0] is given by

int[param // Last]
(* 0.940888 *)

Let us turn now to reproducing a curve like those in the second figure in the question, for which a[0] is specified to be zero. Then, C[_] definitely are not the same above and below a[s] = - Pi/2. The four seemingly undetermined constants are reduced to two as follow. At s = 0.

c2p = Last[eqaext /. {a[s] -> 0, r -> 9.9/8}]
(* 2.32873 *)

Next, note that s must be continuous at a[s] = - Pi/2, which in turn requires that c1p*c2p = c1m*c2m. (Constants with names ending in p are for - Pi/2 < a[s], and with m for - Pi/2 > a[s].) This reduces the number of free constants to two, as desired. As before, determine them by requiring that {x[s], y[s]} both are equal to zero at s = 0 and s = 1.

f0[c1p_, c1m_] := {NIntegrate[Piecewise[{{Cos[int[s/c1p + c2p]], s < -c1p*c2p}, 
    {Cos[int[s/c1m + c2p*c1p/c1m]], s > -c1p*c2p}}, 0], {s, 0, 1}], 
    NIntegrate[Piecewise[{{Sin[int[s/c1p + c2p]], s < -c1p*c2p}, 
    {Sin[int[s/c1m + c2p*c1p/c1m]], s > -c1p*c2p}}, 0], {s, 0, 1}]}

param = FindRoot[Quiet@f0[c1p0, c1m0], {{c1p0, -.1}, {c1m0, -.01}}, 
    Evaluated -> False] // Values
Quiet[f0 @@ %]
(* {-0.21471, -0.0133781} *)
(* {4.17224*10^-17, -9.19403*10^-17} *)

Finally, a plot of x[t] and y[t] is obtained from

ps0 = ParametricNDSolveValue[{x'[s] == Piecewise[{{Cos[int[s/c1p + c2p]], 
    s < -c1p*c2p}, {Cos[int[s/c1m + c2p*c1p/c1m]], s > -c1p*c2p}}, 0], 
    y'[s] == Piecewise[{{Sin[int[s/c1p + c2p]], 
    s < -c1p*c2p}, {Sin[int[s/c1m + c2p*c1p/c1m]], s > -c1p*c2p}}, 0], 
    x[0] == 0, y[0] == 0}, {x[s], y[s]}, {s, 0, 1}, {c1p, c1m}];
ps0 @@ param;
ParametricPlot[%, {s, 0, 1}, AxesLabel -> {"x", "y"}, 
    ImageSize -> Large, LabelStyle -> {15, Black, Bold}]

enter image description here

which fits well between the r = 1 and r = 1.5 curves in the second figure of the question. Generating all the curves in the second figure would be straightforward.

Improve this answer
$\endgroup$
4
  • $\begingroup$ This is a good algorithm (+1). But how to reproduce fig. 4 from the article journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.144501 ? $\endgroup$
    – Alex Trounev
    Dec 18, 2019 at 11:50
  • $\begingroup$ See Fig. 2 above with $Dr/W=0.2,0.5, 1, 1.5, 2.5, 5, 10$ and with $\alpha(0)=0$. And give me email for copy. $\endgroup$
    – Alex Trounev
    Dec 18, 2019 at 16:10
  • $\begingroup$ I sent you a copy. $\endgroup$
    – Alex Trounev
    Dec 18, 2019 at 19:35
  • $\begingroup$ @AlexTrounev Thanks again for sharing the article. Having read it, I now see how to proceed to obtain the second plot in the question, which I shall do tomorrow. $\endgroup$
    – bbgodfrey
    Dec 19, 2019 at 5:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.