{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Worked problem: The string problem\n", "\n", "Ignore the following drawing code." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [ "input-hidden" ] }, "outputs": [ { "data": { "image/svg+xml": [ "\n", " \n", "\n", " \n", " \n", " \n", "\n", " \n", " \n", " L\n", " \n", " \n", " \n", " W12\n", " \n", " \n", " \n", " W23\n", " \n", " 1\n", " 2\n", " 3\n", "" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import SVG\n", "from collections import namedtuple\n", "\n", "V = namedtuple(\"V\", [\"x\", \"y\"])\n", "\n", "p0 = V(10, 10)\n", "p1 = V(50, 80)\n", "p2 = V(150, 100)\n", "p3 = V(290, 10)\n", "\n", "\n", "def sub(x):\n", " return f'{x}'\n", "\n", "\n", "SVG(\n", " f\"\"\"\n", "\n", " \n", "\n", " \n", " \n", " \n", "\n", " \n", " \n", " L\n", " \n", " \n", " \n", " W{sub(12)}\n", " \n", " \n", " \n", " W{sub(23)}\n", " \n", " 1\n", " 2\n", " 3\n", "\n", "\"\"\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll define $T_n$ to be the tension, $L_n$ to be the length, and $\\theta_n$ to be the angle from horizantal for each string labeled $n$.\n", "\n", "First, we assume the horizontal lengths add up to L:\n", "\n", "$$\n", "L_1 \\cos \\theta_1 + L_2 \\cos \\theta_2 + L_3 \\cos \\theta_3 = L \\tag{1}\n", "$$\n", "\n", "Then, we can assume that the vertical lengths cancel:\n", "\n", "$$\n", "L_1 \\sin \\theta_1 + L_2 \\sin \\theta_2 - L_3 \\sin \\theta_3 = 0\n", "\\tag{2}\n", "$$\n", "\n", "(We are explicilty defining theta to be the positive angle from the $x$ axis)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also use trigometric identites:\n", "\n", "$$\n", "\\begin{align}\n", "\\sin^2 \\theta_1 + \\cos^2 \\theta_1 &=& 1 \\tag{3} \\\\\n", "\\sin^2 \\theta_2 + \\cos^2 \\theta_2 &=& 1 \\tag{4} \\\\\n", "\\sin^2 \\theta_3 + \\cos^2 \\theta_3 &=& 1 \\tag{5}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can use physics to get 4 force equations, two for each weight.\n", "\n", "$$\n", "\\begin{align}\n", "T_1 \\sin \\theta_1 - T_2 \\sin \\theta_2 - W_{12} &=& 0 \\tag{6} \\\\\n", "T_1 \\cos \\theta_1 - T_2 \\cos \\theta_2 &=& 0 \\tag{7} \\\\\n", "T_2 \\sin \\theta_2 + T_3 \\sin \\theta_3 - W_{23} &=& 0 \\tag{8} \\\\\n", "T_2 \\cos \\theta_2 - T_3 \\cos \\theta_3 &=& 0 \\tag{9} \\\\\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we have our unknowns:\n", "\n", "$$\n", "\\mathbf{x} = \\left(\n", "\\begin{matrix}\n", "\\sin{\\theta_1} \\\\\n", "\\sin{\\theta_2} \\\\\n", "\\sin{\\theta_3} \\\\\n", "\\cos{\\theta_1} \\\\\n", "\\cos{\\theta_2} \\\\\n", "\\cos{\\theta_3} \\\\\n", "T_1 \\\\\n", "T_2 \\\\\n", "T_3 \\\\\n", "\\end{matrix}\n", "\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But, we have the problem that our solution is non-linear:\n", "\n", "$$\n", "f(\\mathbf{x}) = \\left(\n", "\\begin{matrix}\n", "3 x_5 + 4 x_4 + 4 x_5 - 8 \\\\\n", "3 x_0 + 4 x_1 + 4 x_2 \\\\\n", "x_6 x_0 - x_7 x_1 - 10 \\\\\n", "x_6 x_3 - x_7 x_4 \\\\\n", "x_7 x_1 + x_8 x_2 - 20 \\\\\n", "x_7 x_4 + x_8 x_6 \\\\\n", "x_0^2 + x_3^2 - 1 \\\\\n", "x_1^2 + x_4^2 - 1 \\\\\n", "x_2^2 + x_5^2 - 1 \\\\\n", "\\end{matrix}\n", "\\right) = 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Unlike the book, I'm using 0 based indexing. Let's use SymPy to give us some symbolic manipulation abilities:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import sympy as s\n", "\n", "s.init_printing()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "x = s.Matrix(s.symbols(\"x:10\"))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}3 x_{3} + 4 x_{4} + 4 x_{5} - 8\\\\3 x_{0} + 4 x_{1} - 4 x_{2}\\\\x_{0} x_{6} - x_{1} x_{7} - 10\\\\x_{3} x_{6} - x_{4} x_{7}\\\\x_{1} x_{7} + x_{2} x_{8} - 20\\\\x_{4} x_{7} - x_{5} x_{8}\\\\x_{0}^{2} + x_{3}^{2} - 1\\\\x_{1}^{2} + x_{4}^{2} - 1\\\\x_{2}^{2} + x_{5}^{2} - 1\\end{matrix}\\right]$" ], "text/plain": [ "⎡3⋅x₃ + 4⋅x₄ + 4⋅x₅ - 8⎤\n", "⎢ ⎥\n", "⎢ 3⋅x₀ + 4⋅x₁ - 4⋅x₂ ⎥\n", "⎢ ⎥\n", "⎢ x₀⋅x₆ - x₁⋅x₇ - 10 ⎥\n", "⎢ ⎥\n", "⎢ x₃⋅x₆ - x₄⋅x₇ ⎥\n", "⎢ ⎥\n", "⎢ x₁⋅x₇ + x₂⋅x₈ - 20 ⎥\n", "⎢ ⎥\n", "⎢ x₄⋅x₇ - x₅⋅x₈ ⎥\n", "⎢ ⎥\n", "⎢ 2 2 ⎥\n", "⎢ x₀ + x₃ - 1 ⎥\n", "⎢ ⎥\n", "⎢ 2 2 ⎥\n", "⎢ x₁ + x₄ - 1 ⎥\n", "⎢ ⎥\n", "⎢ 2 2 ⎥\n", "⎣ x₂ + x₅ - 1 ⎦" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f = s.Matrix(\n", " [\n", " 3 * x[3] + 4 * x[4] + 4 * x[5] - 8,\n", " 3 * x[0] + 4 * x[1] - 4 * x[2],\n", " x[6] * x[0] - x[7] * x[1] - 10,\n", " x[6] * x[3] - x[7] * x[4],\n", " x[7] * x[1] + x[8] * x[2] - 20,\n", " x[7] * x[4] - x[8] * x[5],\n", " x[0] ** 2 + x[3] ** 2 - 1,\n", " x[1] ** 2 + x[4] ** 2 - 1,\n", " x[2] ** 2 + x[5] ** 2 - 1,\n", " ]\n", ")\n", "f" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 3 & 4 & 4 & 0 & 0 & 0\\\\3 & 4 & -4 & 0 & 0 & 0 & 0 & 0 & 0\\\\x_{6} & - x_{7} & 0 & 0 & 0 & 0 & x_{0} & - x_{1} & 0\\\\0 & 0 & 0 & x_{6} & - x_{7} & 0 & x_{3} & - x_{4} & 0\\\\0 & x_{7} & x_{8} & 0 & 0 & 0 & 0 & x_{1} & x_{2}\\\\0 & 0 & 0 & 0 & x_{7} & - x_{8} & 0 & x_{4} & - x_{5}\\\\2 x_{0} & 0 & 0 & 2 x_{3} & 0 & 0 & 0 & 0 & 0\\\\0 & 2 x_{1} & 0 & 0 & 2 x_{4} & 0 & 0 & 0 & 0\\\\0 & 0 & 2 x_{2} & 0 & 0 & 2 x_{5} & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "⎡ 0 0 0 3 4 4 0 0 0 ⎤\n", "⎢ ⎥\n", "⎢ 3 4 -4 0 0 0 0 0 0 ⎥\n", "⎢ ⎥\n", "⎢ x₆ -x₇ 0 0 0 0 x₀ -x₁ 0 ⎥\n", "⎢ ⎥\n", "⎢ 0 0 0 x₆ -x₇ 0 x₃ -x₄ 0 ⎥\n", "⎢ ⎥\n", "⎢ 0 x₇ x₈ 0 0 0 0 x₁ x₂ ⎥\n", "⎢ ⎥\n", "⎢ 0 0 0 0 x₇ -x₈ 0 x₄ -x₅⎥\n", "⎢ ⎥\n", "⎢2⋅x₀ 0 0 2⋅x₃ 0 0 0 0 0 ⎥\n", "⎢ ⎥\n", "⎢ 0 2⋅x₁ 0 0 2⋅x₄ 0 0 0 0 ⎥\n", "⎢ ⎥\n", "⎣ 0 0 2⋅x₂ 0 0 2⋅x₅ 0 0 0 ⎦" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = s.Matrix([f.diff(x[i]).T for i in range(9)]).T\n", "df" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of iterations = 8\n", "Final Solution:\n", "x_arr[0] = 0.7610026921032024\n", "x_arr[1] = 0.2649538102798601\n", "x_arr[2] = 0.8357058293572619\n", "x_arr[3] = 0.648748720702048\n", "x_arr[4] = 0.9642611048977284\n", "x_arr[5] = 0.5491773545757356\n", "x_arr[6] = 17.16020978458048\n", "x_arr[7] = 11.54527968431142\n", "x_arr[8] = 20.271528044627093\n" ] } ], "source": [ "x_arr = np.array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0, 1.0])\n", "eps = 1e-3\n", "\n", "for it in range(15):\n", "\n", " # Compute f and its derivative\n", " y = f.evalf(subs={a: b for a, b in zip(x, x_arr)})\n", " M = df.evalf(subs={a: b for a, b in zip(x, x_arr)})\n", "\n", " # Convert these SymPy contraptions into Numpy\n", " y = np.asarray(y).astype(np.float64).flatten()\n", " M = np.asarray(M).astype(np.float64)\n", "\n", " # Solve for Δx\n", " Δx = np.linalg.solve(M, -y)\n", "\n", " # Compute new x\n", " x_arr += Δx\n", "\n", " errX = abs(Δx)\n", " errX[x_arr != 0] /= abs(x_arr[x_arr != 0]) # Relitave error only if x is not 0\n", " errF = abs(y)\n", "\n", " if np.all(errX <= eps) and np.all(errF <= eps):\n", " break\n", "\n", "print(\"Number of iterations = \", it + 1)\n", "print(\"Final Solution:\")\n", "for i in range(9):\n", " print(f\"x_arr[{i}] = {x_arr[i]}\")" ] } ], "metadata": { "kernelspec": { "display_name": "compclass", "language": "python", "name": "compclass" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }