NumPy-Algebre_lineaire.ipynb

{
 "metadata": {
  "name": "NumPy - Alg\u00e8bre lin\u00e9aire"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Alg\u00e8bre lin\u00e9aire avec NumPy."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "On commence par importer tout ce qui est contenu dans le NumPy, et par d\u00e9finir une matrice m inversible."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = array( [[1, 2, 3],\n",
      "            [2, 3, 4],\n",
      "            [3, 4, 6]] )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Les op\u00e9rations classiques li\u00e9es \u00e0 la r\u00e9solution de syst\u00e8mes lin\u00e9aires et \u00e0 la r\u00e9duction d'endomorphimes sont quant \u00e0 elles disponibles dans la biblioth\u00e8qe linalg."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "linalg.det(m) # d\u00e9terminant"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "-1.0000000000000004"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "minv = linalg.inv(m) # inverse\n",
      "print minv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[-2.  0.  1.]\n",
        " [ 0.  3. -2.]\n",
        " [ 1. -2.  1.]]\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "On peut v\u00e9rifier que minv est bien l'inverse (num\u00e9rique) de m."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print dot(minv, m) # M^{-1} . M = I"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  1.00000000e+00  -1.33226763e-15  -1.77635684e-15]\n",
        " [  0.00000000e+00   1.00000000e+00   0.00000000e+00]\n",
        " [  0.00000000e+00   4.44089210e-16   1.00000000e+00]]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Pour r\u00e9soudre une syst\u00e8me lin\u00e9aire, on utilise la fonction linalg.solve(). Notez que la matrice doit \u00eatre carr\u00e9e et inversible."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "b = array([2, 2, 3])\n",
      "x = linalg.solve(m, b) # r\u00e9solution de syst\u00e8me\n",
      "print x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[-1.  0.  1.]\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print dot(m, x) # on v\u00e9rifie que la solution est correcte"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 2.  2.  3.]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "On peut aussi calculer les valeurs propres d'une matrice, et les vecteurs propres associ\u00e9s."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "valp, vecp = linalg.eig(m) # valeurs et vecteurs propres\n",
      "print \"Valeurs propres de m :\", valp\n",
      "print \"Vecteurs propres associ\u00e9s :\"\n",
      "print vecp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Valeurs propres de m : [ 10.18669767  -0.42027581   0.23357814]\n",
        "Vecteurs propres associ\u00e9s :\n",
        "[[ 0.36530658  0.92694664  0.08556306]\n",
        " [ 0.52826911 -0.13074871 -0.83894966]\n",
        " [ 0.7664743  -0.35167415  0.53744064]]\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print dot(m, vecp[:,0])   # m * 1er vecteur propre\n",
      "print valp[0] * vecp[:,0] # 1ere valeur propre * 1er vecteur propre"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 3.72126769  5.38131768  7.80784197]\n",
        "[ 3.72126769  5.38131768  7.80784197]\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "P = vecp\n",
      "D = diag(valp)\n",
      "Pinv = linalg.inv(P)\n",
      "print dot(P, dot(D, Pinv)) # principe de diagonalisation : m = P . D . Pinv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 1.  2.  3.]\n",
        " [ 2.  3.  4.]\n",
        " [ 3.  4.  6.]]\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Illustration d'un probl\u00e8me d'instabilit\u00e9 num\u00e9rique."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "On consid\u00e8re le syst\u00e8me lin\u00e9aire A.x = b d\u00e9fini comme suit."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = array( [[10., 7., 8., 7.],\n",
      "            [7., 5., 6., 5.],\n",
      "            [8., 6., 10., 9.],\n",
      "            [7., 5., 9., 10.]] )\n",
      "b = array( [32, 23, 33, 31 ] )\n",
      "print \"A = \", A\n",
      "print \"b = \", b.T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "A =  [[ 10.   7.   8.   7.]\n",
        " [  7.   5.   6.   5.]\n",
        " [  8.   6.  10.   9.]\n",
        " [  7.   5.   9.  10.]]\n",
        "b =  [32 23 33 31]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Cherchons une solution x \u00e0 ce probl\u00e8me."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = linalg.solve(A, b)\n",
      "print x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1.  1.  1.  1.]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Nous trouvons bien la solution attendue."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Appliquons maintenant une l\u00e9g\u00e8re pertubation aux coefficients de A, puis essayons de r\u00e9soudre le nouveau syst\u00e8me lin\u00e9aire ainsi obtenu."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "E = array( [[0., 0., 0.1, 0.2],\n",
      "            [0.08, 0.04, 0., 0.],\n",
      "            [0., -0.02, -0.02, 0.],\n",
      "            [-0.01, -0.01, 0., -0.02]] )\n",
      "A2 = A + E\n",
      "x2 = linalg.solve(A2, b)\n",
      "print x2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ -5.79224365  12.01880932  -1.56555429   2.56552237]\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "On constate qu'une petite pertubation sur les coefficients de A provoque un changement assez important sur le r\u00e9sultat du syst\u00e8me lin\u00e9aire."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "On peut quantifier le changement sur la matrice en calculant la norme 2 de E = A - A2."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def matrix_norm2(M): # calcule la norme 2 de la matrice M\n",
      "    valp, _ = linalg.eig( dot(M.conj(), M) )\n",
      "    return sqrt( max(valp) )\n",
      "\n",
      "matrix_norm2(E)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "(0.050463861837858059+0j)"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "De m\u00eame, on peut quantifier le changement au niveau de la solution du syst\u00e8me, en calculant la norme 2 de x-x2."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def vect_norm2(v): # calcule la norme 2 du vecteur v\n",
      "    return sqrt( sum(v**2) )\n",
      "\n",
      "vect_norm2( x-x2 )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "13.288403278162946"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Pour justifier l'\u00e9cart important entre les vecteurs x et x2, il faut regarder le conditionnement de la matrice A. En effet, en d\u00e9finissant cond_2(A) = |||A|||_2 . |||A^{-1}|||_2, on peut montrer que :\n",
      "    || x - x2 ||_2 <= cond_2(A) . ||| A - A2 |||_2.\n",
      "Ainsi, plus le conditionnement de A est \u00e9lev\u00e9, et plus on a un risque d'instabilit\u00e9 num\u00e9rique."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cond = matrix_norm2(A) * matrix_norm2( linalg.inv(A) )\n",
      "cond"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "2984.0927016755245"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "On constate que le condition de A est tr\u00e8s \u00e9lev\u00e9. On en d\u00e9duit une borne sur || x - x2 ||_2 qui est tr\u00e8s \u00e9lev\u00e9e, qui laisse largement place \u00e0 une instabilit\u00e9 num\u00e9rique."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cond * matrix_norm2(E)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "(150.58884180871425+0j)"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}