{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On cherche à fabriquer un toboggan partant du point $A=(x_A,y_A)$ et arrivant\n",
    "au point $B=(x_B,y_B)$ de telle sorte que le temps mis par une bille qui part\n",
    "sans vitesse du point $A$ arrive le plus vite possible au point $B$. On suppose\n",
    "que cette bille n'est soumise qu'à la gravitation et à la réaction\n",
    "du toboggan, sans frottement. Si la forme du toboggan est donnée par la courbe d'équation $y=f(x)$,\n",
    "le temps mis par la boule pour aller de $A$ vers $B$ est\n",
    "\n",
    "$$\n",
    "T(f)\\,=\\, \\int_{x_A}^{x_B} \\sqrt{\\frac{1+(f'(x))^2}{2\\,g\\,(y_A-f(x))}} \\,dx\n",
    "$$\n",
    "\n",
    "o\\`u $g=9,8 m\\,s^{-2}$ est l'accélération de la pesanteur. Dans la suite on suppose $x_A=0$ et $B=(1,0)$, \n",
    "on cherche donc une solution du problème\n",
    "\n",
    "$$\n",
    "(P) \\quad \\quad \\inf\\left\\{ T(f):\\,\\,\\, f: [0,1] \\to \\mathbb{R}, \\,f(0)=y_A, \\,f(1)=0 \\right\\}.\n",
    "$$\n",
    "\n",
    "Afin de résoudre ce problème de manière approchée, on le discrétise par\n",
    "la méthode des éléments finis.\n",
    "Pour $N \\geq 1$ fixé, on considère l'ensemble $\\mathcal{A}_N$ des fonctions\n",
    "affines par morceaux $f: [0,1] \\to \\mathbb{R}$ telles que $f(0)=y_A$, $f(1)=0$ et\n",
    "\n",
    "$$\n",
    "\\forall n \\in \\{0,\\ldots,N\\}, \\quad \\quad\n",
    "\\mbox{$f$ affine sur }\\left[\\frac{n}{N+1},\\frac{n+1}{N+1}\\right].\n",
    "$$\n",
    "\n",
    "Toute fonction $f \\in \\mathcal{A}_N$ est alors caractérisée par le vecteur $y \\in \\mathbb{R}^{N+2}$\n",
    "donné par\n",
    "\n",
    "$$\n",
    "\\forall n \\in \\{0,\\ldots,N+1\\}, \\quad \\quad y_n = f\\left(\\frac{n}{N+1}\\right).\n",
    "$$\n",
    "\n",
    "L'ensemble des vecteurs ainsi obtenus est alors\n",
    "\n",
    "$$\n",
    "\\mathcal{B}_N = \\{ y \\in \\mathbb{R}^{N+2}\\,:\\, y_0 = y_A, \\, y_{N+1} = 0\\}. \n",
    "$$\n",
    "\n",
    "On doit alors résoudre le problème :\n",
    "\n",
    "$$\n",
    "(P_N) \\quad \\quad \\inf\\{ J(y) \\,:\\, y \\in \\mathcal{B}_N \\}\n",
    "$$\n",
    "\n",
    "où $J(y) = T(f)$ pour la fonction affine $f \\in \\mathcal{A}_N$ associée au vecteur $y$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour commencer, on va donc déclarer les variables du problème discrétisé sous Octave, et afficher les valeurs pour l'utilisateur. On définit donc $N$, $y_A$ et $y_B$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nombre de points de discretisation : N = 9 \n",
      "valeurs de y_A et y_B : y_A = 0.200000 , y_B = 0.000000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "N=9;\n",
    "printf(\"nombre de points de discretisation : N = %d \\n\",N);\n",
    "yA=0.2;\n",
    "yB=0;\n",
    "printf(\"valeurs de y_A et y_B : y_A = %f , y_B = %f\\n\\n\",yA, yB);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On aura aussi besoin d'un vecteur initial pour la méthode de la plus grande pente, on prend l'interpolation affine entre les valeurs $y_A$ et $y_B$ sur les $N+2$ points $\\frac{i}{N+1}$ pour $i \\in \\{0,\\ldots,N+1\\}$. Ainsi le vecteur $y$ est de taille $N+2$, avec $y_1 =y_A$ et $y_{N+2}=y_B$ fixés, et les autres points sont obtenus par interpolation : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y =\n",
      "\n",
      "   0.20000\n",
      "   0.18000\n",
      "   0.16000\n",
      "   0.14000\n",
      "   0.12000\n",
      "   0.10000\n",
      "   0.08000\n",
      "   0.06000\n",
      "   0.04000\n",
      "   0.02000\n",
      "   0.00000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "y=zeros(N+2,1);\n",
    "i=1;\n",
    "while (i <= N+2)\n",
    "  y(i) = ((N+2-i)*yA+(i-1)*yB)/(N+1);\n",
    "  i=i+1;\n",
    "endwhile\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Désormais on définit les paramètres $\\rho$ (le pas) et $\\varepsilon$ (la précision) pour la méthode de la plus grande pente :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho=10^(-2);\n",
    "epsilon=10^(-2);\n",
    "eps = epsilon^2;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On va maintenant déclarer une fonction funcJ pour calculer le temps $T_N(y)$ associé au toboggan correspondant à la fonction affine dont les valeurs aux points de discrétisation sont contenues dans $y$. Pour cela, on s'appuie sur la formule\n",
    "\n",
    "$$\n",
    "J(y) \\,\\, = \\,\\, T(f)\\,\\,=\\,\\, \\sum_{n=0}^N \\int_{\\frac{n}{N+1}}^{\\frac{n+1}{N+1}}\n",
    "\\sqrt{\\frac{1+(f'(x))^2}{y_A-f(x)}} \\,dx,\n",
    "$$\n",
    "\n",
    "où on oublie volontairement la constante $2g$.\n",
    "Chaque terme de cette somme\n",
    "ne dépend que de $y_n$ et $y_{n+1}$ et peut être calculé par la formule :\n",
    "\n",
    "\\begin{align*}\n",
    "\\int_{\\frac{n}{N+1}}^{\\frac{n+1}{N+1}} \\sqrt{\\frac{1+(f'(x))^2}{(y_A-f(x))}} \\,dx\n",
    "&= \\frac{2\\,\\sqrt{1+((N+1)\\,(y_{n+1}-y_n))^2}}{(N+1)\\,(y_{n+1}-y_n)}\\left(\\sqrt{y_A-y_n} - \\sqrt{y_A-y_{N+1}} \\right)\\\\\n",
    "&= \\frac{2\\,\\sqrt{1+((N+1)\\,(y_{n+1}-y_n))^2}}{(N+1)\\,(\\sqrt{y_A-y_n} + \\sqrt{y_A-y_{N+1}})}\\,.\n",
    "\\end{align*}\n",
    "\n",
    "Notons que cette dernière formule est bien définie même dans le cas $y_n = y_{n+1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On obtient donc la fonction suivante en sommant ces $N$ termes :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "function r=funcJ(x)\n",
    "  r=0.;\n",
    "  n=length(x);\n",
    "  yA=x(1);\n",
    "  i=1;\n",
    "  while (i <= n-1)\n",
    "    a = x(i);\n",
    "    b = x(i+1);\n",
    "    r = r + 2*sqrt(1+((b-a)*(n-1))**2)/((n-1)*(sqrt(yA-a) + sqrt(yA-b)));\n",
    "    i = i+1;\n",
    "  endwhile\n",
    "endfunction\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On utilise ensuite une fonction gradJ pour la différenciation de $J$, dans laquelle on fait attention à ce que les dérivées partielles pour la première et la dernière coordonnée soient nulles : il ne faut pas faire varier les valeurs $y_A$ et $y_B$ dans le vecteur y !\n",
    "\n",
    "Les dérivées partielles sont calculées grâce à la formule de différence finie :\n",
    "\n",
    "$$\n",
    "\\frac{\\partial J}{\\partial x_i} \\,=\\, \\frac{J(x+h \\, e_i) - J(x-h \\,e_i)}{2\\,}\n",
    "$$\n",
    "\n",
    "où $e_i$ est le $i$-ème vecteur de la base canonique de $\\mathbb{R}^{N+2}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "function z=gradJ(x)\n",
    "  h=10^(-8);\n",
    "  n=length(x); \n",
    "  i=2;\n",
    "  z=zeros(n,1);\n",
    "  while (i <= (n-1))\n",
    "    vecbase = zeros(n,1);\n",
    "    vecbase(i)=1;\n",
    "    z(i)=(funcJ(x+h*vecbase)-funcJ(x-h*vecbase))/(2*h);\n",
    "    i=i+1;\n",
    "  endwhile\n",
    "endfunction\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enfin on applique la méthode de la plus grande pente et on affiche les résultats : le vecteur $y$ et le nombre d'itérations $k$ :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solution calculéee :\n",
      "y =\n",
      "\n",
      "   0.20000\n",
      "   0.00008\n",
      "  -0.07206\n",
      "  -0.11881\n",
      "  -0.14759\n",
      "  -0.16143\n",
      "  -0.16148\n",
      "  -0.14774\n",
      "  -0.11901\n",
      "  -0.07227\n",
      "   0.00000\n",
      "\n",
      "nombre d'iteration k = 211\n"
     ]
    }
   ],
   "source": [
    "k=0; #compteur\n",
    "while (gradJ(y)'*gradJ(y) > eps)\n",
    "  y=y-rho*gradJ(y);\n",
    "  k=k+1;\n",
    "  #printf(\"test k=%d : %f\\n\", k,sqrt(gradJ(y)'*gradJ(y))) #Afficher pour chaque iteration la valeur du test\n",
    "  endwhile\n",
    "printf(\"solution calculéee :\\n\");\n",
    "y\n",
    "printf(\"nombre d'iteration k = %d\\n\", k);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et pour finir, on trace le toboggan obtenu. Pour cela, il faut définir le vecteur $x$ des abscisses, puis tracer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAJNmlDQ1BkZWZhdWx0X3JnYi5pY2MAAHiclZFnUJSHFobP933bCwvssnRYepMqZQHpvUmvogJL7yxLEbEhYgQiiog0RZCggAGjUiRWRLEQFBSxoFkkCCgxGEVUUPLDOxPn3vHHfX49884755yZA0ARBQBARQFSUgV8Pxd7TkhoGAe+IZKXmW7n4+MJ3+X9KCAAAPdWfb/zXSjRMZk8AFgGgHxeOl8AgOQCgGaOIF0AgBwFAFZUUroAADkLACx+SGgYAHIDAFhxX30cAFhRX30eAFj8AD8HABQHQKLFfeNR3/h/9gIAKNvxBQmxMbkc/7RYQU4kP4aT6ediz3FzcOD48NNiE5Jjvjn4/yp/B0FMrgAAwCEtfRM/IS5ewPmfoUYGhobw7y/e+gICAAh78L//AwDf9NIaAbgLANi+f7OoaoDuXQBSj//NVI8CMAoBuu7wsvjZXzMcAAAeKMAAFkiDAqiAJuiCEZiBJdiCE7iDNwRAKGwAHsRDCvAhB/JhBxRBCeyDg1AD9dAELdAOp6EbzsMVuA634S6MwhMQwhS8gnl4D0sIghAROsJEpBFFRA3RQYwQLmKNOCGeiB8SikQgcUgqkoXkIzuREqQcqUEakBbkF+QccgW5iQwjj5AJZBb5G/mEYigNZaHyqDqqj3JRO9QDDUDXo3FoBpqHFqJ70Sq0ET2JdqFX0NvoKCpEX6ELGGBUjI0pYboYF3PAvLEwLBbjY1uxYqwSa8TasV5sALuHCbE57COOgGPiODhdnCXOFReI4+EycFtxpbga3AlcF64fdw83gZvHfcHT8XJ4HbwF3g0fgo/D5+CL8JX4Znwn/hp+FD+Ff08gENgEDYIZwZUQSkgkbCaUEg4TOgiXCcOEScICkUiUJuoQrYjexEiigFhErCaeJF4ijhCniB9IVJIiyYjkTAojpZIKSJWkVtJF0ghpmrREFiWrkS3I3uRo8iZyGbmJ3Eu+Q54iL1HEKBoUK0oAJZGyg1JFaadco4xT3lKpVGWqOdWXmkDdTq2inqLeoE5QP9LEado0B1o4LYu2l3acdpn2iPaWTqer023pYXQBfS+9hX6V/oz+QYQpoifiJhItsk2kVqRLZETkNYPMUGPYMTYw8hiVjDOMO4w5UbKouqiDaKToVtFa0XOiY6ILYkwxQzFvsRSxUrFWsZtiM+JEcXVxJ/Fo8ULxY+JXxSeZGFOF6cDkMXcym5jXmFMsAkuD5cZKZJWwfmYNseYlxCWMJYIkciVqJS5ICNkYW53txk5ml7FPsx+wP0nKS9pJxkjukWyXHJFclJKVspWKkSqW6pAalfokzZF2kk6S3i/dLf1UBiejLeMrkyNzROaazJwsS9ZSlidbLHta9rEcKqct5ye3We6Y3KDcgryCvIt8uny1/FX5OQW2gq1CokKFwkWFWUWmorVigmKF4iXFlxwJjh0nmVPF6efMK8kpuSplKTUoDSktKWsoByoXKHcoP1WhqHBVYlUqVPpU5lUVVb1U81XbVB+rkdW4avFqh9QG1BbVNdSD1Xerd6vPaEhpuGnkabRpjGvSNW00MzQbNe9rEbS4Wklah7XuaqPaJtrx2rXad3RQHVOdBJ3DOsOr8KvMV6Wualw1pkvTtdPN1m3TndBj63nqFeh1673WV9UP09+vP6D/xcDEINmgyeCJobihu2GBYa/h30baRjyjWqP7q+mrnVdvW92z+o2xjnGM8RHjhyZMEy+T3SZ9Jp9NzUz5pu2ms2aqZhFmdWZjXBbXh1vKvWGON7c332Z+3vyjhamFwOK0xV+WupZJlq2WM2s01sSsaVozaaVsFWnVYCW05lhHWB+1Ftoo2UTaNNo8t1WxjbZttp2207JLtDtp99rewJ5v32m/6GDhsMXhsiPm6OJY7DjkJO4U6FTj9MxZ2TnOuc153sXEZbPLZVe8q4frftcxN3k3nluL27y7mfsW934Pmoe/R43Hc09tT75nrxfq5e51wGt8rdra1LXd3uDt5n3A+6mPhk+Gz6++BF8f31rfF36Gfvl+A/5M/43+rf7vA+wDygKeBGoGZgX2BTGCwoNaghaDHYPLg4Uh+iFbQm6HyoQmhPaEEcOCwprDFtY5rTu4bircJLwo/MF6jfW5629ukNmQvOHCRsbGyI1nIvARwRGtEcuR3pGNkQtRblF1UfM8B94h3qto2+iK6NkYq5jymOlYq9jy2Jk4q7gDcbPxNvGV8XMJDgk1CW8SXRPrExeTvJOOJ60kByd3pJBSIlLOpYqnJqX2pymk5aYNp+ukF6ULMywyDmbM8z34zZlI5vrMHgFLkC4YzNLM2pU1kW2dXZv9ISco50yuWG5q7uAm7U17Nk3nOef9tBm3mbe5L18pf0f+xBa7LQ1bka1RW/u2qWwr3Da13WX7iR2UHUk7fiswKCgveLczeGdvoXzh9sLJXS672opEivhFY7std9f/gPsh4YehPav3VO/5UhxdfKvEoKSyZLmUV3rrR8Mfq35c2Ru7d6jMtOzIPsK+1H0P9tvsP1EuVp5XPnnA60BXBaeiuOLdwY0Hb1YaV9YfohzKOiSs8qzqqVat3le9XBNfM1prX9tRJ1e3p27xcPThkSO2R9rr5etL6j8dTTj6sMGloatRvbHyGOFY9rEXTUFNAz9xf2pplmkuaf58PPW48ITfif4Ws5aWVrnWsja0Latt9mT4ybs/O/7c067b3tDB7ig5BaeyTr38JeKXB6c9Tved4Z5pP6t2tq6T2VnchXRt6prvju8W9oT2DJ9zP9fXa9nb+aver8fPK52vvSBxoewi5WLhxZVLeZcWLqdfnrsSd2Wyb2Pfk6shV+/3+/YPXfO4duO68/WrA3YDl25Y3Th/0+LmuVvcW923TW93DZoMdv5m8lvnkOlQ1x2zOz13ze/2Dq8ZvjhiM3LlnuO96/fd7t8eXTs6/CDwwcOx8DHhw+iHM4+SH715nP146cn2cfx48VPRp5XP5J41/q71e4fQVHhhwnFi8Ln/8yeTvMlXf2T+sTxV+IL+onJacbplxmjm/Kzz7N2X615OvUp/tTRX9KfYn3WvNV+f/cv2r8H5kPmpN/w3K3+XvpV+e/yd8bu+BZ+FZ+9T3i8tFn+Q/nDiI/fjwKfgT9NLOcvE5arPWp97v3h8GV9JWVn5By6ikLxSF1/9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAAHXRFWHRTb2Z0d2FyZQBHUEwgR2hvc3RzY3JpcHQgOS4xOJQFEHMAABlcSURBVHic7d09buxIuibgqEG5Y7CtOi4bdwECyxrgYhxugY27AmkFhdQK7ijRK7i5BNFuq2m0O0ATuYKiW+URg3HG0xihw8rOo9TJfwbJ57FSqvzhYSnzzYjvi+APb29vAQDG9t/GPgAACEEgAZAIgQRAEgQSAEkQSAAkQSABkASBBEASBBIASfjx8qdo27bv+xBClmVFUXx4n7qu4x0+uQ8AS3aFEVLXdWVZlmXZdd2h+2RZVlVVWZYxugBgzxUCKcuyvRvfKssyhND3fdM0l78iAPNzhSm7I/V9v16vV6vVt/+p60JV/a//+I/3g/ny5ctPP/10twOblu12+/DwMPZRTIATdTzn6khO1Cd+//333377Ld7+xz/+8be//e2MJ7lCIA2zcLvTcUNVKf7Ytm1d1y8vLx8+Q56Hn376H7/88j8vP5gliMNNvsuJOp5zdSQn6qauMGWX53nTNE3T5Hk+/PLp6enp6Snebtt2vV6XZdk0zXq9vvwVAZifK4yQPuyae3193b3D8OOh7xd/+tP/6bqwk2gALEsq65D+/d//n3aHY+RC+zhO1PGcqyM5UUf68uXLeQ9MJZDy/E8awo/hLXEkJ+p4ztWRnKgjnd2VlkogAbBwCQVSnofDK2sBmLmEAqksgzISwGIlFEhZFpSRABYroUACYMnSCiRlJIDFSiuQqirU9dgHAcAY0gokABZLIAGQhOQCqShC2459EADcXXKBZDUSwDIlF0gALJNAAiAJKQaSMhLAAqUYSMpIAAuUYiCFELJs7CMA4L4SDSS7rAIsTaKBZNYOYGkSDSR9DQBLk2ggAbA0AgmAJKQbSMpIAIuSbiApIwEsSrqBBMCiCCQAkpB0IJWlK5oDLEXSgVQUoevGPggA7iLpQAJgOVIPpCyzrx3AIqQeSFYjASxE6oGU58pIAIuQeiABsBATCCRlJIAlmEAgKSMBLMEEAkkZCWAJJhBIACzBNALJIAlg9qYRSMpIALM3jUDSaAcwe9MIJABm78dTH9C2bd/3IYQsy4qi+PA+dV23bVuWZVmW8TfN1xm3Tx71uVhGyvMzHgrABJw8Quq6LiZNd7jNoKqql5eXvV/GR52XRkEZCWDuTh4hZVm2d+PIR202mzzPD42QttvtcDvP8/yboZAyEkCyuq4bRinb7XaYHjvJyYF0nqIoYg41B4Y5Dw8P5/0DABjdhwOJU508Zdd/Haf0OwOWvu/7249frEYCmLGTR0h5nsdRzm4YPj09hRBeX1/jj7HxoW3bYYJuvV7H+bpLcqssQ12Hx8eznwCAdJ0cSB9WgIYo2r3P7hTcarU6/dj2KSMBzJh1SAAkYWKBVBShbcc+CABuYGKBZDUSwFxNLJAAmCuBBEASphdIykgAszS9QFJGApil6QUSALM0yUA6ZVtXAKZhkoFUFGbtAOZmqoGkrwFgZiYZSADMj0ACIAlTDSTN3wAzM9VAUkYCmJmpBhIAMyOQAEjChANJGQlgTiYcSMpIAHMy4UACYE6mHUhZFvp+7IMA4BqmHUjKSACzMe1AyvPQdWMfBADXMO1AAmA2Jh9IykgA8zD5QFJGApiHyQeSMhLAPEw+kACYhzkEkkESwAzMIZCUkQBmYA6BpNEOYAbmEEgAzMBMAkkZCWDqZhJIykgAUzeTQFJGApi6mQQSAFM3n0BSRgKYtPkEkjISwKTNJ5CUkQAmbT6BBMCkzSqQiiK07dgHAcBZfrz8Kdq27fs+hJBlWVEUH96nruu2bcuyLMvy8lc8pCzDeh0OHAIASbvCCKnrupg03eEut6qqXl5eLn8tAObqCoGUZdneDQA41RWm7K5iu90Ot/M8z/P8vOeJZSSzdgD31HXdMEm23W7Pq85cIZD6r93W/U7b9VBVOvJJHh4erlJeUkYCuL9LBhKDKwRSnudN08Qbwy+fnp5CCK+vr/HH2PjQtu0njQ8ALNkVAunDgBmiaPc+N22xA2DSZrUOKbKHEMAUzTCQLI8FmKIZBhIAUySQAEjCPANJGQlgcuYZSMpIAJMzz0ACYHIEEgBJmG0gKSMBTMtsA0kZCWBaZhtIAEzLnAMpy8LO/uMAJG3OgaSMBDAhcw6kPA+HL6oOQFrmHEgATMjMA0kZCWAqZh5IykgAUzHzQFJGApiKmQcSAFMx/0BSRgKYhPkHUlWFuh77IAD4nvkHkhESwCTMP5AAmIRFBJJeO4D0LSKQrEYCSN8iAkkZCSB9iwgkANK3lEBSRgJI3FICSRkJIHFLCSRlJIDELSWQAEjcggKpKELbjn0QABywoEBSRgJI2YICCYCUCSQAkrCsQFJGAkjWsgJJGQkgWcsKJACSJZAASMLiAqkslZEAUrS4QCoKZSSAFC0ukABI04+nPqBt277vQwhZlhVFceR9mq+jkk8eBcCSnTxC6rquLMuyLLvD1xf68D7xNymkkeZvgASdPELKsmzvxjH3ybJss9nkeX5ohLTdbofbeZ7neX7qgR2vKMJ6Hcrydq8AsCxd1w0jkO12W571CXtyIJ2nKIq9ubs9Dw8P5/0DABjdVQYSJ0/Z9V+vc9fvXPCu7/u9H7+9DwB84uQRUp7ncZSzG4ZPT08hhNfX10P3Wa/Xcb4ukYiKZSRDMoB0/PD29jb2MYQQQtM0d56yW6/DanXPFwRYhLM/z61DAiAJyw2kLAtpTB8CEMKSA8lqJICkLDeQ8jwcXtoLwL0tN5AASMqiA0kZCSAdiw4kZSSAdCw6kJSRANKx6EACIB1LDyRlJIBELD2QqirU9dgHAYBAMkICSMTSAwmARAgkvXYASRBIViMBJEEgKSMBJEEgAZAEgRSCMhJAAgRSCMpIAAkQSCEoIwEkQCC9E0gA4xJI78oytO3YBwGwYALpnTISwLgEEgBJEEgAJEEg/aEolJEARiOQ/qCMBDAigQRAEgQSAEkQSP9CGQlgLALpX1geCzAWgbTPHkIAoxBIACRBIO3T/A0wCoG0T18DwCgEEgBJEEgAJEEgfUAZCeD+BNIHlJEA7k8gAZAEgfSxLLNCFuCufrzPy7Rt2/d9CCHLsqIo7vOil4hlpKoa+zgAFuNOI6Su68qyLMuy67r7vOKF8jxM5EgBZuJOgZRl2d4NANh1pym779put8PtPM/zPB/xYKKyDHVt1g7g+7quGybA/vM//3tZnvMkdwqk/muHQH+gVeDh4aE8719wM7H5u+tCAuEIkLRhIFHXoar+93lPcqcpuzzPm6ZpmiaFoc/xHh/DZjP2QQBMRNeFvg//9m//97yH32mENInOug+tVmG9DqvV2McBkLzNJry8nL/TjXVI35FlIc9t3ADwHZd/dxdI31dVoa6tkwU4qG1DnocL26gF0lFWK8UkgI/1/XV2EhBIR8my9y5wAPZsNuHx8QrPI5COVRTvDSQADOo6lOWlk3WRQDpB7LgDIOq60HXhWm3UAuk0ViYBDDaba66KEUiniW0kusABrlU6Ggikk1VVaBrFJGDR2vZ9meYVCaRzmLgDluxafd57BNI5siwUhS5wYKGuPlkXCaQzlaUucGCJrtjnvUcgnU8XOLA0fX/NPu89AukiiknAotz06gcC6SK6wIHluFHpaCCQLhX3AgeYt1v0ee8RSFegmATMXl1fv897j0C6gtgFfvZFEgESd58LZwuk6yjL0La6wIEZappQFDfp894jkK7GxB0wP30f2jaU5T1eSyBdU1XpAgdm5T6TdZFAuqa4WKzrxj4OgGu4QyPDLoF0ZZbKAvMQd0e70aYMHxJI16eYBMzArZfBfksgXZ8ucGDq7lk6Ggikm9AFDkxX07zvi3ZnAulWTNwBUxT7vO/ZyzAQSDdkmztgcu5fOhoIpBsqivdrhwBMwu0uvncMgXRbusCBqbh/n/cegXRziknAJIw4WRcJpJuLVxBxET8gZaP0ee8RSPdQVaFpdIEDiWrbcfq89wikO1FMAtLU96Fpxunz3iOQ7iTLQlnqAgeSM3rpaCCQ7kcXOJCacfu89wikuzJxB6QjfkUesc97j0C6N13gQCJS6KzbJZDuTRc4kIJ0SkcDgTSCuMedLnBgLPE7cZ6PfRz/SiCNY7VSTAJGU9fJDY+CQBqLLnBgLKmVjgY/nvqAtm37vg8hZFlWHGjO+PY+zdfrp37yqKWJV5Xt+1QaLoElaJpQFIl+7Jw8Quq6rizLsiy7wwtqPrxP/I002qXjDrinePG9shz7OA44eYSUfQ3W7HDCfnufLMs2m02e54dGSNvtdrid53meWq3tZuLKpAQnc4H5Wa/Dy8tNnrnrumEEst1uy7NC77NA6vu+3WlPLorikxD6XFEUe3N3ex4eHs77B0xd3NCwbRNamwbM0k2/+15lIPFZIGVZ9m1I9F+7lfudtuWhYvTJfTikqsJ6ncRWu8BcxdFL4nNPJ0/Z5XkeRzm7Yfj09BRCeH19PXSf9Xod5+tE1IfixF2afS/ADGw2t5qsu6If3t7exj6GEEJommaZU3aD2HGXwg7wwMys1+Hx8X5zMGd/nluHlIqyfL+gPcAVpdznvUcgJUQXOHBdifd57xFIaXF9CuCKkt2U4UMCKS1DFzjAhep6YmVpgZScuBc4wCViTXpaCxwFUooUk4ALTXELGIGUoix733oV4AzTKh0NBFKiyjK0rS5w4GRtO9WdXwRSukzcAafq+9A0E+tlGAikpFWVLnDgBHFThokSSEmLHTKHrzwF8IfY5z3FybpIIKXOUlngGH0fum5ifd57BNIEKCYB3zXRzrpdAmkCdIEDn5viqqNvCaRp0AUOHNK2IctSv/jeMQTSZJi4A7416T7vPQJpSmxzB+yZx2RdJJCmJF5lSyYB0WYTynLCfd57BNLElGXIc5kEhPU6FMW0+7z3CKTpKYqQ5+H5eezjAMYTd2SYUxoFgTRRRaHHARaq78Pzc3h8nM9M3UAgTVWWhcfHsF7rBYcF6fuw2YTVaoZpFATSpMkkWJR5p1EQSFOXZeHlJWw2Mglmrm3nsDnQ5wTSHMR6UtuOfRzAbbRt6Lrw8jL2cdyYQJqJl5fQtjIJZqiuQ9fNZC+Gzwmk+Xh8DF0nk2BW4qLDJaRREEgzU1WhbS2bhZmo65DnS0mjIJDm5/HRVg4wB8/PIc/ntvT1cwJphuJWDpbNwnTFhrpFpVEQSHNVFKGqZBJMT9+/bws018VGnxBIs5Xn78tmgalYchoFgTRvcSuHpyfLZmEC4kYMLy8LTaMgkGYvy8J//ZetHCB1S9iI4bsE0iKsVmGzCV039nEAH4mr2me/EcN3CaSlWK1CXVs2C8mJ2wLN5jLklxBIC7Ja2coB0hK/Ji5n6evnBNKyVFXourDZjH0cwNeNGIyNBgJpcarKVg4wvvV6cRsxfJdAWqKylEkwpvU6lKU02ieQFipuL/T8PPZxwPLEpa/S6Fs/Xv4Ubdv2fR9CyLKsOHCO67pu27Ysy7IsL39FrmLY8m7hSx/gbuJGDDO+BvmFrjBC6rouJk13eJ1LVVUveuzTE7dyeH62bBZuzkYM33WFQMq+nt3MaZ6gLHu/ArpMgtuJ3a1mIz532pRd3/ftzjKWoiiuFULb7Xa4ned5nudXeVqOkWXh5WXRWzrCTcXLZs57kqjrumGSbLvdnledOS2Qsiz79mX6r1+t+53v2ENV6chnfnh4UF4a12oVnp9DVam1wjXFjRjmnUbhSgOJKzQ15HneNE28Mfzy6ekphPD6+hp/jI0Pbdt+0vjA6F5e3nvB/S+Cq9hsQpbZiOFYP7y9vY19DCGE0DSNEVIi6jpkWfB/Ay602LfS2Z/n1iGxr6pC31s2CxeJ2wItMI0uIZD4gO2F4BLPz7YFOodA4mO2coDzxKWv0ugMAomDiiI8Pob1euzjgIno+/D8bPnE+QQSn4l748sk+C7bAl1OIPEdtheC77It0FUIJL4vbuWw2cgk+EDb2qT4OgQSx1qtwmbjCujwLxayEcN9CCROsFqFppFJ8K5pQtfZiOFqBBKnWa1C18kkCHUtja5MIHGyqnrfSx8WK27E8Pg49nHMi0DiHHFTcFs5sEzrtY0YbkIgcaa4lYNxEkuzXrtKy60IJM5XFKEobC/EgsTrWLqA6I1c4XpILNmw5V1Z2tiYOWua0DQ2YrgtIyQuFZfNhhCen8PXSxjDfHTd+zSAjRhuzQiJ64gjpLoOdW1zSWYibgiU59a93olA4prixf3ie9j6DCYtLjPy7eqeBBJXlmXvi2cVlpioWC56fPSl6t7UkLiJYZZDYYkJ2S0XaaW7P4HEDZVleHkJbevqFaQuXluvbcPLi2H9aEzZcXNVFarq/Sp/tugnNbHqOTSLMiKBxJ2sVu+X1IxX/IMUxKt86VxIhEDifvQ7kI54IZWqUitKiBoS9zb0O6zX+h0YwdC5sFpJo7QYITGOOEKKe7NWlQkT7kG5KHECiTE9Pv7xGaGwxE0pF6VPIDEyhSVura5D29oXdQLUkEhCLCxl2ftaELiKuAZu+OsicQKJhBRFeHl5Hy1ZSMslhr+ilxcX05sMU3YkJy6kNePPeWzRPV0CiUTpd+AMvsdMmkAiXUO/w3odikK/A5+x0HUG1JBIXZ6/74BnIS0fstB1NoyQmIZhIa0JGQZxd8TYC8MMCCSmRGGJQdw/XhTNiUBiYmJhKV69pihc03OJDJTnSiAxSXE7srjssaosNFkKnQvzpqmBCYvFgzha0u8wbzoXlsAIicmL/Q51HeraNM4M9X2o6xCUixZAIDETVaXfYYaUixZFIDEfFtLOSdOEpgmPjyboFuTkQGrbtu/7EEKWZcWBUnJd123blmVZfv1IaJom3vjkUXAVcSFt07y3Bed5KEvfryejrt/LgWVpjm5xTg6kruuqqgoh1HV9KFqqqqqqagihqPR9lTsaLq3UdaGuQ9+Hvne9pUS1bRg+LeLWuizTyYGUff2qmZ3ynTPLss1mk+f5oRHSdrsdbud5nhulcyV5/kdJqW3Deh2y7D2cjNVH1HWhad4vMlKW77tDMV1d13VfW1232+15I5DPAqnv+3bnWmlFUZwUQruKoog5tDdsGjw8PBhCcWtF8UcIDXN6WRbKUqHiHvo+NM37jFyeh6oylTofVxlIfBZIWZZ9GxL91+um9TsXUBuqShceDdzNMH0Xu4rjn3P8lOS6YntClr1nvzPMISdP2eV5Hkc5u2H49PQUQnh9fY0/xsaHtm2HCbr1eh3n63rXASUxu23isUMv/nJ3OMWpYllomB3VnsAxfnh7exv7GEIIoWkaU3YkJe5SE8zpHW23LKTtfsnO/jy3Dgk+tjun1zTvS241ke8ZykJ9/77XrZPD2QQSfEeW/dGLPDSRh7DoPr294aOyEFchkOAE3zaRh8XM6cUwjjRqcwsCCc40dD3Eaas4coq1k9lMW+21IAohbkogwaXinF60O6c33cL+sH9P/KfNJl9JnECCa9qb04ubVYcpFJzaNrTtH0erLMT9CSS4ld2VTHX9vi4nwVJTbE8oCpftYGQCCe4h5QHHROcVmR+XMAcgCQIJgCQIJACSIJAASIJAAiAJAgmAJAgkAJIgkABIgkACIAkCCYAkCCQAkiCQAEiCQAIgCQIJgCQIJACSIJAASIJAAiAJAgmAJAgkAJIgkABIgkACIAkCCYAkCCQAkiCQAEiCQAIgCQIJgCQIJACSIJAASIJAAiAJAgmAJAgkAJKQSiD9/vvvYx/CNHRdN/YhTIMTdTzn6khO1JHO/jxPJZB+++23sQ9hGrwljuREHc+5OpITdaSzP89/vPy127bt+z6EkGVZURQf3qeu63iHT+4DwJJdYYTUdV1ZlmVZfvL1IcuyqqrKsozRBQB7rhBIWZbt3fhWWZYhhL7vm6a5/BUBmJ/Tpuz6vm/bdvixKIpPQujbx67X69VqdegOf/3rX+ONL1++/PTTTycd2KII9SM5Ucdzro7kRB3y+++/X94KcFogZVkWxzq7hlm43em4oaoUf2zbtq7rl5eXQ8/8yy+/nHQkAMzMD29vbxc+xYdNDX/5y19CCK+vr/EO6/X68fEx3v5kkATAYl0hkADgcqmsQwJg4QQSAEkQSAAk4Qo7NVzomI0elik2mPZ9X1XVh3fo+76u6zzPw9eVXsv03RMVxeaa2GizWHVdx97XT/5ghrekP6rP/6g2m02e513XxY4torqu27aNuyWc/OC3sb2+vu7d4O3t7ddff/3111/f3t7++c9/xhvf+vvf//7h7aWJ//bhjB2yWq2WfJbedk7R5+fBO9G773LnnZPxp+yO2ehhgbqui0OfoigO7ck0fAHZbDYGl/G76qH/ul6vPx8/LcHwR/WJ+N22aZol7/J1zLuvKIrNZtM0jXffFY0fSFzo+fn5pC0zFih+tvrUOEbf9/ETduGZdIy4UYC33hWNH0gfbvTA8H2/bdvdb7V739een5+rqvJRG775+r93ouInbNu2S76CwKFB5O4vi6KI32+qqtrdJ2xRjnn3tW0bx9xLPlFXN34g5XneNE3TNN+dTFiU+JZommb3c7Zt2z//+c/DX//z83OWZXHL2oVvsRXDZvdP6Oeff16v1/F2/CbrDyzP87Zt9/5U1uv1zz//PPwYrxETZ6IW29RwzLuvLMu6rpumqet6sSfqQ/FvrG3bM3LaTg0AJGH8ERIABIEEQCIEEgBJEEgAJEEgAZAEgQRAEv4/OL6FwllWID4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=zeros(N+2,1);\n",
    "i=1;\n",
    "while (i <= N+2)\n",
    "  x(i) = (i-1)/(N+1);\n",
    "  i=i+1;\n",
    "endwhile\n",
    "\n",
    "plot(x,y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On pourra tester les expériences suivantes :\n",
    "\n",
    "- faire varier $N$ dans $\\{1,5,10,20\\}$, \n",
    "\n",
    "- faire varier $\\rho$ dans $\\{10^{-1}, 5 \\times 10^{-2}, 10^{-2}, 5 \\times 10^{-3}, 10^{-3} \\}$\n",
    "\n",
    "- faire varier $\\varepsilon$ dans $\\{ 10^{-2}, 10^{-4}, 10^{-6}\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour pouvoir prendre des valeurs plus grandes de $N$, il faut faire appel à la recherche linéaire de Wolfe afin de\n",
    "calculer un meilleur pas que le pas constant à chaque iteration de la boucle principale.\n",
    "\n",
    "Rappelons que cette méthode fait appel à deux fonctions : une fonction principale Wolfe_search qui cherche un\n",
    "intervalle $[t_g, t_d]$ dans lequel il existe au moins un pas qui satisfait les conditions de Wolfe, et une fonction\n",
    "Zoom_search qui cherche un pas $\\rho$ dans dans cet intervalle.\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rappelons que dans la méthode de Wolfe, il s'agit de trouver un pas $\\rho$ tel que\n",
    "\n",
    "$$\n",
    "q(\\rho) \\leq q(0) + c_2 \\,q'(0)\\,\\rho \\quad \\quad \\text{et} |q'(\\rho)| \\leq - c_1 \\,q'(0) \n",
    "$$\n",
    "\n",
    "où $0 < c_1 < c_2 < 1$ sont des paramètres fixés, et $q$ est la fonction définie sur $\\mathbb{R}_+$ par\n",
    "\n",
    "$$\n",
    "q(t) \\, = \\, J( y - t\\, \\nabla J (y)).\n",
    "$$\n",
    "\n",
    "On rappelle que dans ce cas $q'(t) = - \\nabla J( y - t\\, \\nabla J (y)) \\cdot \\nabla J (y)$, et en particulier\n",
    "$q'(0) = - \\lVert \\nabla J (y) \\rVert^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans les fonctions suivantes, on doit passer les paramètres suivants :\n",
    "- le pas $\\rho$ utilisé à l'itération précédente;\n",
    "- le vecteur $y$ de l'itération précédente;\n",
    "- $g = \\nabla J (y)$ ;\n",
    "- $qp_0 = q'(0) = - \\lVert \\nabla J (y) \\rVert^2$ ;\n",
    "- les paramètres $c_1$ et $c_2$, ainsi que le nombre maximal d'itérations dans la recherche de Wolfe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "function t=Wolfe_search(rho,y,g,qp_0,c1,c2,nmax)\n",
    "  q_0=funcJ(y);\n",
    "  c1qp=c1*qp_0;\n",
    "  c2qp=c2*qp_0;\n",
    "  test=0;\n",
    "  tg=0.;\n",
    "  td=2*rho;\n",
    "  t=rho;\n",
    "  while (test==0)\n",
    "    qt=funcJ(y-t*g);\n",
    "    if (qt>(q_0+c2qp*t))\n",
    "      test =1;\n",
    "      t = Zoom_search(tg,t,y,g,q_0,c1qp,c2qp,nmax);\n",
    "    else\n",
    "      qp=-gradJ(y-t*g)'*g;\n",
    "      if (abs(qp) <= -c1qp)\n",
    "        test=2;\n",
    "      elseif (qp >=0)\n",
    "        test=3;\n",
    "        t = Zoom_search(tg,t,y,g,q_0,c1qp,c2qp,nmax);\n",
    "      else\n",
    "        tg=t;\n",
    "        td=10*td;\n",
    "        t=0.5*(tg+td);\n",
    "      endif\n",
    "    endif\n",
    "  endwhile\n",
    "endfunction\n",
    "\n",
    "function t=Zoom_search(tg,td,y,g,q_0,c1qp,c2qp,nmax)\n",
    "  n=0;\n",
    "  testZ=0;\n",
    "  while and(testZ==0,n<nmax)\n",
    "    t=0.5*(tg+td);\n",
    "    qt=funcJ(y-t*g);\n",
    "    if (qt>(q_0+c2qp*t))\n",
    "      td=t;\n",
    "    else    \n",
    "      qp=-gradJ(y-t*g)'*g;\n",
    "      if (abs(qp) <= -c1qp)\n",
    "        testZ=1;\n",
    "      elseif (qp >=0)\n",
    "        td=t;\n",
    "      else\n",
    "        tg=t; \n",
    "      endif\n",
    "    endif\n",
    "    n=n+1;\n",
    "  endwhile\n",
    "endfunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A vous d'intégrer cette méthode dans une méthode de plus grande pente !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "4.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}