Project 2 - Mechanics with SymPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project 2 - Mechanics with SymPy\n",
    "\n",
    "![SymPy logo](https://www.sympy.org/static/images/logo.png)\n",
    "\n",
    "**SymPy** is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. It is available under a permissive, commercial-friendly open source license (BSD 3-Clause).\n",
    "\n",
    "Some of its capabilities are:\n",
    "\n",
    "- Arithmetic, simplification, expansion, substitution\n",
    "- Limits, differentiation, integration, Taylor series\n",
    "- Equation solving\n",
    "- Matrix algebra\n",
    "- Physics\n",
    "\n",
    "SymPy is one of the most active Python projects in terms of contributors.\n",
    "\n",
    "- **Documentation**: https://docs.sympy.org/latest/index.html\n",
    "- **Source code**: https://github.com/sympy/sympy/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basics of symbolic computation\n",
    "\n",
    "(_From the_ [_SymPy tutorial_](https://docs.sympy.org/latest/tutorial/intro.html))\n",
    "\n",
    "Symbolic computation deals with the computation of mathematical objects symbolically. This means that the mathematical objects are represented exactly, not approximately, and mathematical expressions with unevaluated variables are left in symbolic form.\n",
    "\n",
    "Let’s take an example. Say we wanted to use the built-in Python functions to compute square roots. We might do something like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "math.sqrt(9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "9 is a perfect square, so we got the exact answer, 3. But suppose we computed the square root of a number that isn’t a perfect square:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "math.sqrt(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But suppose we want to go further. Recall that $\\sqrt{8} = \\sqrt{4 \\cdot 2} = 2 \\sqrt{2}$. We would have a hard time deducing this from the above result. This is where symbolic computation comes in. With a symbolic computation system like SymPy, square roots of numbers that are not perfect squares are left unevaluated by default:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.sqrt(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore—and this is where we start to see the real power of symbolic computation—symbolic results can be symbolically simplified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.sqrt(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A more interesting example\n",
    "\n",
    "Symbolic computation systems (which by the way, are also often called computer algebra systems, or just CASs) such as SymPy are capable of computing symbolic expressions with variables.\n",
    "\n",
    "Let us define a symbolic expression, representing the mathematical expression $x+2y$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = sp.symbols(\"x y\")\n",
    "expr = x + 2 * y\n",
    "expr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we wrote `x + 2*y` just as we would if `x` and `y` were ordinary Python variables. But in this case, instead of evaluating to something, the expression remains as just `x + 2*y`. Now let us play around with it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expr - x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x * expr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we might have expected $x(x+2y)$ to transform into $x^2+2xy$, but instead we see that the expression was left alone. This is a common theme in SymPy. Aside from obvious simplifications like $x−x=0$ and $\\sqrt{8} = 2\\sqrt{2}$, most simplifications are not performed automatically. This is because we might prefer the factored form $x(x+2y)$, or we might prefer the expanded form $x^2+2xy$. Both forms are useful in different circumstances. In SymPy, there are functions to go from one form to the other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expanded_expr = sp.expand(x * expr)\n",
    "expanded_expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.factor(expanded_expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many more things that can be done: differentiating and integrating expressions with respect to variables, computing limits, solving algebraic and differential equations... Check out the [SymPy tutorial](https://docs.sympy.org/latest/tutorial/index.html) for more!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mechanics\n",
    "\n",
    "### Reference frames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To study Mechanics and Kinematics using SymPy, we need to start by defining our own reference frames. Then we will use them to express our geometrical relations and transform vectors between them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()\n",
    "\n",
    "from sympy.physics.mechanics import ReferenceFrame\n",
    "\n",
    "# We define our own class so our unit vectors are (i, j, k)\n",
    "class IJKReferenceFrame(ReferenceFrame):\n",
    "    def __init__(self, name):\n",
    "        super().__init__(name, latexs=['\\mathbf{%s}_{%s}' % (idx, name) for idx in (\"i\", \"j\", \"k\")])\n",
    "        self.i = self.x\n",
    "        self.j = self.y\n",
    "        self.k = self.z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define our own reference frame like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = IJKReferenceFrame(\"1\")\n",
    "A.i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And to write vectors in this frame we **multiply each component by the corresponding unit vector**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2 * A.i - 1 * A.j"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vector algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our vectors work like symbols, and we can perform dot and cross products with them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R, V = sp.symbols(\"R, V\", positive=True)\n",
    "r1 = R * (A.i + A.j + A.k)\n",
    "v1 = V * (A.i - 2 * A.k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r1.dot(v1)\n",
    "r1 & v1  # Alternative to .dot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r1.cross(v1)\n",
    "r1 ^ v1  # Alternative to .cross"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute their `.magnitude()` and `.normalize()` them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(r1 ^ v1).magnitude()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(r1 ^ v1).normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Using directly the Transport Theorem, Basic Kinematic Equation or formula for time derivative in a rotating frame:\n",
    "\n",
    "$$\\left(\\frac{\\operatorname{d}\\!\\mathbf{a}}{\\operatorname{d}\\!t}\\right)_1 = \\left(\\frac{\\operatorname{d}\\!\\mathbf{a}}{\\operatorname{d}\\!t}\\right)_0 + \\mathbf{\\omega}_{01}\\! \\times \\mathbf{a}$$\n",
    "\n",
    "Compute the derivative of this vector:\n",
    "\n",
    "$$\n",
    "\\mathbf{a} = R \\mathbf{i}_0\n",
    "$$\n",
    "\n",
    "where $A_0$ is a reference frame that rotates with respect to the inertial frame with angular velocity $\\mathbf{\\omega}_{01}=\\Omega \\mathbf{k}_0$. **What's the magnitude of the derivative?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your solution here!\n",
    "# Remember to insert extra cells if you need it\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dynamic symbols\n",
    "\n",
    "There is a straightforward way of creating symbols that have a time dependence using `dynamicsymbols`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.physics.mechanics import dynamicsymbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "α = dynamicsymbols(\"α\")\n",
    "α"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And compute its derivative using `.diff()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "α.diff()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relative motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SymPy provides powerful methods to compute rotation matrices between two reference frames. For that we need to specify their relative orientation using `.orient`, and then we can recover the direct cosine matrix using `.dcm`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A0 = IJKReferenceFrame(\"0\")\n",
    "A1 = IJKReferenceFrame(\"1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ϕ = dynamicsymbols(\"ϕ\")\n",
    "A1.orient(A0, \"Axis\", [ϕ, A0.z])  # Rotation ϕ around A1.z axis\n",
    "A1.dcm(A0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can obtain its angular velocity using `.ang_vel_in`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A1.ang_vel_in(A0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using value `Axis` we specified the rotation around an axis. Other methods are:\n",
    "\n",
    "* `Body`: three Euler angles.\n",
    "* `Space`: like `Body`, but the rotations are applied in the reverse order.\n",
    "* `Quaternion`: rotating around a unit vector $\\mathbf{\\lambda}$ a quantity $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To express a vector in a different reference frame, we can use `.express` or `.to_matrix`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A0.i.express(A1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A0.i.to_matrix(A1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time derivative in rotating frame\n",
    "\n",
    "Rather than applying the Basic Kinematic Equation manually as we did above, SymPy can do that automatically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v1 = A1.x\n",
    "v1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dv1 = v1.diff(sp.symbols(\"t\"), A0)\n",
    "dv1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dv1.to_matrix(A1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(dv1 & A1.j).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "![Reference frames of a rotor blade](img/rotor-blade.jpg)\n",
    "\n",
    "(Cuerva et al. \"Teoría de los Helicópteros\")\n",
    "\n",
    "1. Obtain the direct cosine matrix of the blade $B$ with respect to $A1$.\n",
    "2. Obtain the angular velocity of $B$ with respect to $A1$, expressed in $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your solution here!\n",
    "# Remember to insert extra cells if you need it\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Points and velocities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last step in our kinematics study is the possibility to define points in solids and compute their velocity field. For this, SymPy provides the `Point` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.physics.mechanics import Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "O = Point(\"O\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To set a point $O$ as the origin of $A$, we will define its velocity to be zero with `.set_vel`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = IJKReferenceFrame(\"A\")\n",
    "A1 = IJKReferenceFrame(\"A1\")\n",
    "ψ = dynamicsymbols(\"ψ\")\n",
    "A1.orient(A, \"Axis\", [ψ, A.z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "O.set_vel(A, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, to define new points, we use the method `.locate_new`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e_b = sp.symbols(\"e_b\")\n",
    "E_b = O.locatenew(\"E_b\", e_b * A1.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain the position vector of one point with respect to the other, we use `.pos_from`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_b.pos_from(O)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, the **velocity field of a rigid body** is obtained using `.v2pt_theory`:\n",
    "\n",
    "$$v^P_A = v^O_A + \\omega_{A_1 A} \\times \\mathbf{OP}$$\n",
    "\n",
    "This method belongs to _the point we want to study_ and receives three parameters:\n",
    "\n",
    "* `O`: point with known velocity in $A$\n",
    "* `A`: reference frame where we want to compute the velocity\n",
    "* `A1`: reference frame where both points are fixed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, to compute the velocity of the point we just created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_b.v2pt_theory(O, A, A1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The no-slip disc\n",
    "\n",
    "![No-slip disc](img/no-slip-disc.png)\n",
    "\n",
    "(Notes from Óscar López Rebollal, Polytechnic University of Madrid)\n",
    "\n",
    "**Compute the velocity and acceleration of point $P$**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A1 = IJKReferenceFrame(\"1\")\n",
    "A0 = IJKReferenceFrame(\"0\")\n",
    "A2 = IJKReferenceFrame(\"2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ξ, θ = dynamicsymbols(\"ξ, θ\")\n",
    "ξ, θ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A0.orient(A1, \"Axis\", [0, A1.k])  # A0 does not rotate with respect to A1\n",
    "A2.orient(A0, \"Axis\", [θ, A0.k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A2.dcm(A1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = Point(\"C\")  # Center of the disc\n",
    "C.set_vel(A1, ξ.diff() * A1.i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We locate P, fixed point on the disc, with respect to C, in A2 (rotating with the disc)\n",
    "R = sp.symbols(\"R\", positive=True)\n",
    "P = C.locatenew(\"P\", -R * A2.j)\n",
    "P.pos_from(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Velocity of P in A1\n",
    "# With this function call, we express that C and P are fixed in A2\n",
    "P.v2pt_theory(C, A1, A2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P.v2pt_theory(C, A1, A2).express(A0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Mission accomplished 🚀**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "- Vector documentation (kinematics) https://docs.sympy.org/latest/modules/physics/vector/index.html\n",
    "- Rigid bodies and more https://docs.sympy.org/latest/modules/physics/mechanics/masses.html\n",
    "- Dynamics, linearization... https://docs.sympy.org/latest/modules/physics/mechanics/\n",
    "\n",
    "(_This material is based on [Curso AeroPython by Juan Luis Cano Rodriguez and Alejandro Sáez Mollejo](https://github.com/AeroPython/Curso_AeroPython) (CC-BY)_)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}