diff --git a/2.1 Interpolation de donnee.ipynb b/2.1 Interpolation de donnee.ipynb
new file mode 100644
index 0000000..d7e1923
--- /dev/null
+++ b/2.1 Interpolation de donnee.ipynb
@@ -0,0 +1,363 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 2 Interpolation et approximation de données\n",
+ "\n",
+ "## 2.1 Position du problème"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Interpolation de données\n",
+ "\n",
+ "Soit $n \\geq 0$ un nombre entier. Etant donnés $(n+1)$ noeuds \n",
+ "distincts $t_0$, $t_1$,$\\dots$ $t_n$ et $(n+1)$ valeurs $y_0$,\n",
+ "$y_1$,$\\dots$ $y_n$, on cherche un polynôme $p$\n",
+ " de degré $n$, tel que\n",
+ "\n",
+ "$$p(t_j)=y_j \\qquad \\text{ pour } \\, 0\\leq j\\leq n.$$\n",
+ "\n",
+ "**Exemple** On cherche le polynôme $\\Pi_n$ de degré $n=4$ tel que $\\Pi_n(t_j) = y_j, j =1,...,5$ avec \n",
+ "les données suivantes : \n",
+ "\n",
+ "| $t_k$ | $y_k$ |\n",
+ "| --- | --- |\n",
+ "| 1 | 3 |\n",
+ "| 1.5 | 4 |\n",
+ "| 2 | 2 |\n",
+ "| 2.5 | 5 |\n",
+ "| 3 | 1 |\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2.1.1 Interpolation\n",
+ "\n",
+ "\n",
+ "## 2.2.1 Unicité\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# importing libraries used in this book\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "\n",
+ "# Some data given: t=1, 1.5, 2, 2.5, 3 and y = 3,4,2,5,1 \n",
+ "t = np.linspace(1, 3, 5) # equivalent to np.array([ 1, 1.5, 2, 2.5, 3 ])\n",
+ "y = np.array([3, 4, 2, 5, 1])\n",
+ "\n",
+ "# Plot the points using matplotlib\n",
+ "plt.plot(t, y, 'ro')\n",
+ "\n",
+ "plt.xlabel('t'); plt.ylabel('y'); #plt.title('data')\n",
+ "plt.legend(['data'])\n",
+ "plt.show()\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Si ce polynôme existe, on le note $p=\\Pi_n$. On appelle $\\Pi_n$ le\n",
+ "polynôme d'interpolation des valeurs $y_j$ aux noeuds $t_j,\\, j=0,\\dots,n$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plot the interpolating function\n",
+ "\n",
+ "# Defining the polynomial function \n",
+ "def p(t):\n",
+ " # coefficients of the interpolating polynomial\n",
+ " a = np.array([-140.0, 343.0, -872./3., 104.0, -40./3.])\n",
+ " \n",
+ " # value of the polynomial in all the points t\n",
+ " return a[0] + a[1]*t + a[2]*(t**2) + a[3]*(t**3) + a[4]*(t**4)\n",
+ "\n",
+ "\n",
+ "# points used to plot the graph \n",
+ "z = np.linspace(1, 3, 100)\n",
+ "\n",
+ "plt.plot(t, y, 'ro', z, p(z))\n",
+ "plt.xlabel('t'); plt.ylabel('y'); #plt.title('data')\n",
+ "plt.legend(['data','$\\Pi_2(t)$'])\n",
+ "plt.show()\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Interpolation de fonctions\n",
+ "\n",
+ "\n",
+ "Soit $f\\in C^0(I)$ et $t_0,\\ldots,t_n\\in I$. \n",
+ "Si on prend $$y_j=f(t_j),\\quad 0\\leq j\\leq n,$$ \n",
+ "alors le polynôme d'interpolation\n",
+ "$\\Pi_n(t)$ est noté $\\Pi_n f(t)$ et est appelé l'interpolant de $f$ aux\n",
+ "noeuds $t_0$,$\\dots$ $t_n$.\n",
+ "\n",
+ "**Exemple** Soient $$t_1=1, t_2=1.75, t_3=2.5, t_4=3.25, t_5=4$$ les points d'interpolation et $$f(t) = t \\sin(2\\pi t).$$ On cherche l'interpolant $\\Pi_n f$ de degré $n=4$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# defining the fonction that we want to interpolate\n",
+ "def f(t):\n",
+ " return t*np.sin(t*2.*np.pi)\n",
+ "\n",
+ "# The interpolation must occour at points t=1, 1.75, 2.5, 3.25, 4\n",
+ "t = np.linspace(1, 4, 5)\n",
+ "\n",
+ "# points used to plot the graph \n",
+ "z = np.linspace(1, 4, 100)\n",
+ "\n",
+ "# (Complete the code below)\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plot the interpolating function\n",
+ "\n",
+ "# Defining the polynomial function \n",
+ "def p(t):\n",
+ " # coefficients of the interpolating polynomial\n",
+ " a = np.array([0, 7.9012, -13.037, 5.9259, -0.79012])\n",
+ " \n",
+ " # value of the polynomial in all the points t\n",
+ " return a[0] + a[1]*t + a[2]*(t**2) + a[3]*(t**3) + a[4]*(t**4)\n",
+ "\n",
+ "\n",
+ "# points where to evaluate the polynomial\n",
+ "z = np.linspace(1, 4, 100)\n",
+ "\n",
+ "# (Complete the code below)\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Une méthode naïve\n",
+ "\n",
+ "Il est possible d'écrire un système d'équations et de trouver les coefficients de manière directe.\n",
+ "Comme expliqué dans le MOOC, ce n'est pas toujours la meilleure solution.\n",
+ "\n",
+ "Nous cherchons les coefficients du polynôme $p(t) = a_0 + a_1 t + ... + a_n t^n$ qui satisfait les $(n+1)$ équations $p(t_k) = y_k, k=0,...,n$, c'est-à-dire\n",
+ "\n",
+ "$$a_0 + a_1 t_k + ... + a_n t_k^n = y_k, k=0,...,n$$\n",
+ "\n",
+ "Ce système s'écrit sous forme matricielle\n",
+ "\n",
+ "$$\\begin{pmatrix}\n",
+ "1 & t_0 & t_0^2 & \\cdots & t_0^n \\\\\n",
+ "\\vdots & & & & \\vdots \\\\\n",
+ "1 & t_n & t_n^2 & \\cdots & t_n^n\n",
+ "\\end{pmatrix}\n",
+ "\\begin{pmatrix} a_0 \\\\ \\vdots \\\\ a_n \\end{pmatrix}\n",
+ "=\\begin{pmatrix} y_0 \\\\ \\vdots \\\\ y_n \\end{pmatrix}$$\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2.1.3 Matrice de Vandermonde\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Exemple** On cherche les coefficients du polynôme d'interpolation de degré $n=4$ \n",
+ "des valeurs suivantes \n",
+ "\n",
+ "| $t_k$ | $y_k$ |\n",
+ "| --- | --- |\n",
+ "| 1 | 3 |\n",
+ "| 1.5 | 4 |\n",
+ "| 2 | 2 |\n",
+ "| 2.5 | 5 |\n",
+ "| 3 | 1 |\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Some data given: t=1, 1.5, 2, 2.5, 3 and y = 3,4,2,5,1 \n",
+ "t = np.linspace(1, 3, 5)\n",
+ "y = np.array([3, 4, 2, 5, 1])\n",
+ "n = t.size - 1\n",
+ "\n",
+ "# The following is a trick to put toghether the matrix A. Don't need to learn by heart ...\n",
+ "# np.tile(t, (n+1, 1)) : repete n+1 times of the array t \n",
+ "# .T : transpose\n",
+ "# np.linspace(0,n,n+1) : [0,...,n+1]\n",
+ "# np.tile : again: repete n+1 times the array [0,...,n+1]\n",
+ "# np.power : element by element power funcion\n",
+ "A = np.power( np.tile(t, (n+1, 1)).T , np.tile(np.linspace(0,n,n+1), (n+1, 1)))\n",
+ "# print(A)\n",
+ "\n",
+ "# (Complete the code below)\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Now we can define the polynomial\n",
+ "p = lambda t : a[0] + a[1]*t + a[2]*(t**2) + a[3]*(t**3) + a[4]*(t**4)\n",
+ "\n",
+ "# points used to plot the graph \n",
+ "z = np.linspace(1, 3, 100)\n",
+ "\n",
+ "# (Complete the code below)\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Alternatives : polyfit et polyval\n",
+ "\n",
+ "Les fonctions `polyfit` et `polyval` de `numpy` font essentiellement la même chose que les paragraphes ci-dessous. Plus tard nous verrons des méthodes plus performantes.\n",
+ "\n",
+ "`a = numpy.polyfit(t, y, n, ... )` :\n",
+ "\n",
+ " * input : $t,y$ les données à interpoler, $n$ le degré du polynôme recherché\n",
+ " * outut : les coefficients du polynôme _dans l'ordre inverse_ de ce que nous avons vu !\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Some data given: t=1, 1.5, 2, 2.5, 3 and y = 3,4,2,5,1 \n",
+ "t = np.linspace(1, 3, 5)\n",
+ "y = np.array([3, 4, 2, 5, 1])\n",
+ "n = t.size - 1\n",
+ "\n",
+ "a = np.polyfit(t,y,n)\n",
+ "\n",
+ "# Now we can define the polynomial, with coeffs in the reverse order !\n",
+ "p = lambda t : a[4] + a[3]*t + a[2]*(t**2) + a[1]*(t**3) + a[0]*(t**4)\n",
+ "\n",
+ "# We can also use polyval instead !\n",
+ "# np.polyval(a,t)\n",
+ "\n",
+ "# points used to plot the graph \n",
+ "z = np.linspace(1, 3, 100)\n",
+ "\n",
+ "plt.plot(t, y, 'ro', z, p(z), '.', z, np.polyval(a,z))\n",
+ "plt.xlabel('t'); plt.ylabel('y'); #plt.title('data')\n",
+ "plt.legend(['data','p(t)','polyval'])\n",
+ "plt.show()\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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.7.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/2.2 Interpolation de Lagrange.ipynb b/2.2 Interpolation de Lagrange.ipynb
index c42f6af..4faa0ff 100644
--- a/2.2 Interpolation de Lagrange.ipynb
+++ b/2.2 Interpolation de Lagrange.ipynb
@@ -1,190 +1,200 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## 1.2 Interpolation de Lagrange\n",
+ "## 2.2 Interpolation de Lagrange\n",
"\n",
"### Base de Lagrange\n",
"\n",
"On considère les polynômes $\\varphi_k$, $k=0,\\dots, n$ de\n",
"degré $n$ tels que\n",
"\n",
"$$\\varphi_k(t_j)=\\delta_{jk}, \\qquad k,j=0,\\dots, n,$$\n",
"\n",
"où $\\delta_{jk}=1$ si $j=k$ et $\\delta_{jk}=0$ si $j\\neq k$.\n",
"Explicitement, on a\n",
"\n",
"$$\\varphi_k(t)=\\prod_{j=0,j\\ne k}^n\\dfrac{(t-t_j)}{(t_k-t_j)}.$$"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2.2.1 Base de Lagrange\n",
+ ""
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# importing libraries used in this book\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercice (théorique)** Vérifiez que\n",
"\n",
"1. $B = \\{\\varphi_k, k=0,\\dots, n\\}$ est une base de $P_n(\\mathbb R)$\n",
"2. Chaque polynôme $\\varphi_k$ est de degré $n$\n",
"3. $\\varphi_k(t_j)=\\delta_{jk}, \\qquad k,j=0,\\dots, n$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Defining the Lagrange basis functions\n",
"def phi(t,k,z):\n",
" # the input variables are:\n",
" # t : the interpolatory points\n",
" # k : which basis function\n",
" # z : where to evaluate the function\n",
" \n",
"\n",
" # careful, there are n+1 interpolation points!\n",
" n = t.size - 1 \n",
"\n",
" # init result to one, of same type and size as z\n",
" result = np.zeros_like(z) + 1\n",
"\n",
" # first few checks on k:\n",
" if (type(k) != int) or (t.size < 1) or (k > n) or (k < 0):\n",
" raise ValueError('Lagrange basis needs a positive integer k, smaller than the size of t')\n",
" \n",
" # loop on n to compute the product\n",
" for j in range(0,n+1) :\n",
" if (j == k) :\n",
" continue\n",
" if (t[k] == t[j]) :\n",
" raise ValueError('Lagrange basis: all the interpolation points need to be distinct')\n",
"\n",
"# (Complete the code below)\n",
" result = \n",
"\n",
" return result\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exemple\n",
"\n",
"Pour $n=2$, $t_0=-1$, $t_1=0$, $t_2=1$, les polynômes de la base de Lagrange\n",
"sont\n",
"\n",
"\\begin{equation*}\n",
"\\begin{array}{lcl}\n",
"\\varphi_0(t)&=&\\displaystyle\\dfrac{(t-t_1)(t-t_2)}{(t_0-t_1)(t_0-t_2)}=\n",
"\\dfrac{1}{2}t(t-1),\\\\[2mm]\n",
"\\varphi_1(t)&=&\\displaystyle\\dfrac{(t-t_0)(t-t_2)}{(t_1-t_0)(t_1-t_2)}=\n",
"-(t+1)(t-1),\\\\[2mm]\n",
"\\varphi_2(t)&=&\\displaystyle\\dfrac{(t-t_0)(t-t_1)}{(t_2-t_0)(t_2-t_1)}=\n",
"\\dfrac{1}{2}t(t+1).\n",
"\\end{array}\n",
"\\end{equation*}\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the Lagrange Basis functions \n",
"t = np.linspace(-1., 1, 3)\n",
"z = np.linspace(-1.1, 1.1, 100)\n",
"\n",
"plt.plot(z, phi(t,0,z), 'g', z, phi(t,1,z), 'r', z, phi(t,2,z),':')\n",
"\n",
"plt.xlabel('t'); plt.ylabel('$\\\\varphi_{k}(t)$'); plt.title('Lagrange basis functions')\n",
"plt.legend(['$\\\\varphi_{0}$','$\\\\varphi_{1}$','$\\\\varphi_{2}$'])\n",
"plt.grid(True)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exercice\n",
"\n",
"Visualisez la base de Lagrange associée aux points $t_k = 1, 1.5, 2, 2.5, 3$.\n",
"Evaluez sur le graphique les valeurs de $\\varphi_k(t_j)$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the Lagrange Basis functions \n",
"t = np.linspace(1., 3, 5)\n",
"z = np.linspace(0.9, 3.1, 100)\n",
"\n",
"# (Complete the code below)\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Polynôme d'interpolation\n",
"\n",
"Le polynôme d'interpolation $\\Pi_n$ des\n",
"valeurs $y_j$ aux noeuds $t_j$, $j=0,\\dots,n$, s'écrit\n",
"\\begin{equation}\\label{e:int_lagr}\n",
" \\Pi_n(t)=\\sum_{k=0}^n y_k \\varphi_k(t),\n",
"\\end{equation}\n",
"car il vérifie $\\Pi_n(t_j)=\\sum_{k=0}^n y_k \\varphi_k(t_j)=y_j$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
}
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "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.7.7"
+ "version": "3.7.16"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
diff --git a/2.3 Interpolation de fonction.ipynb b/2.3 Interpolation de fonction.ipynb
index e0a6467..de4de3c 100644
--- a/2.3 Interpolation de fonction.ipynb
+++ b/2.3 Interpolation de fonction.ipynb
@@ -1,401 +1,417 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## 1.3 Interpolation d'une fonction continue\n",
+ "## 2.3 Interpolation d'une fonction continue\n",
"\n",
"Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $t_0,\\ldots,t_n\\in [a,b]$ des noeuds distincts. Le polynôme d'interpolation\n",
"$\\Pi_n(t)$ est noté $\\Pi_n f(t)$ et est appelé l'interpolant de $f$ aux\n",
"noeuds $t_0,\\dots,t_n$.\n",
"\n",
"\n",
"Si on prend $$y_k=f(t_k), k= 0,..., n,$$ \n",
"alors on aura\n",
"\\begin{equation*}\n",
"\\Pi_nf(t)=\\sum_{k=0}^n\n",
"f(t_k)\\varphi_k(t).\n",
"\\end{equation*}"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 2.3.1 Interpolation de fonction\n",
+ "\n",
+ "\n",
+ "## 2.3.2 Interpolation de Lagrange\n",
+ "\n",
+ "\n"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exercice\n",
"\n",
"Ecrivez une fonction Python qui a la définition suivante, en utilisant la fonction `phi` définie plus haut.\n",
"Ecrivez aussi un petit test sur la base de l'exercice précédent.\n",
"\n",
"```python\n",
"# Lagrange Interpolation of data (t,y), evaluated at ordinate(s) z\n",
"def LagrangePi(t,y,z):\n",
" # the input variables are:\n",
" # t : the interpolatory points\n",
" # y : the corresponding data at the points t\n",
" # z : where to evaluate the function\n",
" \n",
"```\n",
"\n",
"Utilisez le fait que $\\{\\varphi_k, k = 0,...,n\\}$ est une base des polynômes de degré $\\leq n$\n",
"et que le vecteur $y$ représente les coordonnéee du polynôme d'interpolation recherché par rapport à cette base,\n",
"c'est-à-dire\n",
"$$\\Pi_n (z) = y_0 \\varphi_0 + ... + y_n \\varphi_n$$\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# importing libraries used in this book\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from InterpolationLib import LagrangeBasis as phi"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Lagrange Interpolation of data (t,y), evaluated at ordinate(s) z\n",
"def LagrangePi(t,y,z):\n",
" # the input variables are:\n",
" # t : the interpolatory points\n",
" # y : the corresponding data at the points t\n",
" # z : where to evaluate the function\n",
" \n",
" # (Complete the code below)\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
" return result "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# vecteur des points d'interpolation\n",
"t = np.linspace(0, 1, 5)\n",
"\n",
"# vecteur des valeurs\n",
"y = np.array([3.38, 3.86, 3.85, 3.59, 3.49])\n",
"\n",
"z = np.linspace(-0.1, 1.1, 100)\n",
"plt.plot(t, y, 'ro', z, LagrangePi(t,y,z))\n",
"\n",
"plt.xlabel('t'); plt.ylabel('y');\n",
"plt.legend(['data','p(t)'])\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Erreur d'interpolation\n",
"\n",
"Soient $t_0$, $t_1$, $\\ldots$, $t_n$,\n",
"$(n+1)$ nœuds *équirépartis* dans $I=[a,b]$ et soit $f\\in C^{n+1}(I)$.\n",
"Alors\n",
"\\begin{equation}\n",
" \\max_{x\\in I} |f(t) - \\Pi_n f(t)| \\leq \\dfrac{1}{2(n+1)}\n",
"\\left(\\frac{b-a}{n}\\right)^{n+1}\\max_{x\\in I}|f^{(n+1)}(t)|.\n",
"\\end{equation}\n",
"\n",
"On remarque que l'erreur d'interpolation dépend de la dérivée\n",
"$(n+1)$-ième de $f$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exercice\n",
"\n",
"On considère les points d'interpolation $$t_0=1, t_1=1.75, t_2=2.5, t_3=3.25, t_4=4$$ \n",
"et la fonction $$f(t) = t \\sin(2\\pi t)$$.\n",
"\n",
"1. Calculez la base de Lagrange associée à ces points. D'abord sur papier, ensuite utilisez Python pour en dessiner le graphique.\n",
"\n",
"2. Calculez le polynôme d'interpolation $\\Pi_n$ à l'aide de la base de Lagrange. D'abord sur papier, ensuite avec Python.\n",
"\n",
"4. Quelle est l'erreur théorique d'interpolation ? \n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Comportement pour $n$ grand: eg, la fonction de Runge.\n",
"\n",
"Le fait que\n",
"$$\\lim_{n\\rightarrow\\infty} \\dfrac{1}{2(n+1)}\\left(\\dfrac{b-a}{n}\\right)^{n+1}=0$$\n",
"n'implique pas forcément que $\\max_{t \\in I} |E_nf(t)|$ tende vers zéro quand $n\\rightarrow\\infty$.\n",
"\n",
"Soit $f(t)=\\dfrac{1}{1+t^2}$, $t\\in [-5,5]$. Si on l'interpole dans des\n",
"noeuds équirépartis, l'interpolant présente des oscillations au voisinage des extrémités de l'intervalle.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Runge fonction\n",
"f = lambda t : 1./(1+t**2)\n",
"\n",
"# Values of N to use\n",
"Nrange = [3,5,10]\n",
"# plotting points\n",
"z = np.linspace(-5, 5, 100)\n",
"\n",
"# (Complete the code below)\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`sympy` est une librairie pour le calcul symbolique en Python. Nous allons l'utiliser pour étudier le comportement de l'erreur d'interpolation de la fonction de Runge."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Using symbolic python to compute derivatives\n",
"import sympy as sp\n",
"\n",
"# define x as symbol\n",
"x = sp.symbols('x')\n",
"\n",
"# define Runge function\n",
"f = 1/(1+x**2)\n",
"\n",
"# pretty print the 5th derivative of f\n",
"f5 = sp.diff(f, x,5)\n",
"# display f5\n",
"sp.init_printing(use_unicode=True)\n",
"display(f5)\n",
"\n",
"\n",
"# evalf can be used to compute the value of a function at a given point\n",
"print('5th derivative evaluated at 3 :')\n",
"print( f5.evalf(subs={x: 3}) )\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Pour définir une fonction qui accepte un array de valeur, il faut utiliser les lignes suivantes.\n",
"\n",
"Ensuite on peut aussi dessiner le graphe ..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# to evaluate a function at many given points, we need the following trick\n",
"diff_f_func = lambda t: float(sp.diff(f,x,k).evalf(subs={x: t}))\n",
"diff_f = np.vectorize(diff_f_func)\n",
"\n",
"# the derivative can be set with k (not very elegant...)\n",
"k = 4\n",
"print(diff_f(4.5))\n",
"\n",
"# plotting points\n",
"z = np.linspace(-5, 5, 100)\n",
"# plot the derivative between -5 and 5\n",
"\n",
"plt.plot(z,diff_f(z), 'b')\n",
"plt.xlabel('t'); plt.ylabel('y'); plt.title('Derivatives of Runge function')\n",
"plt.legend(Nrange)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"... ou évaluer le maximum pour plusieurs $n$ et voir le comportement de $\\max |f^{(n)}|$ en fonction de $n$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot max(abs(fn)) in the range -5,5\n",
"z = np.linspace(-5, 5, 100)\n",
"\n",
"Nmax = 10\n",
"maxValFn = np.zeros(Nmax)\n",
"for k in range(Nmax):\n",
" maxValFn[k] = np.max(np.abs(diff_f(z)))\n",
" \n",
"plt.plot(range(10), maxValFn)\n",
"\n",
"plt.yscale('log')\n",
"plt.xlabel('n'); plt.ylabel('$\\max|\\partial f|$');\n",
"plt.title('Max of $|\\partial^n f|$');\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Interpolation de Chebyshev\n",
"\n",
"\n",
"Pour chaque entier positif $n\\geq 1$, pour $i=0,\\dots n$, on note\n",
"$$\\hat{t}_i=-\\cos(\\pi i/n)\\in[-1,1]$$\n",
"**les points de Chebyshev** et on définit\n",
"\n",
"$$t_i=\\dfrac{a+b}{2}+\\dfrac{b-a}{2}\\hat{t}_i\\in[a,b],$$\n",
"\n",
"pour un intervalle arbitraire $[a,b]$. Pour une fonction continue $f\\in\n",
"C^1([a,b])$, le polynôme d'interpolation $\\Pi_n f$ de degré $n$ aux noeuds\n",
"$\\{t_i,i=0,\\ldots,n\\}$ converge uniformément vers $f$ quand\n",
"$n\\rightarrow \\infty$.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Chebichev points on the interval [-1,1]\n",
"\n",
"z = np.linspace(0,1, 100)\n",
"plt.plot(np.cos(np.pi*z), np.sin(np.pi*z))\n",
"\n",
"n =5\n",
"z = np.linspace(0,1, n+1)\n",
"plt.plot(np.cos(np.pi*z), np.sin(np.pi*z), 'o')\n",
"plt.plot(np.cos(np.pi*z), 0*z, 'x')\n",
"\n",
"for k in range(0,n+1) :\n",
" plt.plot([np.cos(np.pi*z[k]),np.cos(np.pi*z[k])],[0,np.sin(np.pi*z[k])],':')\n",
"\n",
"plt.axis('equal')\n",
"\n",
"plt.xlabel('t'); \n",
"plt.title('Chebyshev points')\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exemple\n",
"\n",
"On reprend le\n",
" même exemple mais on interpole la fonction de Runge dans les points de\n",
"Chebyshev. La figure montre les polynômes de\n",
"Chebyshev de degré $n=5$ et $n=10$. On remarque que les oscillations\n",
"diminuent lorsqu'on augmente le degré du polynôme.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Runge fonction\n",
"f = lambda t : 1./(1+t**2)\n",
"\n",
"# Values of N to use\n",
"Nrange = [3,5,10]\n",
"# plotting points\n",
"[a,b] = [-5,5]\n",
"z = np.linspace(a,b, 100)\n",
"\n",
"for n in Nrange :\n",
"\n",
" # Chebyshev points on [-1,1]\n",
" ht = -np.cos(np.pi*np.linspace(0,n,n+1)/n)\n",
" # mapped to [a,b]\n",
" t =(a+b)/2 + (b-a)/2*ht\n",
"\n",
" y = f(t);\n",
"\n",
" plt.plot(z, LagrangePi(t,y,z), ':')\n",
"\n",
"plt.plot(z,f(z), 'b')\n",
"plt.xlabel('t'); plt.ylabel('y'); plt.title('Chebyshev interpolation')\n",
"plt.legend(Nrange)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "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.7.7"
+ "version": "3.7.16"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
diff --git a/2.4 Interpolation par intervalles.ipynb b/2.4 Interpolation par intervalles.ipynb
index 073791f..1449885 100644
--- a/2.4 Interpolation par intervalles.ipynb
+++ b/2.4 Interpolation par intervalles.ipynb
@@ -1,292 +1,308 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## 1.4 Interpolation par intervalles\n",
+ "## 2.4 Interpolation par intervalles\n",
"\n",
"### Interpolation linéaire par morceaux\n",
"\n",
"Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $a=x_0<\\ldots\n",
+ "\n",
+ "## 2.4.2 Convergence\n",
+ "\n"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# importing libraries used in this book\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from InterpolationLib import PiecewiseLinearInterpolation as PiH1\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"La fonction `PiH1` implémente l'interpolation linéaire par morceaux pour des points équidistribués\n",
"\n",
"```python\n",
"def PiecewiseLinearInterpolation(a,b,N,f,z):\n",
" # the input variables are:\n",
" # a,b : x[0] = a, x[n] = b\n",
" # f : the corresponding data at the points x\n",
" # z : where to evaluate the function\n",
"```\n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercice\n",
"\n",
"Soient $x_k, k=0,...,4$ des points équidistribués sur l'intervalle $[1,5]$ et \n",
"$y_k = (3.38, 3.86, 3.85, 3.59, 3.49)$\n",
"les valeurs d'une fonction en ces points. \n",
"\n",
"* Dessinez le graphe de l'interpolateur par morceaux de cette fonction\n",
"* Calculez numériquement (sur papier) la valeur de $\\Pi^H_1(4.5)$ et vérifiez le résultat sur le graphique\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# intervalle d'interpolation\n",
"a = 1; b = 5\n",
"\n",
"# (Complete the code below)\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercice\n",
"\n",
"Dessinez le graphe de la fonction de Runge et l'interpolateur linéaire par morceaux \n",
"sur l'intervalle $[-5,5]$ pour $N=3,5,10$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# interval and function\n",
"a = -5; b = 5\n",
"f = lambda x : 1/(1+x**2)\n",
"\n",
"# (Complete the code below)\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Erreur d'interpolation linéaire par morceaux \n",
"\n",
"#### Théorème \n",
"\n",
"Soient $f\\in C^{2}([a,b])$, $H = \\frac{b-a}{N}$, $x_i = a + iH$ pour $i=0,...,N$.\n",
" \n",
"Soit $E^H_1f(t) = f(t) - \\Pi_1^{H} f(t)$, alors\n",
"$$\\max_{x\\in I}\\mid E^H_1f(t) \\mid\\leq \\dfrac{H^2}{4}\\max_{x \\in I} | f''(t)|.$$\n",
"\n",
"**Preuve**\n",
"D'après la formule, sur chaque intervalle $I_i=[x_i,x_{i+1}]$, on a\n",
" \n",
"$$\\max_{x\\in[x_i,x_{i+1}]}\\mid E^H_1f(t)\\mid \\leq\n",
" \\dfrac{H^2}{2(1+1)}\\max_{x\\in I_i}\\mid f''(t)\\mid.$$\n",
"\n",
"**Remarque**\n",
"On peut aussi montrer que, si l'on utilise un polynôme\n",
"de degré $n$ ($\\geq 1$) et si l'on dénote $E^H_nf(t) = f(t) - \\Pi_n^{H} f(t)$,\n",
" dans chaque sous-intervalle $I_i$, on\n",
"trouve\n",
"\n",
"$$\\max_{x\\in I}\\mid E^H_nf(t) \\mid \\leq \\dfrac{H^{n+1}}{2(n+1)} \\max_{x \\in I} |f^{(n+1)} (t)| \\, .$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Interpolation quadratique par morceaux\n",
"\n",
"Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $a=x_0<\\ldots$(m+1)$ points \n",
"~~distincts~~ $t_0$, $t_1$,$\\dots$ $t_m$ et $(m+1)$ valeurs $y_0$,\n",
"$y_1$,$\\dots$ $y_n$, on cherche un polynôme $p$\n",
" de degré $n , tel que\n",
"\n",
"$$\\sum_{j=0}^m |y_j - p(t_j)|^2 \\qquad \\text{ soit le plus petit possible}.$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On appelle **polynôme aux moindres carrés de degré $n$**\n",
"le polynôme $\\hat{p}_n(t)$ de degré $n$ tel que\n",
"\\begin{equation}\n",
" \\sum_{j=0}^m |y_j - \\hat{p}_n(t_j)|^2 \\leq \\sum_{j=0}^m |y_j - p_n(t_j)|^2\n",
" \\qquad \\forall p_n(t) \\in \\mathbb{P}_n\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nous cherchons les coefficients du polynôme $p(t) = a_0 + a_1 t + ... + a_n t^n$ qui satisfait *au mieux* les $(n+1)$ équations $p(t_k) = y_k, k=0,...,n$, c'est-à-dire\n",
"\n",
"$$a_0 + a_1 t_k + ... + a_n t_k^n = y_k, k=0,...,m$$\n",
"\n",
"Ce système s'écrit sous forme matricielle \n",
"\n",
"$$\\begin{pmatrix}\n",
"1 & t_0 & \\cdots & t_0^n \\\\\n",
"\\vdots & & & \\vdots \\\\\n",
"1 & t_m & \\cdots & t_m^n\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix} a_0 \\\\ \\vdots \\\\ a_n \\end{pmatrix}\n",
"=\\begin{pmatrix} y_0 \\\\ \\vdots \\\\ y_m \\end{pmatrix}$$\n",
"\n",
"Puisque $m\n",
+ "\n"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exemple\n",
"\n",
"On considère un test mécanique pour\n",
"établir le lien entre contraintes ($MPa=100 N/cm^2$)\n",
"et déformations relatives (cm/cm) d'un échantillon de tissu biologique (disque intervertébral, selon P. Komarek,\n",
"Ch. 2 de *Biomechanics of Clinical Aspects of Biomedicine*, 1993, J. Valenta ed., Elsevier).\n",
"\n",
"![title](Figures/disque.png)\n",
"\n",
"\n",
"On cherche à approximer au sens des moindres carrés avec un polynôme $\\hat p_n$ de degré $n=1,2,3$.\n",
"Les mesures effectuées sont les suivantes \n",
"```Python\n",
"sigma = np.array([0.00, 0.06, 0.14, 0.25, 0.31, 0.47, 0.50, 0.70]);\n",
"epsilon = np.array([0.00, 0.08, 0.14, 0.20, 0.22, 0.26, 0.27, 0.29]);\n",
"```\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# importing libraries used in this book\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# data given:\n",
"sigma = np.array([0.00, 0.06, 0.14, 0.25, 0.31, 0.47, 0.50, 0.70]);\n",
"epsilon = np.array([0.00, 0.08, 0.14, 0.20, 0.22, 0.26, 0.27, 0.29]);\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"m = sigma.size - 1\n",
"\n",
"# degree of the polynomial\n",
"n = 3;\n",
"\n",
"t = sigma\n",
"# The following is a trick to put toghether the matrix B.\n",
"# np.tile(t, (n+1, 1)) : repete n+1 times of the array t \n",
"# .T : transpose\n",
"# np.linspace(0,n,n+1) : [0,...,n+1]\n",
"# np.tile : again: repete m+1 times the array [0,...,n]\n",
"# np.power : element by element power funcion\n",
"B = np.power( np.tile(t, (n+1, 1)).T , np.tile(np.linspace(0,n,n+1), (m+1, 1)))\n",
"# print(A)\n",
"\n",
"# (Complete the code below)\n",
"\n",
"\n",
"\n",
"# print the coefficients on screen\n",
"print('The coefficients a_0, ..., a_n are')\n",
"# print(a)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def polynomial(a,t):\n",
" n = a.size-1\n",
" # \\hat p = a_0 + a_1 t + ... + a_n t^n\n",
" # is equal to the scalar product between the vectors a and (1, t, ..., t^n) :\n",
" return np.power( np.tile(t, (n+1, 1)).T , np.linspace(0,n,n+1)).dot(a)\n",
"\n",
"\n",
"# (Complete the code below)\n",
"\n",
"\n",
"plt.xlabel('$\\sigma$'); plt.ylabel('$\\epsilon$');\n",
"plt.legend(['data','$\\hat p_n$'])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Exercice\n",
"\n",
"Depuis https://hsso.ch/fr/2012/b/14 on a téléchargé les données de la population suisse dans\n",
"le fichier `Data/PopulationSuisse.csv`.\n",
"\n",
"Approximer l'évolution de la population avec un polynôme de degré $n=1,2,3,7$. \n",
"\n",
"Ensuite faites l'hypothèse de croissance exponentielle de la population, c'est-à-dire $p(x) = e^{a_1 x+a_0}$ où $a_0$ et $a_1$ sont des paramètres. Comment utiliser l'approximation polynômiale dans ce contexte ? "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Data of the population has been dowloaded from https://hsso.ch/fr/2012/b/14\n",
"# into the file Data/PopulationSuisse.csv\n",
"import pandas as pd \n",
"\n",
"# Read data from file 'PopulationSuisse.csv' \n",
"data = pd.read_csv(\"Data/PopulationSuisse.csv\") \n",
"# Preview the first 5 lines of the loaded data \n",
"data.head()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "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.7.7"
+ "version": "3.7.16"
}
},
"nbformat": 4,
"nbformat_minor": 4
}