Statics (11/15/2019)

@Karan Chatrath has given a nice summary of the Newton Raphson technique. This method converges very powerfully, but the disadvantage is that it requires a great deal of symbolic manipulation. An alternative approach is a hill-climbing (or valley-descent) algorithm. This algorithm requires a great deal of processing, but it involves only elementary math operations. The broad details are:

1) The unknown quantities are $(x_3,y_3, \theta)$ , where $\theta$ is the angle of the rod with respect to the positive $x$ axis
2) Randomly mutate the $(x_3,y_3, \theta)$ coordinates in order to minimize a "residual"
3) The residual is formed as $|F_x| + |F_y| + |\tau|$ where ( $F_x, F_y, \tau)$ are the net $x$ force, net $y$ force, and net torque
4) After many iterations and mutations, the residual is brought nearly to zero, yielding a stasis point

I have attached graphs of three different stasis configurations. The first one is the "natural" configuration that one would find in a laboratory if one physically did the experiment. This is a stable equilibrium point. The other two are likely to be unstable equilibria, which would be disrupted and undone by any small perturbation. These would consequently not be easily observable in a laboratory. Coordinates for these configurations are included within the attached code file.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

import math
import random

#################################

L = 3.0
M = 1.0

L1 = 1.0
k1 = 5.0

L2 = 2.0
k2 = 3.0

x1 = 0.0
y1 = 0.0

x2 = 3.0
y2 = 5.0

g = 10.0

#################################

delta = 0.001

minres = 99999999.0

x3 = -5.0 + 10.0 * random.random()
y3 = -5.0 + 10.0 * random.random()
theta = 2.0*math.pi * random.random()

for j in range(0,3*10**6):

    if j > 2*10**6:
        delta = 0.0001

    x3m = x3 + delta * (-1.0 + 2.0 * random.random())
    y3m = y3 + delta * (-1.0 + 2.0 * random.random())
    thetam = theta + delta * (-1.0 + 2.0 * random.random())

    x4 = x3m + L*math.cos(thetam)
    y4 = y3m + L*math.sin(thetam)

    v1x = x1 - x3m
    v1y = y1 - y3m

    v2x = x2 - x4
    v2y = y2 - y4

    v1 = math.hypot(v1x,v1y)
    v2 = math.hypot(v2x,v2y)

    u1x = v1x/v1
    u1y = v1y/v1

    u2x = v2x/v2
    u2y = v2y/v2

    s1 = v1 - L1
    s2 = v2 - L2

    Fs1 = k1*s1
    Fs2 = k2*s2

    Fs1x = Fs1 * u1x
    Fs1y = Fs1 * u1y

    Fs2x = Fs2 * u2x
    Fs2y = Fs2 * u2y

    Fgx = 0.0
    Fgy = -M*g

    Fx = Fgx + Fs1x + Fs2x
    Fy = Fgy + Fs1y + Fs2y

    x5 = (x3m+x4)/2.0
    y5 = (y3m+y4)/2.0

    rgx = x5 - x3m
    rgy = y5 - y3m

    rs2x = x4 - x3m
    rs2y = y4 - y3m

    Tg = rgx*Fgy - rgy*Fgx
    Ts2 = rs2x*Fs2y - rs2y*Fs2x

    T = Tg + Ts2

    res = math.fabs(Fx) + math.fabs(Fy) + math.fabs(T)

    if res < minres:

        minres = res

        x3 = x3m
        y3 = y3m
        theta = thetam

        x4_store = x4
        y4_store = y4

        x5_store = x5
        y5_store = y5

        s1_store = s1
        s2_store = s2

#################################

if minres < 0.001:

    print "minres"
    print minres
    print ""
    print "x3,y3"
    print x3
    print y3
    print ""
    print "x4,y4"
    print x4_store
    print y4_store
    print ""
    print "x5,y5"
    print x5_store
    print y5_store
    print ""
    print "theta"
    print (theta*180.0/math.pi)
    print ""
    print "coordsum"
    print (x3 + y3 + x4_store + y4_store)
    print ""
    print "s1"
    print s1_store
    print ""
    print "s2"
    print s2_store


print ""
print "############################"
print ""

#################################

# Solution 1

#minres
#1.19864241774e-05

#x3,y3
#0.547534058653
#-1.62491176637

#x4,y4
#2.2796078788
#0.824561703903

#x5,y5
#1.41357096873
#-0.400175031233

#theta
#54.7350720305

#coordsum
#2.02679187499

#s1
#0.71468125139

#s2
#2.23712753794

#################################

# Solution 2

#minres
#1.62878861332e-05

#x3,y3
#2.00250143645
#1.18250689221

#x4,y4
#0.312444855334
#-1.29614369886

#x5,y5
#1.15747314589
#-0.0568184033262

#theta
#235.712057229

#coordsum
#2.20130948513

#s1
#1.32558262659

#s2
#4.84575621333

#################################

# Solution 3

#minres
#3.50611294441e-05

#x3,y3
#-0.00252373567027
#0.0318641434243

#x4,y4
#2.72529991218
#1.2804546072

#x5,y5
#1.36138808825
#0.656159375313

#theta
#384.594708655

#coordsum
#4.03509492714

#s1
#-0.968036069111

#s2
#1.72967533002

#################################

Consider the given system:

The total potential energy of the system is:

$V = V_{grav} + V_{k_1} + V_{k_2}$

$V = mg\left(\frac{y_3+y_4}{2}\right) + \frac{1}{2}k_1\left(\sqrt{x_3^2 + y_3^2}-L_1\right)^2 + \frac{1}{2}k_2\left(\sqrt{(x_4-3)^2 + (y_4-5)^2}-L_2\right)^2$

Also, the length of the rod is a constraint equation which is:

$C = \sqrt{(x_4-x_3)^2 + (y_4-y_3)^2}-L = 0$

The system is in equilibrium when the potential energy of the system is a minimum. In other words, a constrained optimization problem must be solved to obtain the equilibrium configuration.

The optimization problem is: Minimise:

$F = V + \lambda C$

Subject to $x_3>0$ and $y_3 <0$ . Here $\lambda$ is an unknown Lagrange multiplier.

The next step involves computing partial derivatives and equating them to zero. This gives five equations:

$\frac{\partial F}{\partial x_3} = 0$ $\frac{\partial F}{\partial y_3} = 0$ $\frac{\partial F}{\partial x_4} = 0$ $\frac{\partial F}{\partial y_4} = 0$ $\frac{\partial F}{\partial \lambda} = 0$

These five equations yield a system of nonlinear equations that have multiple solutions. The required solution would satisfy the additional constraints $x_3>0$ and $y_3 <0$ . Of course, the optimization problem could have been formulated accounting for the inequality constraints and thereby dealing with two additional Lagrange multipliers. The optimal solution would then need to satisfy the KKT conditions. However, this has not been done.

The manner in which these equations are solved is iterative. The algorithm used is the Newton-Raphson method the details of which are omitted as well. An appropriate start point must be chosen. Different start points are likely to yield different solutions. The solution satisfying the additional constraints are:

$\boxed{x_3 = 0.5475 \ ; \ y_3 = -1.6249 \ ; \ x_4 = 2.2796 \ ; \ y_4 = 0.8246}$

Of course, a second derivative test would need to be carried out to see whether the optimum point corresponds to a minimum. The system of equations would yield other solutions which are non-intuitive as well. These most likely correspond to unstable equilibria. In other words, the potential energy attains a local maximum value at those points.

If there happens to be an easier way of solving this problem, I am keen on knowing. For convenience, I have left out the gory details of the derivative crunching.