{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Le groupe orthogonal de $\\mathbb R^3$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marc Lorenzi - juillet 2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0 Préliminaires"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On se place dans l'espace euclidien orienté $E=\\mathbb R^3$. Les endomorphismes orthogonaux de $E$ sont de l'un des types suivants :\n",
    "1. Les réflexions, symétries orthogonales par rapport à un plan, caractérisées par un vecteur normal au plan de la réflexion.\n",
    "2. Les rotations. L'identité mise à part, elles sont caractérisées par un axe, orienté par un vecteur $e$, et un angle $\\theta\\not\\equiv 0[2\\pi]$.\n",
    "3. Les composées d'une réflexion et d'une rotation différente de l'identité, dont l'axe est orthogonal au plan de la réflexion. L'axe et l'angle de la rotation caractérisent complètement ces endomorphismes orthogonaux.\n",
    "\n",
    "On se propose de répondre aux deux questions suivantes :\n",
    "1. Connaissant les éléments caractéristiques d'un endomorphisme orthogonal, déterminer sa matrice dans la base canonique.\n",
    "2. Connaissant la matrice d'un endomorphisme orthogonal, déterminer sa nature et ses éléments caractéristiques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Définissons tout d'abord les vecteurs de la base canonique de $\\mathbb R^3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e1 = Matrix([1, 0, 0])\n",
    "e2 = Matrix([0, 1, 0])\n",
    "e3 = Matrix([0, 0, 1])\n",
    "e1, e2, e3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et voici une petite fonction renvoyant un vecteur unitaire colinéaire à un vecteur non nul $u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normaliser(u):\n",
    "    return u / u.norm()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normaliser(Matrix([1, 2, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Les réflexions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Image d'un vecteur par une réflexion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $f$ la réflexion de plan $P$, où $P$ est donné par un vecteur normal $e$. On a alors pour tout $u\\in\\mathbb R^3$, $f(u)=u-2\\frac{<u,e>}{||e||^2}e$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reflexion(e, u):\n",
    "    return u - 2 * u.dot(e) / e.norm() ** 2 * e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reflexion(Matrix([1,1,1]), Matrix([1,2,1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Matrice d'une réflexion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour obtenir la matrice d'une réflexion, il suffit de calculer les images des vecteurs de la base canonique par celle-ci, puis de les ranger dans une matrice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrice_reflexion(e):\n",
    "    u1 = reflexion(e, e1)\n",
    "    u2 = reflexion(e, e2)\n",
    "    u3 = reflexion(e, e3)\n",
    "    return u1.col_insert(1, u2).col_insert(2, u3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A0 = matrice_reflexion(Matrix([1, 2, 3]))\n",
    "A0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelques vérifications ... la matrice $A_0$ est symétrique, ce qui est une bonne chose. Ensuite,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A0 * A0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "et donc $A_0$ est une matrice de symétrie orthogonale. Enfin,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "det(A0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "donc $A_0$ est une matrice de réflexion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Caractéristiques d'une réflexion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Donnons-nous une matrice $A$. La matrice $A$ est-elle une matrice de réflexion ? Et si oui, par rapport à quel plan ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def est_reflexion(A):\n",
    "    return A != - eye(3) and A == A.T and A.T * A == eye(3) and det(A) == -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_reflexion(A0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyse_reflexion(A):\n",
    "    return (A + eye(3)).nullspace()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse_reflexion(A0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Rotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Image d'un vecteur par une rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Une rotation $f$ de l'espace différente de l'identité est caractérisée par son axe, orienté par un vecteur $\\omega$ et son angle $\\theta\\not\\equiv 0[2\\pi]$. En prenant $\\omega$ unitaire, on a pour tout vecteur $u$ de l'espace la formule d'Euler-Rodrigues :\n",
    "\n",
    "$f(u)=\\cos\\theta u + <\\omega, u> (1 - \\cos\\theta)\\omega + \\sin\\theta\\ \\omega\\land u$\n",
    "\n",
    "La fonction __rotation__ renvoie l'image du vecteur u par la rotation caractérisée par $\\omega$ unitaire et l'angle $\\theta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotation(omega, theta, u):\n",
    "    c = cos(theta)\n",
    "    s = sin(theta)\n",
    "    v1 = c * u\n",
    "    v2 = s * omega.cross(u)\n",
    "    v3 = omega.dot(u) * (1 - c) * omega\n",
    "    return v1 + v2 + v3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, z = symbols('x y z')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation(normaliser(Matrix([1, 1, 1])), pi/2, Matrix([x, y, z]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simplify(_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour obtenir la matrice de la rotation, il suffit de calculer les images des vecteurs de la base canonique par celle-ci."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrice_rotation(omega, theta):\n",
    "    omega = normaliser(omega)\n",
    "    u = rotation(omega, theta, e1)\n",
    "    v = rotation(omega, theta, e2)\n",
    "    w = rotation(omega, theta, e3)\n",
    "    return u.col_insert(1, v).col_insert(2, w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = matrice_rotation(Matrix([1, 1, 1]), pi / 2)\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 2.2 Déterminer les caractéristiques d'une rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est facile de savoir si une matrice $A$ est la matrice d'une rotation différente de l'identité."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def est_rotation(A):\n",
    "    return simplify(A.T * A) == eye(3) and A != eye(3) and simplify(det(A)) == 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'axe de la rotation est l'ensemble des invariants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def invariants(A):\n",
    "    return (A - eye(3)).nullspace()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "invariants(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e=simplify(_)\n",
    "e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le cosinus de l'angle est donné par la trace de la rotation : $Tr\\ A = 1 + 2\\cos\\theta$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cos_angle(A):\n",
    "    return (A.trace() - 1) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cos_angle(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour obtenir le signe du sinus, on prend un vecteur $u$ qui n'est pas sur l'axe et on calcule le déterminant de $u,f(u),e$ où $e$ est un vecteur qui dirige l'axe (pas forcément unitaire). Le signe de $\\sin\\theta$ est alors le signe de ce daterminant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def signe_sin_angle(A, e):\n",
    "    if e.cross(Matrix([1, 0, 0])) != 0:\n",
    "        u = Matrix([1, 0, 0])\n",
    "    else:\n",
    "        u = Matrix([0, 1, 0])\n",
    "    v = A * u\n",
    "    M = Matrix([u.T, v.T, e.T]).T\n",
    "    return M.det()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signe_sin_angle(A, e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On regroupe le tout en une seule fonction. La fonction __analyse_rotation__ prend en paramètre une matrice $A$ censée être une matrice de rotation différente de l'identité. Elle renvoie le couple $(e,\\theta)$ où $e$ est un vecteur de l'axe (pas nécessairement unitaire) et $\\theta$ est l'angle de la rotation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyse_rotation(A):\n",
    "    e = simplify(invariants(A))\n",
    "    c = cos_angle(A)\n",
    "    theta = acos(c)\n",
    "    s = signe_sin_angle(A, e)\n",
    "    if s >= 0: return (e, theta)\n",
    "    else: return (e, -theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse_rotation(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Essayons deux autres exemples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A2 = Matrix([[3, 1, sqrt(6)],[1, 3, -sqrt(6)],[-sqrt(6),sqrt(6), 2]])/4\n",
    "A2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_rotation(A2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse_rotation(A2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A3 = Matrix([[8,1,-4],[-4,4,-7],[1,8,4]])/9\n",
    "A3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_rotation(A3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse_rotation(A3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquons pour terminer que l'identité, qui est une rotation, n'a pas d'axe défini de façon unique. Cela dit, notre fonction renvoie des résultats tout à fait logiques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse_rotation(eye(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Composée d'une réflexion et d'une rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $f$ un endomorphisme orthogonal de l'espace qui n'est ni ni une rotation ni une réflexion et qui est différent de $-id$. Alors $f$ s'crit de façon unique $g\\circ h$ où $g$ est une rotation différente de l'identité, $h$ est une réflexion, et le plan de la réflexion $h$ est orthogonal à l'axe de la rotation $g$. De plus, $g$ et $h$ commutent.\n",
    "\n",
    "Pour obtenir les caractéristiques de $f$, on cherche tout d'abord l'ensemble $D$ des vecteurs changés en leur opposé. Ce sera l'axe de $g$. On oriente $D$ par un vecteur $e$. Le plan de $h$ est alors l'orthogonal de $e$ et $h$ est complètement déterminée. Comme $f = g\\circ h$, on a aussi $g = f\\circ h$ et la détermination de $g$ est alors facile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyse_reflex_rot(A):\n",
    "    e = (A + eye(3)).nullspace()[0]\n",
    "    B = matrice_reflexion(e) * A\n",
    "    (e1, theta) = analyse_rotation(B)\n",
    "    return (e, theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons l'une des matrices de rotation vues ci-dessus et prenons l'opposé de sa première colonne."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A4 = Matrix([[-3, 1, sqrt(6)],[-1, 3, -sqrt(6)],[sqrt(6),sqrt(6), 2]])/4\n",
    "A4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e, theta = analyse_reflex_rot(A4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vérifications ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrice_reflexion(e) * matrice_rotation(e, theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "matrice_rotation(e, theta) * matrice_reflexion(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On retrouve bien la matrice $A_4$, ceci quel que soit le sens dans lequel on effectue les produits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse_rotation(Matrix([[1,0,0],[0,1,0],[0,0,1]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquons pour terminer que $-id$ se décompose aussi comme un produit de réflexion et de rotation, mais on n'a plus unicité."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse_reflex_rot(-eye(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 4 Regroupons le tout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyse(A):\n",
    "    if simplify(A * A.T) != eye(3):\n",
    "        return (0, [])\n",
    "    elif est_reflexion(A):\n",
    "        return (1, analyse_reflexion(A))\n",
    "    elif est_rotation(A):\n",
    "        return (2, analyse_rotation(A))\n",
    "    else:\n",
    "        return (3, analyse_reflex_rot(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse(A4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse(A0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse(2*eye(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}