Two masses $M_1$ and $M_2$ are hung from a horizontal rod with length $L$ in such a way that a rope of length $L_1$ connects the left end of the rod to $M_1$, a rope of length $L_2$ connects $M_1$ and $M_2$, and a rope of length $L_3$ connects $M_2$ to the right end of the rod. The system is at rest (in equilibrium under gravity).
Find the angles that the ropes make with the rod and the tension forces in the ropes.
Treat $\sin\theta_i$ and $\cos\theta_i$ together with $T_i$, $1\leq i \leq 3$, as unknowns that have to simultaneously fulfill the nine equations \begin{align} -T_1 \cos\theta_1 + T_2\cos\theta_2 &= 0\\ T_1 \sin\theta_1 - T_2\sin\theta_2 - W_1 &= 0\\ -T_2\cos\theta_2 + T_3\cos\theta_3 &= 0\\ T_2\sin\theta_2 + T_3\sin\theta_3 - W_2 &= 0\\ L_1\cos\theta_1 + L_2\cos\theta_2 + L_3\cos\theta_3 - L &= 0\\ -L_1\sin\theta_1 -L_2\sin\theta_2 + L_3\sin\theta_3 &= 0\\ \sin^2\theta_1 + \cos^2\theta_1 - 1 &= 0\\ \sin^2\theta_2 + \cos^2\theta_2 - 1 &= 0\\ \sin^2\theta_3 + \cos^2\theta_3 - 1 &= 0 \end{align}
Consider the nine equations a vector function $\mathbf{f}$ that takes a 9-vector $\mathbf{x}$ of the unknowns as argument: \begin{align} \mathbf{f}(\mathbf{x}) &= 0\\ \mathbf{x} &= \left(\begin{array}{c} x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \\ x_6 \\ x_7 \\ x_8 \end{array}\right) = \left(\begin{array}{c} \sin\theta_1 \\ \sin\theta_2 \\ \sin\theta_3 \\ \cos\theta_1 \\ \cos\theta_2 \\ \cos\theta_3 \\ T_1 \\ T_2 \\ T_3 \end{array}\right) \\ \mathbf{L} &= \left(\begin{array}{c} L \\ L_1 \\ L_2 \\ L_3 \end{array}\right), \quad \mathbf{W} = \left(\begin{array}{c} W_1 \\ W_2 \end{array}\right) \end{align}
Using the unknowns from above, our system of 9 coupled equations is:
\begin{align} -x_6 x_3 + x_7 x_4 &= 0\\ x_6 x_0 - x_7 x_1 - W_1 &= 0\\ -x_7x_4 + x_8 x_5 &= 0\\ x_7x_1 + x_8 x_2 - W_2 &= 0\\ L_1x_3 + L_2 x_4 + L_3 x_5 - L &= 0\\ -L_1x_0 - L_2 x_1 + L_3 x_2 &= 0\\ x_{0}^{2} + x_{3}^{2} - 1 &= 0\\ x_{1}^{2} + x_{4}^{2} - 1 &= 0\\ x_{2}^{2} + x_{5}^{2} - 1 &= 0 \end{align}Solve the root-finding problem $\mathbf{f}(\mathbf{x}) = 0$ with the generalized Newton-Raphson algorithm:
$$ \mathsf{J}(\mathbf{x}) \Delta\mathbf{x} = -\mathbf{f}(\mathbf{x}) $$and
$$ \mathbf{x} \leftarrow \mathbf{x} + \Delta\mathbf{x}. $$Set the problem parameters and the objective function $\mathbf{f}(\mathbf{x})$
import numpy as np
# problem parameters
W = np.array([10, 20])
L = np.array([8, 3, 4, 4])
def f_2masses(x, L, W):
return np.array([
-x[6]*x[3] + x[7]*x[4],
x[6]*x[0] - x[7]*x[1] - W[0],
-x[7]*x[4] + x[8]*x[5],
x[7]*x[1] + x[8]*x[2] - W[1],
L[1]*x[3] + L[2]*x[4] + L[3]*x[5] - L[0],
-L[1]*x[0] - L[2]*x[1] + L[3]*x[2],
x[0]**2 + x[3]**2 - 1,
x[1]**2 + x[4]**2 - 1,
x[2]**2 + x[5]**2 - 1,
])
def fLW(x, L=L, W=W):
return f_2masses(x, L, W)
Guess some initial values (they don't have to fullfil the equations!):
# initial parameters
#theta0 = np.deg2rad([45, 45, 90])
#T0 = np.array([1, 1, 2])
#x0 = np.concatenate([np.sin(theta0), np.cos(theta0), T0])
x0 = np.array([1.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1, 1])
x0
array([1.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1. , 1. , 1. ])
f_2masses(x0, L, W)
array([ 0. , -9. , 0. , -19. , -2.5, -4.5, 1.5, -0.5, -0.5])
Plot the positions of the 2 masses and the 3 strings for any solution vector $\mathbf{x}$:
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
def plot_2masses(x, L, W, **kwargs):
"""Plot 2 mass/3 string problem for parameter vector x and parameters L and W"""
kwargs.setdefault('linestyle', '-')
kwargs.setdefault('marker', 'o')
kwargs.setdefault('linewidth', 1)
ax = kwargs.pop('ax', None)
if ax is None:
ax = plt.subplot(111)
r0 = np.array([0, 0])
r1 = r0 + np.array([L[0], 0])
rod = np.transpose([r0, r1])
L1 = r0 + np.array([L[1]*x[3], -L[1]*x[0]])
L2 = L1 + np.array([L[2]*x[4], -L[2]*x[1]])
L3 = L2 + np.array([L[3]*x[5], L[3]*x[2]])
strings = np.transpose([r0, L1, L2, L3])
ax.plot(rod[0], rod[1], color="black", marker="d", linewidth=4)
ax.plot(strings[0], strings[1], **kwargs)
ax.set_aspect(1)
return ax
What does the initial guess look like?
plot_2masses(x0, L, W)
<Axes: >
Write a function Jacobian(f, x, h=1e-5)
that computes the Jacobian matrix numerically (use the central difference algorithm).
def Jacobian(f, x, h=1e-5):
"""df_i/dx_j with central difference (fi(xj+h/2)-fi(xj-h/2))/h"""
J = np.zeros((len(f(x)), len(x)), dtype=np.float64)
hvec = np.zeros_like(x, dtype=np.float64)
for j in range(len(x)):
hvec *= 0
hvec[j] = 0.5*h
J[:, j] = (f(x + hvec) - f(x - hvec))/h
return J
Test Jacobian()
on
$$
\mathbf{g}(\mathbf{x}) = \begin{pmatrix}
x_0^2 - x_1 \\ x_0
\end{pmatrix}
$$
with analytical result
Given a test vector $\mathbf{x}_\text{test} = (1, 0)$, what is the numerical answer for $\mathsf{J}(\mathbf{x}_\text{test})$?
$$ \mathsf{J}(\mathbf{x}_\text{test}) = \begin{pmatrix} 2 & -1\\ 1 & 0 \end{pmatrix} $$Test your Jacobian()
function with $\mathbf{x}_\text{test}$ and check that you get the same answer:
def g(x):
return np.array([
x[0]**2 - x[1],
x[0]
])
x_test = np.array([1, 0])
J = Jacobian(g, x_test)
print(J)
[[ 2. -1.] [ 1. 0.]]
Test that it also works for our starting vector:
J0 = Jacobian(fLW, x0)
J0
array([[ 0. , 0. , 0. , -1. , 1. , 0. , -0.5, 0.5, 0. ], [ 1. , -1. , 0. , 0. , 0. , 0. , 1.5, -0.5, 0. ], [ 0. , 0. , 0. , 0. , -1. , 1. , 0. , -0.5, 0.5], [ 0. , 1. , 1. , 0. , 0. , 0. , 0. , 0.5, 0.5], [ 0. , 0. , 0. , 3. , 4. , 4. , 0. , 0. , 0. ], [-3. , -4. , 4. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 3. , 0. , 0. , 1. , 0. , 0. , 0. , 0. , 0. ], [ 0. , 1. , 0. , 0. , 1. , 0. , 0. , 0. , 0. ], [ 0. , 0. , 1. , 0. , 0. , 1. , 0. , 0. , 0. ]])
J0.shape
(9, 9)
Write a function newton_raphson(f, x, Nmax=100, tol=1e-8, h=1e-5)
to find a root for a vector function f(x)=0
. (See also 12 Root-finding by trial-and-error and the 1D Newton-Raphson algorithm in 12-Root-finding.ipynb.) As a convergence criterion we demand that the length of the vector f(x)
(the norm --- see np.linalg.norm
) be less than the tolerance.
def newton_raphson(f, x, Nmax=100, tol=1e-8, h=1e-5):
"""n-D Newton-Raphson: solves f(x) = 0.
Iterate until |f(x)| < tol or nmax steps.
"""
x = x.copy()
for istep in range(Nmax):
fx = f(x)
if np.linalg.norm(fx) < tol:
break
J = Jacobian(f, x, h=h)
Delta_x = np.linalg.solve(J, -fx)
x += Delta_x
else:
raise RuntimeError(f"Newton-Raphson: no root found after {istep} iterations"
f"(eps={tol}); best guess is {x} with error {np.linalg.norm(fx)}")
return x
x = newton_raphson(fLW, x0)
print(x0)
print(x)
[1.5 0.5 0.5 0.5 0.5 0.5 1. 1. 1. ] [ 0.76100269 0.26495381 0.83570583 0.64874872 0.9642611 0.54917735 17.16020978 11.54527968 20.27152804]
Plot the starting configuration and the solution:
ax = plot_2masses(x0, L, W)
ax = plot_2masses(x, L, W, ax=ax)
Pretty-print the solution (angles in degrees):
def pretty_print(x):
theta = np.rad2deg(np.arcsin(x[0:3]))
tensions = x[6:]
print("theta1 = {0[0]:.1f} \t theta2 = {0[1]:.1f} \t theta3 = {0[2]:.1f}".format(theta))
print("T1 = {0[0]:.1f} \t T2 = {0[1]:.1f} \t T3 = {0[2]:.1f}".format(tensions))
print("Starting values")
pretty_print(x0)
print()
print("Solution")
pretty_print(x)
Starting values theta1 = nan theta2 = 30.0 theta3 = 30.0 T1 = 1.0 T2 = 1.0 T3 = 1.0 Solution theta1 = 49.6 theta2 = 15.4 theta3 = 56.7 T1 = 17.2 T2 = 11.5 T3 = 20.3
/var/folders/sm/37rm_wm16tq9qsf4n5md98ph0000gp/T/ipykernel_82281/3554523550.py:2: RuntimeWarning: invalid value encountered in arcsin theta = np.rad2deg(np.arcsin(x[0:3]))
Create a new function newton_raphson_intermediates()
based on newton_raphson()
that returns all trial x
values including the last one.
def newton_raphson_intermediates(f, x, Nmax=100, tol=1e-8, h=1e-5):
"""n-D Newton-Raphson: solves f(x) = 0.
Iterate until |f(x)| < tol or nmax steps.
Returns all intermediates.
"""
intermediates = []
x = x.copy()
for istep in range(Nmax):
fx = f(x)
if np.linalg.norm(fx) < tol:
break
J = Jacobian(f, x, h=h)
Delta_x = np.linalg.solve(J, -fx)
intermediates.append(x.copy())
x += Delta_x
else:
raise RuntimeError(f"Newton-Raphson: no root found after {istep} iterations"
f"(eps={tol}); best guess is {x} with error {np.linalg.norm(fx)}")
return np.array(intermediates)
Visualize the intermediate configurations:
x_series = newton_raphson_intermediates(fLW, x0)
ax = plt.subplot(111)
ax.set_prop_cycle("color", [plt.cm.viridis_r(i) for i in np.linspace(0, 1, len(x_series))])
for iteration, x in enumerate(x_series):
plot_2masses(x, L, W, label=str(iteration), ax=ax)
ax.legend(loc="best");
It's convenient to turn the above plotting code into a function that we can reuse:
def plot_series(x_series, L, W):
"""Plot all N masses/strings solution vectors in x_series (N, 9) array"""
ax = plt.subplot(111)
ax.set_prop_cycle("color", [plt.cm.viridis_r(i) for i in np.linspace(0, 1, len(x_series))])
for iteration, x in enumerate(x_series):
plot_2masses(x, L, W, label=str(iteration), ax=ax)
ax.legend(loc="best")
return ax
Try different masses, e.g. M1 = M2 = 10, or M1= 0 , M2 = 10.
Use nicer starting parameters that already fulfill the angle equations (7) - (9) (but it works with pretty much any guess):
# initial parameters
theta0 = np.deg2rad([45, 45, 90])
T0 = np.array([1, 1, 2])
x0 = np.concatenate([np.sin(theta0), np.cos(theta0), T0])
W_2 = np.array([10, 10])
def fLW_2(x, L=L, W=W_2):
return f_2masses(x, L, W)
x_series_2 = newton_raphson_intermediates(fLW_2, x0)
pretty_print(x_series_2[-1])
theta1 = 57.9 theta2 = 8.8 theta3 = 52.1 T1 = 13.1 T2 = 7.0 T3 = 11.3
plot_series(x_series_2, L, W_2)
<Axes: >
W_3 = np.array([0, 10])
def fLW_3(x):
return f_2masses(x, L=L, W=W_3)
x_series_3 = newton_raphson_intermediates(fLW_3, x0)
pretty_print(x_series_3[-1])
theta1 = 30.0 theta2 = 30.0 theta3 = 61.0 T1 = 4.8 T2 = 4.8 T3 = 8.7
plot_series(x_series_3, L, W_3)
<Axes: >