arthur/ArtStudies
1
0
Fork 0
mirror of https://github.com/ArthurDanjou/ArtStudies.git synced 2026-01-14 13:54:06 +01:00
Code Issues Packages Projects Releases Wiki Activity

Refactor code for improved readability and consistency across notebooks

Browse Source 
- Standardized spacing around operators and function arguments in TP7_Kmeans.ipynb and neural_network.ipynb.
- Enhanced the formatting of model building and training code in neural_network.ipynb for better clarity.
- Updated the pyproject.toml to remove a specific TensorFlow version and added linting configuration for Ruff.
- Improved comments and organization in the code to facilitate easier understanding and maintenance.
This commit is contained in:
Arthur Danjou
2025-07-01 20:46:08 +02:00
parent e273cf90f7
commit f94ff07cab
34 changed files with 5713 additions and 5047 deletions
Download Patch File Download Diff File Expand all files Collapse all files

68
L3/Analyse Matricielle/TP1_Methode_de_Gauss.ipynb
View File

@@ -76,24 +76,40 @@
],
"source": [
"import numpy as np\n",
"\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"u = np.array([1,2,3,4,5])\n",
"v = np.array([[1,2,3,4,5]])\n",
"su=u.shape\n",
"sv=v.shape\n",
"u = np.array([1, 2, 3, 4, 5])\n",
"v = np.array([[1, 2, 3, 4, 5]])\n",
"su = u.shape\n",
"sv = v.shape\n",
"ut = np.transpose(u)\n",
"vt = np.transpose(v)\n",
"vt2 = np.array([[1],[2],[3],[4],[5]])\n",
"A = np.array([[1,2,0,0,0],[0,2,0,0,0],[0,0,3,0,0],[0,0,0,4,0],[0,0,0,0,5]])\n",
"B = np.array([[1,2,3,4,5],[2,3,4,5,6],[3,4,5,6,7],[4,5,6,7,8],[5,6,7,8,9]])\n",
"d=np.diag(A)\n",
"dd=np.array([np.diag(A)])\n",
"dt=np.transpose(d)\n",
"ddt=np.transpose(dd)\n",
"Ad=np.diag(np.diag(A))\n",
"vt2 = np.array([[1], [2], [3], [4], [5]])\n",
"A = np.array(\n",
" [\n",
" [1, 2, 0, 0, 0],\n",
" [0, 2, 0, 0, 0],\n",
" [0, 0, 3, 0, 0],\n",
" [0, 0, 0, 4, 0],\n",
" [0, 0, 0, 0, 5],\n",
" ]\n",
")\n",
"B = np.array(\n",
" [\n",
" [1, 2, 3, 4, 5],\n",
" [2, 3, 4, 5, 6],\n",
" [3, 4, 5, 6, 7],\n",
" [4, 5, 6, 7, 8],\n",
" [5, 6, 7, 8, 9],\n",
" ]\n",
")\n",
"d = np.diag(A)\n",
"dd = np.array([np.diag(A)])\n",
"dt = np.transpose(d)\n",
"ddt = np.transpose(dd)\n",
"Ad = np.diag(np.diag(A))\n",
"\n",
"print(np.dot(np.linalg.inv(A), A))"
]
@@ -138,11 +154,11 @@
" x = 0 * b\n",
" n = len(b)\n",
" if np.allclose(A, np.triu(A)):\n",
" for i in range(n-1, -1, -1):\n",
" x[i] = (b[i] - np.dot(A[i,i+1:], x[i+1:])) / A[i,i]\n",
" for i in range(n - 1, -1, -1):\n",
" x[i] = (b[i] - np.dot(A[i, i + 1 :], x[i + 1 :])) / A[i, i]\n",
" elif np.allclose(A, np.tril(A)):\n",
" for i in range(n):\n",
" x[i] = (b[i] - np.dot(A[i,:i], x[:i])) / A[i,i]\n",
" x[i] = (b[i] - np.dot(A[i, :i], x[:i])) / A[i, i]\n",
" else:\n",
" raise ValueError(\"A est ni triangulaire supérieure ni triangulaire inférieure\")\n",
" return x"
@@ -171,7 +187,7 @@
"b = np.dot(A, xe)\n",
"x = remontee_descente(A, b)\n",
"\n",
"print(np.dot(x - xe, x-xe))"
"print(np.dot(x - xe, x - xe))"
]
},
{
@@ -263,9 +279,9 @@
" U = A\n",
" n = len(A)\n",
" for j in range(n):\n",
" for i in range(j+1, n):\n",
" beta = U[i,j]/U[j,j]\n",
" U[i,j:] = U[i,j:] - beta * U[j, j:]\n",
" for i in range(j + 1, n):\n",
" beta = U[i, j] / U[j, j]\n",
" U[i, j:] = U[i, j:] - beta * U[j, j:]\n",
" return U"
]
},
@@ -282,14 +298,16 @@
" if n != m:\n",
" raise ValueError(\"Erreur de dimension : A doit etre carré\")\n",
" if n != b.size:\n",
" raise valueError(\"Erreur de dimension : le nombre de lignes de A doit être égal au nombr ede colonnes de b\")\n",
" U = np.zeros((n, n+1))\n",
" raise valueError(\n",
" \"Erreur de dimension : le nombre de lignes de A doit être égal au nombr ede colonnes de b\"\n",
" )\n",
" U = np.zeros((n, n + 1))\n",
" U = A\n",
" V = b\n",
" for j in range(n):\n",
" for i in range(j+1, n):\n",
" beta = U[i,j]/U[j,j]\n",
" U[i,j:] = U[i,j:] - beta * U[j, j:]\n",
" for i in range(j + 1, n):\n",
" beta = U[i, j] / U[j, j]\n",
" U[i, j:] = U[i, j:] - beta * U[j, j:]\n",
" V[i] = V[i] - beta * V[j]\n",
" return remontee_descente(U, V)"
]

141
L3/Calculs Numériques/DM1.ipynb
View File

@@ -9,9 +9,9 @@
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'retina'\n",
"\n",
"import numpy as np # pour les numpy array\n",
"import matplotlib.pyplot as plt # librairie graphique\n",
"from scipy.integrate import odeint # seulement odeint"
"import numpy as np # pour les numpy array\n",
"import matplotlib.pyplot as plt # librairie graphique\n",
"from scipy.integrate import odeint # seulement odeint"
]
},
{
@@ -103,13 +103,14 @@
"source": [
"# Initialisation des variables\n",
"T = 130\n",
"t = np.arange(0, T+1)\n",
"t = np.arange(0, T + 1)\n",
"\n",
"K = 50 # capacité d'accueil maximale du milieu\n",
"K = 50 # capacité d'accueil maximale du milieu\n",
"K_star = 50 # capacité minimale pour maintenir l'espèce\n",
"r = 0.1 # taux de corissance de la capacité d'accueil du milieu.\n",
"t_fl = 30 # la find ed la période de formation.\n",
"K0 = 1 # Valeur initiale de la capacité d'accueil du milieu.\n",
"r = 0.1 # taux de corissance de la capacité d'accueil du milieu.\n",
"t_fl = 30 # la find ed la période de formation.\n",
"K0 = 1 # Valeur initiale de la capacité d'accueil du milieu.\n",
"\n",
"\n",
"def C(t):\n",
" \"\"\"\n",
@@ -117,11 +118,12 @@
" \"\"\"\n",
" return K_star + K / (1 + (K / K0 - 1) * np.exp(-r * (t - t_fl)))\n",
"\n",
"\n",
"# On trace le graphique de la solution exacte\n",
"plt.plot(t, C(t), label=\"C(t)\")\n",
"plt.hlines(K_star, 0, T, linestyle='dotted', label=\"C = K*\", color='red')\n",
"plt.hlines(K + K_star, 0, T, linestyle='dotted', label=\"C = K + K*\", color='green')\n",
"plt.plot(t_fl, K0 + K_star, 'o', label=\"(t_fl, K0 + K*)\")\n",
"plt.hlines(K_star, 0, T, linestyle=\"dotted\", label=\"C = K*\", color=\"red\")\n",
"plt.hlines(K + K_star, 0, T, linestyle=\"dotted\", label=\"C = K + K*\", color=\"green\")\n",
"plt.plot(t_fl, K0 + K_star, \"o\", label=\"(t_fl, K0 + K*)\")\n",
"plt.legend()\n",
"plt.xlim(0, 130)\n",
"plt.suptitle(\"Courbe de la solution exacte du problème\")\n",
@@ -133,14 +135,16 @@
"N0 = 10\n",
"r_N = 0.2\n",
"\n",
"\n",
"def dN(N, t, C_sol):\n",
" \"\"\"\n",
" Fonction calculant la dérivée de la solution approchée du problème à l'instant t dépendant de N(t) et de C(t)\n",
" \"\"\"\n",
" return r_N * N * (1 - N / C_sol(t))\n",
"\n",
"\n",
"t = np.linspace(0, T, 200)\n",
"N_sol = odeint(dN, N0, t, args=(C,)) # On calcule la solution a l'aide de odeint\n",
"N_sol = odeint(dN, N0, t, args=(C,)) # On calcule la solution a l'aide de odeint\n",
"\n",
"# On trace le graphique de la solution approchée en comparaison à la solution exacte\n",
"plt.plot(t, N_sol, label=\"Solution approchée\")\n",
@@ -219,42 +223,47 @@
"T, N = 200, 100\n",
"H0, P0 = 1500, 500\n",
"\n",
"\n",
"def F(X, t, a, b, c, d, p):\n",
" \"\"\"Fonction second membre pour le système\"\"\"\n",
" x, y = X\n",
" return np.array([x * (a - p - b*y), y * (-c - p + d*x)])\n",
" return np.array([x * (a - p - b * y), y * (-c - p + d * x)])\n",
"\n",
"t = np.linspace(0, T, N+1)\n",
"sardines, requins = np.meshgrid(\n",
" np.linspace(0.1, 3000, 20),\n",
" np.linspace(0.1, 4500, 30)\n",
")\n",
"\n",
"t = np.linspace(0, T, N + 1)\n",
"sardines, requins = np.meshgrid(np.linspace(0.1, 3000, 20), np.linspace(0.1, 4500, 30))\n",
"fsardines = F((sardines, requins), t, a, b, c, d, 0)[0]\n",
"frequins = F((sardines, requins), t, a, b, c, d, 0)[1]\n",
"n_sndmb = np.sqrt(fsardines**2 + frequins**2) \n",
"n_sndmb = np.sqrt(fsardines**2 + frequins**2)\n",
"\n",
"# On crée une figure à trois graphiques\n",
"fig = plt.figure(figsize=(12, 6))\n",
"ax = fig.add_subplot(1, 2, 2) # subplot pour le champ de vecteurs et le graphe sardines vs requins\n",
"axr = fig.add_subplot(2, 2, 1) # subplot pour le graphe du nombre de requins en fonction du temps\n",
"axs = fig.add_subplot(2, 2, 3) # subplot pour le graphe du nombre de sardines en fonction du temps \n",
"ax.quiver(sardines, requins, fsardines/n_sndmb, frequins/n_sndmb)\n",
"ax = fig.add_subplot(\n",
" 1, 2, 2\n",
") # subplot pour le champ de vecteurs et le graphe sardines vs requins\n",
"axr = fig.add_subplot(\n",
" 2, 2, 1\n",
") # subplot pour le graphe du nombre de requins en fonction du temps\n",
"axs = fig.add_subplot(\n",
" 2, 2, 3\n",
") # subplot pour le graphe du nombre de sardines en fonction du temps\n",
"ax.quiver(sardines, requins, fsardines / n_sndmb, frequins / n_sndmb)\n",
"\n",
"list_p = [0, 0.02, 0.04, 0.06]\n",
"for k, pk in enumerate(list_p):\n",
" couleur = (0, k/len(list_p), 1-k/len(list_p))\n",
" couleur = (0, k / len(list_p), 1 - k / len(list_p))\n",
" X = odeint(F, np.array([H0, P0]), t, args=(a, b, c, d, pk))\n",
" \n",
" # Tracer la courbe parametrée (H(t),P(t)) \n",
"\n",
" # Tracer la courbe parametrée (H(t),P(t))\n",
" ax.plot(X[:, 0], X[:, 1], linewidth=2, color=couleur, label=f\"$p={pk}$\")\n",
" \n",
" # Tracer H en fonction du temps \n",
"\n",
" # Tracer H en fonction du temps\n",
" axs.plot(t, X[:, 0], label=f\"Sardines pour p={pk}\", color=couleur)\n",
" \n",
"\n",
" # Tracer P en fonction du temps\n",
" axr.plot(t, X[:, 1], label=f\"Requins pour p={pk}\", color=couleur)\n",
" \n",
"ax.axis('equal')\n",
"\n",
"ax.axis(\"equal\")\n",
"ax.set_title(\"Champ de vecteur du problème de Lotka-Volterra\")\n",
"ax.set_xlabel(\"Sardines\")\n",
"ax.set_ylabel(\"Requins\")\n",
@@ -266,8 +275,8 @@
"axs.set_xlabel(\"Temps t\")\n",
"axs.set_ylabel(\"Sardines\")\n",
"\n",
"axs.set_title('Evolution des sardines')\n",
"axr.set_title('Evolution des requins')\n",
"axs.set_title(\"Evolution des sardines\")\n",
"axr.set_title(\"Evolution des requins\")\n",
"plt.show()"
]
},
@@ -309,11 +318,11 @@
"source": [
"def crank_nicolson(y0, T, N, r):\n",
" \"\"\"\n",
" schéma de Crank-Nicolson pour le modèle de Malthus \n",
" \n",
" schéma de Crank-Nicolson pour le modèle de Malthus\n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" y0: float\n",
" donnée initiale\n",
" T: float\n",
@@ -325,34 +334,35 @@
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" t: ndarray\n",
" les instants où la solution approchée est calculée\n",
" y: ndarray\n",
" les valeurs de la solution approchée par le theta-schema\n",
" \"\"\"\n",
" \n",
"\n",
" dt = T / N\n",
" t = np.zeros(N+1)\n",
" y = np.zeros(N+1)\n",
" t = np.zeros(N + 1)\n",
" y = np.zeros(N + 1)\n",
" tk, yk = 0, y0\n",
" y[0] = yk\n",
" \n",
"\n",
" for n in range(N):\n",
" tk += dt\n",
" yk *= (2 + dt * r) / (2 - dt * r)\n",
" y[n+1] = yk\n",
" t[n+1] = tk\n",
" \n",
" y[n + 1] = yk\n",
" t[n + 1] = tk\n",
"\n",
" return t, y\n",
"\n",
"\n",
"def euler_explicit(y0, T, N, r):\n",
" \"\"\"\n",
" schéma de d'Euler pour le modèle de Malthus \n",
" \n",
" schéma de d'Euler pour le modèle de Malthus\n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" y0: float\n",
" donnée initiale\n",
" T: float\n",
@@ -364,29 +374,30 @@
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" t: ndarray\n",
" les instants où la solution approchée est calculée\n",
" y: ndarray\n",
" les valeurs de la solution approchée par le theta-schema\n",
" \"\"\"\n",
" dt = T / N\n",
" t = np.zeros(N+1)\n",
" y = np.zeros(N+1)\n",
" t = np.zeros(N + 1)\n",
" y = np.zeros(N + 1)\n",
" tk, yk = 0, y0\n",
" y[0] = yk\n",
" \n",
"\n",
" for n in range(N):\n",
" tk += dt\n",
" yk += dt * r * yk\n",
" y[n+1] = yk\n",
" t[n+1] = tk\n",
" \n",
" y[n + 1] = yk\n",
" t[n + 1] = tk\n",
"\n",
" return t, y\n",
"\n",
"\n",
"def solution_exacte(t):\n",
" \"\"\"\n",
" Fonction calculant la solution exacte du modèle de Malthus à l'instant t \n",
" Fonction calculant la solution exacte du modèle de Malthus à l'instant t\n",
" \"\"\"\n",
" return y0 * np.exp(r * t)"
]
@@ -436,23 +447,25 @@
"# Schéma d'Euler explicite\n",
"ax = fig.add_subplot(1, 2, 1)\n",
"for n in liste_N:\n",
" t, y = euler_explicit(y0, T, n, r) # On calcule la fonction Euler pour chaque n\n",
" t, y = euler_explicit(y0, T, n, r) # On calcule la fonction Euler pour chaque n\n",
" ax.scatter(t, y, label=f\"Solution approchée pour N={n}\")\n",
" \n",
"ax.plot(t_exact, solution_exacte(t_exact), label='Solution exacte')\n",
"\n",
"ax.plot(t_exact, solution_exacte(t_exact), label=\"Solution exacte\")\n",
"ax.legend()\n",
"ax.axis('equal')\n",
"ax.axis(\"equal\")\n",
"ax.set_title(\"Schéma d'Euler explicite\")\n",
"ax.set_xlabel(\"Temps t\")\n",
"ax.set_ylabel(\"y\")\n",
" \n",
"\n",
"\n",
"# Schéma de Crank-Nicolson\n",
"ax = fig.add_subplot(1, 2, 2)\n",
"for n in liste_N:\n",
" t, y = crank_nicolson(y0, T, n, r) # On calcule la fonction Crank-Nicolson pour chaque n\n",
" t, y = crank_nicolson(\n",
" y0, T, n, r\n",
" ) # On calcule la fonction Crank-Nicolson pour chaque n\n",
" ax.scatter(t, y, label=f\"Solution approchée pour N={n}\")\n",
"ax.plot(t_exact, solution_exacte(t_exact), label='Solution exacte')\n",
"ax.plot(t_exact, solution_exacte(t_exact), label=\"Solution exacte\")\n",
"ax.legend()\n",
"ax.set_title(\"Schéma de Crank-Nicolson\")\n",
"ax.set_xlabel(\"Temps t\")\n",
@@ -504,7 +517,7 @@
" t, sol_appr = crank_nicolson(y0, T, n, r)\n",
" sol_ex = solution_exacte(t)\n",
" erreur = np.max(np.abs(sol_appr - sol_ex))\n",
" #erreur = np.linalg.norm(sol_appr - sol_ex, np.inf)\n",
" # erreur = np.linalg.norm(sol_appr - sol_ex, np.inf)\n",
" print(f\"Delta_t = {T / N:10.3e}, e = {erreur:10.3e}\")\n",
" liste_erreur[k] = erreur\n",
"\n",
@@ -514,7 +527,7 @@
"ax.scatter(liste_delta, liste_erreur, color=\"black\")\n",
"for p in [0.5, 1, 2]:\n",
" C = liste_erreur[-1] / (liste_delta[-1] ** p)\n",
" plt.plot(liste_delta, C * liste_delta ** p, label=f\"$p={p}$\")\n",
" plt.plot(liste_delta, C * liste_delta**p, label=f\"$p={p}$\")\n",
"ax.set_title(\"Erreur du schéma de Crank-Nicolson\")\n",
"ax.set_xlabel(r\"$\\Delta t$\")\n",
"ax.set_ylabel(r\"$e(\\Delta t)$\")\n",

135
L3/Calculs Numériques/DM2.ipynb
View File

@@ -153,22 +153,27 @@
"def M(x):\n",
" \"\"\"\n",
" Retourne la matrice du système (2)\n",
" \n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" x: ndarray\n",
" vecteurs contenant les valeurs [x0, x1, ..., xN]\n",
" \n",
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" out: ndarray\n",
" matrice du système (2)\n",
" \"\"\"\n",
" h = x[1:] - x[:-1] # x[i+1] - x[i]\n",
" return np.diag(2*(1/h[:-1] + 1/h[1:])) + np.diag(1/h[1:-1], k=-1) + np.diag(1/h[1:-1], k=1)\n",
" \n",
" h = x[1:] - x[:-1] # x[i+1] - x[i]\n",
" return (\n",
" np.diag(2 * (1 / h[:-1] + 1 / h[1:]))\n",
" + np.diag(1 / h[1:-1], k=-1)\n",
" + np.diag(1 / h[1:-1], k=1)\n",
" )\n",
"\n",
"\n",
"# Test\n",
"print(M(np.array([0, 1, 2, 3, 4])))"
]
@@ -191,10 +196,10 @@
"def sprime(x, y, p0, pN):\n",
" \"\"\"\n",
" Retourne la solution du système (2)\n",
" \n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" x: ndarray\n",
" vecteurs contenant les valeurs [x0, x1, ..., xN]\n",
" y: ndarray\n",
@@ -203,18 +208,18 @@
" première valeur du vecteur p\n",
" pN: int\n",
" N-ième valeur du vecteur p\n",
" \n",
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" out: ndarray\n",
" solution du système (2)\n",
" \"\"\"\n",
" h = x[1:] - x[:-1]\n",
" delta_y = (y[1:] - y[:-1]) / h\n",
" c = 3 * (delta_y[1:]/h[1:] + delta_y[:-1]/h[:-1])\n",
" c[0] -= p0/h[0]\n",
" c[-1] -= pN/h[-1]\n",
" c = 3 * (delta_y[1:] / h[1:] + delta_y[:-1] / h[:-1])\n",
" c[0] -= p0 / h[0]\n",
" c[-1] -= pN / h[-1]\n",
" return np.linalg.solve(M(x), c)"
]
},
@@ -273,52 +278,54 @@
"def f(x):\n",
" \"\"\"\n",
" Retourne la fonction f évaluée aux points x\n",
" \n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" x: ndarray\n",
" vecteurs contenant les valeurs [x0, x1, ..., xN]\n",
" \n",
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" out: ndarray\n",
" Valeur de la fonction f aux points x\n",
" \"\"\"\n",
" return 1 / (1 + x**2)\n",
"\n",
"\n",
"def fprime(x):\n",
" \"\"\"\n",
" Retourne la fonction dérivée de f évaluée aux points x\n",
" \n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" x: ndarray\n",
" vecteurs contenant les valeurs [x0, x1, ..., xN]\n",
" \n",
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" out: ndarray\n",
" Valeur de la fonction dérivée de f aux points x\n",
" \"\"\"\n",
" return -2*x/((1+x**2)**2)\n",
" return -2 * x / ((1 + x**2) ** 2)\n",
"\n",
"# Paramètres \n",
"\n",
"# Paramètres\n",
"xx = np.linspace(-5, 5, 200)\n",
"x = np.linspace(-5, 5, 21)\n",
"pi = sprime(x, f(x), fprime(-5), fprime(5))\n",
"\n",
"# Graphique\n",
"fig, ax = plt.subplots(figsize=(6, 6))\n",
"ax.plot(xx, fprime(xx), label=f'$f\\'$', color='red')\n",
"ax.scatter(x[1:-1], pi, label=f'$p_i$')\n",
"ax.plot(xx, fprime(xx), label=\"$f'$\", color=\"red\")\n",
"ax.scatter(x[1:-1], pi, label=\"$p_i$\")\n",
"ax.legend()\n",
"ax.set_xlabel(f'$x$')\n",
"ax.set_ylabel(f'$f(x)$')\n",
"ax.set_title('Les pentes de la spline cubique')"
"ax.set_xlabel(\"$x$\")\n",
"ax.set_ylabel(\"$f(x)$\")\n",
"ax.set_title(\"Les pentes de la spline cubique\")"
]
},
{
@@ -363,10 +370,10 @@
"def splines(x, y, p0, pN):\n",
" \"\"\"\n",
" Retourne la matrice S de taille (4, N)\n",
" \n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" x: ndarray\n",
" vecteurs contenant les valeurs [x0, x1, ..., xN]\n",
" y: ndarray\n",
@@ -375,20 +382,20 @@
" première valeur du vecteur p\n",
" pN: int\n",
" N-ième valeur du vecteur p\n",
" \n",
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" out: ndarray\n",
" Matrice S de taille (4, N) tel que la i-ième ligne contient les valeurs a_i, b_i, c_i et d_i\n",
" \"\"\"\n",
" h = x[1:] - x[:-1]\n",
" delta_y = (y[1:] - y[:-1]) / h\n",
" \n",
"\n",
" a = y\n",
" b = np.concatenate((np.array([p0]), sprime(x, y, p0, pN), np.array([pN])))\n",
" c = 3/h * delta_y - (b[1:] + 2*b[:-1]) / h\n",
" d = 1/h**2 * (b[1:] + b[:-1]) - 2/h**2 * delta_y\n",
" c = 3 / h * delta_y - (b[1:] + 2 * b[:-1]) / h\n",
" d = 1 / h**2 * (b[1:] + b[:-1]) - 2 / h**2 * delta_y\n",
" return np.transpose([a[:-1], b[:-1], c, d])"
]
},
@@ -412,34 +419,38 @@
},
"outputs": [],
"source": [
"def spline_eval( x, xx, S ):\n",
"def spline_eval(x, xx, S):\n",
" \"\"\"\n",
" Evalue une spline définie par des noeuds équirepartis\n",
" \n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" x: ndarray\n",
" noeuds définissant la spline\n",
" \n",
"\n",
" xx: ndarray\n",
" abscisses des points d'évaluation\n",
" \n",
"\n",
" S: ndarray\n",
" de taille (x.size-1, 4)\n",
" tableau dont la i-ème ligne contient les coéficients du polynome cubique qui est la restriction\n",
" de la spline à l'intervalle [x_i, x_{i+1}]\n",
" \n",
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" ndarray\n",
" ordonnées des points d'évaluation\n",
" \"\"\"\n",
" ind = ( np.floor( ( xx - x[ 0 ] ) / ( x[ 1 ] - x[ 0 ] ) ) ).astype( int )\n",
" ind = np.where( ind == x.size-1, ind - 1 , ind )\n",
" yy = S[ ind, 0 ] + S[ ind, 1 ] * ( xx - x[ ind ] ) + \\\n",
" S[ ind, 2 ] * ( xx - x[ ind ] )**2 + S[ ind, 3 ] * ( xx - x[ ind ] )**3\n",
" ind = (np.floor((xx - x[0]) / (x[1] - x[0]))).astype(int)\n",
" ind = np.where(ind == x.size - 1, ind - 1, ind)\n",
" yy = (\n",
" S[ind, 0]\n",
" + S[ind, 1] * (xx - x[ind])\n",
" + S[ind, 2] * (xx - x[ind]) ** 2\n",
" + S[ind, 3] * (xx - x[ind]) ** 3\n",
" )\n",
" return yy"
]
},
@@ -472,21 +483,21 @@
}
],
"source": [
"# Paramètres \n",
"# Paramètres\n",
"x = np.linspace(-5, 5, 6)\n",
"y = np.random.rand(5+1)\n",
"y = np.random.rand(5 + 1)\n",
"xx = np.linspace(-5, 5, 200)\n",
"s = splines(x, y, 0, 0)\n",
"s_eval = spline_eval(x, xx, s)\n",
"\n",
"# Graphique\n",
"fig, ax = plt.subplots(figsize=(6, 6))\n",
"ax.plot(xx, s_eval, label='spline cubique interpolateur', color='red')\n",
"ax.scatter(x, y, label=f'$(x_i, y_i)$')\n",
"ax.plot(xx, s_eval, label=\"spline cubique interpolateur\", color=\"red\")\n",
"ax.scatter(x, y, label=\"$(x_i, y_i)$\")\n",
"ax.legend()\n",
"ax.set_xlabel(f'$x$')\n",
"ax.set_ylabel(f'$f(x)$')\n",
"ax.set_title('Evaluation de la spline cubique')"
"ax.set_xlabel(\"$x$\")\n",
"ax.set_ylabel(\"$f(x)$\")\n",
"ax.set_title(\"Evaluation de la spline cubique\")"
]
},
{
@@ -525,7 +536,7 @@
}
],
"source": [
"# Paramètres \n",
"# Paramètres\n",
"a, b = -5, 5\n",
"N_list = [4, 9, 19]\n",
"\n",
@@ -533,17 +544,17 @@
"fig, ax = plt.subplots(figsize=(15, 6))\n",
"\n",
"for N in N_list:\n",
" x = np.linspace(a, b, N+1)\n",
" x = np.linspace(a, b, N + 1)\n",
" xx = np.linspace(a, b, 200)\n",
" s = splines(x, f(x), 0, 0)\n",
" s_eval = spline_eval(x, xx, s)\n",
" ax.plot(xx, s_eval, label=f'Spline cubique interpolateur pour N={N}')\n",
" ax.scatter(x, f(x), label=f'f(x) pour N={N}')\n",
" \n",
" ax.plot(xx, s_eval, label=f\"Spline cubique interpolateur pour N={N}\")\n",
" ax.scatter(x, f(x), label=f\"f(x) pour N={N}\")\n",
"\n",
"ax.legend()\n",
"ax.set_xlabel(f'$x$')\n",
"ax.set_ylabel(f'$f(x)$')\n",
"ax.set_title('Evaluation de la spline cubique')"
"ax.set_xlabel(\"$x$\")\n",
"ax.set_ylabel(\"$f(x)$\")\n",
"ax.set_title(\"Evaluation de la spline cubique\")"
]
},
{

109
L3/Calculs Numériques/DM3.ipynb
View File

@@ -60,17 +60,20 @@
},
"outputs": [],
"source": [
"def f0(x): \n",
"def f0(x):\n",
" return np.exp(x)\n",
"\n",
"\n",
"def f1(x):\n",
" return 1 / (1 + 16*np.power(x, 2))\n",
" return 1 / (1 + 16 * np.power(x, 2))\n",
"\n",
"\n",
"def f2(x):\n",
" return np.power(np.abs(x**2 - 1/4), 3)\n",
" return np.power(np.abs(x**2 - 1 / 4), 3)\n",
"\n",
"\n",
"def f3(x):\n",
" return np.power(np.abs(x+1/2), 1/2)"
" return np.power(np.abs(x + 1 / 2), 1 / 2)"
]
},
{
@@ -176,10 +179,12 @@
],
"source": [
"for f in [f0, f1, f2, f3]:\n",
" print(f'Calcule de I(f) par la méthode de gauss et par la formule quadratique pour la fonction {f.__name__}')\n",
" print(\n",
" f\"Calcule de I(f) par la méthode de gauss et par la formule quadratique pour la fonction {f.__name__}\"\n",
" )\n",
" for n in range(1, 11):\n",
" print(f\"Pour n = {n}, gauss = {gauss(f, n)} et quad = {quad(f, -1, 1)[0]}\")\n",
" print('')"
" print(\"\")"
]
},
{
@@ -211,10 +216,10 @@
"def simpson(f, N):\n",
" if N % 2 == 0:\n",
" raise ValueError(\"N doit est impair.\")\n",
" \n",
"\n",
" h = 2 / (2 * (N - 1) // 2)\n",
" fx = f(np.linspace(-1, 1, N))\n",
" \n",
"\n",
" return (h / 3) * (fx[0] + 4 * fx[1:-1:2].sum() + 2 * fx[2:-1:2].sum() + fx[-1])"
]
},
@@ -270,10 +275,12 @@
],
"source": [
"for f in [f0, f1, f2, f3]:\n",
" print(f'Calcule de I(f) par la méthode de simpson et par la formule quadratique pour la fonction {f.__name__}')\n",
" print(\n",
" f\"Calcule de I(f) par la méthode de simpson et par la formule quadratique pour la fonction {f.__name__}\"\n",
" )\n",
" for n in range(3, 16, 2):\n",
" print(f\"Pour n = {n}, simpson = {simpson(f, n)} et quad = {quad(f, -1, 1)[0]}\")\n",
" print('')"
" print(\"\")"
]
},
{
@@ -336,7 +343,7 @@
" elif N == 1:\n",
" return x\n",
" else:\n",
" return 2 * x * poly_tchebychev(x, N-1) - poly_tchebychev(x, N-2)"
" return 2 * x * poly_tchebychev(x, N - 1) - poly_tchebychev(x, N - 2)"
]
},
{
@@ -346,7 +353,7 @@
"outputs": [],
"source": [
"def points_tchebychev(N):\n",
" k = np.arange(1, N+1)\n",
" k = np.arange(1, N + 1)\n",
" return np.cos((2 * k - 1) * np.pi / (2 * N))"
]
},
@@ -449,7 +456,7 @@
" lamk = np.zeros(N)\n",
" for k in range(N):\n",
" s = 0\n",
" for m in range(1, N//2+1):\n",
" for m in range(1, N // 2 + 1):\n",
" T = poly_tchebychev(xk[k], 2 * m)\n",
" s += 2 * T / (4 * np.power(m, 2) - 1)\n",
" lamk[k] = 2 / N * (1 - s)\n",
@@ -529,10 +536,12 @@
],
"source": [
"for f in [f0, f1, f2, f3]:\n",
" print(f'Calcule de I(f) par la méthode de fejer et par la formule quadratique pour la fonction {f.__name__}')\n",
" print(\n",
" f\"Calcule de I(f) par la méthode de fejer et par la formule quadratique pour la fonction {f.__name__}\"\n",
" )\n",
" for n in range(1, 11):\n",
" print(f\"Pour n = {n}, fejer = {fejer(f, n)} et quad = {quad(f, -1, 1)[0]}\")\n",
" print('')"
" print(\"\")"
]
},
{
@@ -568,26 +577,44 @@
"figure = plt.figure(figsize=(15, 10))\n",
"for fi, f in enumerate([f0, f1, f2, f3]):\n",
" error_gauss = np.zeros((N,))\n",
" error_simp = np.zeros(((N-1)//2,))\n",
" error_simp = np.zeros(((N - 1) // 2,))\n",
" error_fejer = np.zeros((N,))\n",
" I_quad, _ = quad(f, -1, 1)\n",
" \n",
" for n in range(1, N+1):\n",
"\n",
" for n in range(1, N + 1):\n",
" I_gauss = gauss(f, n)\n",
" error_gauss[n-1] = np.abs(I_gauss - I_quad)\n",
" error_gauss[n - 1] = np.abs(I_gauss - I_quad)\n",
" I_fejer = fejer(f, n)\n",
" error_fejer[n-1] = np.abs(I_fejer - I_quad)\n",
" \n",
" for n in range( 3, N+1, 2 ):\n",
" error_fejer[n - 1] = np.abs(I_fejer - I_quad)\n",
"\n",
" for n in range(3, N + 1, 2):\n",
" I_simp = simpson(f, n)\n",
" error_simp[(n-2)//2] = np.abs(I_simp - I_quad)\n",
" \n",
" error_simp[(n - 2) // 2] = np.abs(I_simp - I_quad)\n",
"\n",
" ax = figure.add_subplot(2, 2, fi + 1)\n",
" ax.scatter(np.arange(1, N+1), np.log10(error_gauss, out = -16. * np.ones(error_gauss.shape), \n",
" where = (error_gauss > 1e-16)), label = 'Gauss', marker=\"+\")\n",
" ax.scatter(np.arange(3, N+1, 2), np.log10( error_simp ), label = 'Simpson', marker=\"+\")\n",
" ax.scatter(np.arange(1, N+1), np.log10(error_fejer, out = -16. * np.ones(error_fejer.shape), \n",
" where = (error_fejer > 1e-16)), label = 'Fejer', marker=\"+\")\n",
" ax.scatter(\n",
" np.arange(1, N + 1),\n",
" np.log10(\n",
" error_gauss,\n",
" out=-16.0 * np.ones(error_gauss.shape),\n",
" where=(error_gauss > 1e-16),\n",
" ),\n",
" label=\"Gauss\",\n",
" marker=\"+\",\n",
" )\n",
" ax.scatter(\n",
" np.arange(3, N + 1, 2), np.log10(error_simp), label=\"Simpson\", marker=\"+\"\n",
" )\n",
" ax.scatter(\n",
" np.arange(1, N + 1),\n",
" np.log10(\n",
" error_fejer,\n",
" out=-16.0 * np.ones(error_fejer.shape),\n",
" where=(error_fejer > 1e-16),\n",
" ),\n",
" label=\"Fejer\",\n",
" marker=\"+\",\n",
" )\n",
" ax.legend()\n",
" ax.set_title(f\"Erreur de différentes méthodes de quadrature pour {f.__name__}\")\n",
" ax.set_xlabel(\"n\")\n",
@@ -672,8 +699,13 @@
"def f(x, k):\n",
" return x**k\n",
"\n",
"\n",
"print(\"-----------------------------------------------------------------------\")\n",
"print(\"{:>5s} | {:>7s} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s}\".format(\"N\", \"x^0\", \"x^2\", \"x^4\", \"x^6\", \"x^8\", \"x^10\"))\n",
"print(\n",
" \"{:>5s} | {:>7s} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s}\".format(\n",
" \"N\", \"x^0\", \"x^2\", \"x^4\", \"x^6\", \"x^8\", \"x^10\"\n",
" )\n",
")\n",
"print(\"-----------------------------------------------------------------------\")\n",
"\n",
"for N in range(1, 11):\n",
@@ -683,7 +715,10 @@
" I_exact = 2 / (k + 1) if k % 2 == 0 else 0\n",
" approx_error = np.abs(I_approx - I_exact)\n",
" approx_errors.append(approx_error)\n",
" print(\"{:5d} | \".format(N) + \" \".join(\"{:.3f} \".format(e) for e in approx_errors))"
" print(\n",
" \"{:5d} | \".format(N)\n",
" + \" \".join(\"{:.3f} \".format(e) for e in approx_errors)\n",
" )"
]
},
{
@@ -722,8 +757,13 @@
"def f(x, k):\n",
" return x**k\n",
"\n",
"\n",
"print(\"-----------------------------------------------------------------------\")\n",
"print(\"{:>5s} | {:>7s} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s}\".format(\"N\", \"x^0\", \"x^2\", \"x^4\", \"x^6\", \"x^8\", \"x^10\"))\n",
"print(\n",
" \"{:>5s} | {:>7s} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s}\".format(\n",
" \"N\", \"x^0\", \"x^2\", \"x^4\", \"x^6\", \"x^8\", \"x^10\"\n",
" )\n",
")\n",
"print(\"-----------------------------------------------------------------------\")\n",
"\n",
"for N in range(1, 11):\n",
@@ -733,7 +773,10 @@
" I_exact = 2 / (k + 1) if k % 2 == 0 else 0\n",
" approx_error = np.abs(I_approx - I_exact)\n",
" approx_errors.append(approx_error)\n",
" print(\"{:5d} | \".format(N) + \" \".join(\"{:.3f} \".format(e) for e in approx_errors))"
" print(\n",
" \"{:5d} | \".format(N)\n",
" + \" \".join(\"{:.3f} \".format(e) for e in approx_errors)\n",
" )"
]
},
{

52
L3/Calculs Numériques/Interpolation_2.ipynb
View File

@@ -22,8 +22,8 @@
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'retina'\n",
"\n",
"import numpy as np # pour les numpy array\n",
"import matplotlib.pyplot as plt # librairie graphique"
"import numpy as np # pour les numpy array\n",
"import matplotlib.pyplot as plt # librairie graphique"
]
},
{
@@ -72,11 +72,11 @@
" N = len(x)\n",
" M = len(xx)\n",
" L = np.ones((M, N))\n",
" \n",
"\n",
" for i in range(N):\n",
" for j in range(N):\n",
" if i != j:\n",
" L[:, i] *= (xx - x[j])/(x[i]-x[j])\n",
" L[:, i] *= (xx - x[j]) / (x[i] - x[j])\n",
" return L.dot(y)"
]
},
@@ -119,7 +119,7 @@
"y = np.random.rand(N)\n",
"xx = np.linspace(0, 1, 200)\n",
"\n",
"plt.scatter(x,y)\n",
"plt.scatter(x, y)\n",
"plt.plot(xx, interp_Lagrange(x, y, xx))"
]
},
@@ -155,10 +155,16 @@
"outputs": [],
"source": [
"def equirepartis(a, b, N):\n",
" return np.array([a + (b-a) * (i-1)/(N-1) for i in range(1, N+1)])\n",
" \n",
" return np.array([a + (b - a) * (i - 1) / (N - 1) for i in range(1, N + 1)])\n",
"\n",
"\n",
"def tchebychev(a, b, N):\n",
" return np.array([(a+b)/2 + (b-a)/2 * np.cos((2*i-1)/(2*N)*np.pi) for i in range(1, N+1)])"
" return np.array(\n",
" [\n",
" (a + b) / 2 + (b - a) / 2 * np.cos((2 * i - 1) / (2 * N) * np.pi)\n",
" for i in range(1, N + 1)\n",
" ]\n",
" )"
]
},
{
@@ -188,7 +194,7 @@
" L = np.ones_like(xx)\n",
" for j in range(N):\n",
" if i != j:\n",
" L *= (xx-x[j])/(x[i]-x[j]) \n",
" L *= (xx - x[j]) / (x[i] - x[j])\n",
" return L"
]
},
@@ -271,16 +277,16 @@
" ax[0].set_title(f\"Points équi-repartis (N={n})\")\n",
" xeq = equirepartis(a, b, n)\n",
" for i in range(n):\n",
" ax[0].scatter(xeq[i], 0, color='black')\n",
" ax[0].scatter(xeq[i], 1, color='black')\n",
" ax[0].scatter(xeq[i], 0, color=\"black\")\n",
" ax[0].scatter(xeq[i], 1, color=\"black\")\n",
" ax[0].plot(xx, Li(i, xeq, xx))\n",
" ax[0].grid()\n",
" \n",
"\n",
" ax[1].set_title(f\"Points de Tchebychev (N={n})\")\n",
" xchev = tchebychev(a, b, n)\n",
" for i in range(n):\n",
" ax[1].scatter(xchev[i], 0, color='black')\n",
" ax[1].scatter(xchev[i], 1, color='black')\n",
" ax[1].scatter(xchev[i], 0, color=\"black\")\n",
" ax[1].scatter(xchev[i], 1, color=\"black\")\n",
" ax[1].plot(xx, Li(i, xchev, xx))\n",
" ax[1].grid()"
]
@@ -325,20 +331,20 @@
}
],
"source": [
"f = lambda x: 1/(1+x**2)\n",
"f = lambda x: 1 / (1 + x**2)\n",
"a, b = -5, 5\n",
"xx = np.linspace(a, b, 200)\n",
"\n",
"plt.plot(xx, f(xx), label='Courbe de f')\n",
"plt.plot(xx, f(xx), label=\"Courbe de f\")\n",
"for n in [5, 10, 20]:\n",
" xeq = equirepartis(a, b, n)\n",
" for i in range(n):\n",
" plt.scatter(xeq[i], f(xeq[i]))\n",
" plt.plot(xx, Li(i, xeq, xx))\n",
" \n",
"\n",
"plt.ylim(-1, 1)\n",
"plt.legend()\n",
"plt.title('Interpolation de f avec Lagrange pour N points répartis')\n",
"plt.title(\"Interpolation de f avec Lagrange pour N points répartis\")\n",
"plt.grid()"
]
},
@@ -366,20 +372,20 @@
}
],
"source": [
"f = lambda x: 1/(1+x**2)\n",
"f = lambda x: 1 / (1 + x**2)\n",
"a, b = -5, 5\n",
"xx = np.linspace(a, b, 200)\n",
"\n",
"plt.plot(xx, f(xx), label='Courbe de f')\n",
"plt.plot(xx, f(xx), label=\"Courbe de f\")\n",
"\n",
"for n in [5, 10, 20]:\n",
" xchev = tchebychev(a, b, n)\n",
" for i in range(n):\n",
" plt.scatter(xchev[i], f(xchev[i]))\n",
" plt.plot(xx, Li(i, xchev, xx))\n",
" \n",
"\n",
"plt.legend()\n",
"plt.title('Interpolation de f avec Lagrange pour N points de Tchebychev')\n",
"plt.title(\"Interpolation de f avec Lagrange pour N points de Tchebychev\")\n",
"plt.grid()"
]
},
@@ -437,9 +443,11 @@
"source": [
"N = np.arange(5, 101, 5)\n",
"\n",
"\n",
"def n_inf(f, p):\n",
" return np.max([f, p])\n",
"\n",
"\n",
"# Norme inf en fct de N\n",
"for n in N:\n",
" xeq = equirepartis(a, b, n)\n",

30
L3/Calculs Numériques/Point_Fixe.ipynb
View File

@@ -98,16 +98,16 @@
"\n",
"x = np.linspace(-1, 1, 200)\n",
"ax = fig.add_subplot(3, 3, 1)\n",
"ax.plot(x, f1(x), label='Courbe f')\n",
"ax.plot(x, x, label=f'$y=x$')\n",
"ax.plot(x, f1(x), label=\"Courbe f\")\n",
"ax.plot(x, x, label=\"$y=x$\")\n",
"ax.scatter([i for i in x if f1(i) == i], [f1(i) for i in x if f1(i) == i])\n",
"ax.legend()\n",
"\n",
"x = np.linspace(-np.pi / 2, 5*np.pi / 2)\n",
"x = np.linspace(-np.pi / 2, 5 * np.pi / 2)\n",
"for fk, f in enumerate([f2, f3, f4, f5]):\n",
" ax = fig.add_subplot(3, 3, fk+2)\n",
" ax.plot(x, f(x), label='Courbe f')\n",
" ax.plot(x, x, label=f'$y=x$')\n",
" ax = fig.add_subplot(3, 3, fk + 2)\n",
" ax.plot(x, f(x), label=\"Courbe f\")\n",
" ax.plot(x, x, label=\"$y=x$\")\n",
" ax.scatter([i for i in x if f(i) == i], [f(i) for i in x if f(i) == i])\n",
" ax.legend()"
]
@@ -129,13 +129,13 @@
},
"outputs": [],
"source": [
"def point_fixe(f, x0, tol=1.e-6, itermax=5000):\n",
"def point_fixe(f, x0, tol=1.0e-6, itermax=5000):\n",
" \"\"\"\n",
" Recherche de point fixe : méthode brute x_{n+1} = f(x_n)\n",
" \n",
"\n",
" Parameters\n",
" ----------\n",
" \n",
"\n",
" f: function\n",
" la fonction dont on cherche le point fixe\n",
" x0: float\n",
@@ -144,10 +144,10 @@
" critère d'arrêt : |x_{n+1} - x_n| < tol\n",
" itermax: int\n",
" le nombre maximal d'itérations autorisées\n",
" \n",
"\n",
" Returns\n",
" -------\n",
" \n",
"\n",
" x: float\n",
" la valeur trouvée pour le point fixe\n",
" niter: int\n",
@@ -225,14 +225,14 @@
"source": [
"# F1\n",
"\n",
"fig, ax = plt.subplots(figsize=(6,6))\n",
"fig, ax = plt.subplots(figsize=(6, 6))\n",
"\n",
"x, niter, xL = point_fixe(f1, -0.5)\n",
"xx = np.linspace(-1, 1, 200)\n",
"\n",
"ax.plot(xx, f1(xx), label='Courbe f1')\n",
"ax.plot(xx, xx, label=f'$y=x$')\n",
"ax.scatter(x, f1(x), label='Point Fixe')\n",
"ax.plot(xx, f1(xx), label=\"Courbe f1\")\n",
"ax.plot(xx, xx, label=\"$y=x$\")\n",
"ax.scatter(x, f1(x), label=\"Point Fixe\")\n",
"ax.legend()\n",
"ax.set_title(f\"Nombre d'itérations : {niter}\")"
]

198
L3/Equations Différentielles/TP1_EDO_EulerExp.ipynb
View File

@@ -103,38 +103,45 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# Fonction f définissant l'EDO\n",
"def f1(t,y):\n",
" return -t*y\n",
"\n",
"# Solution exacte \n",
"def uex1(t,y0):\n",
" return y0 * np.exp(-t**2 /2)\n",
"# Fonction f définissant l'EDO\n",
"def f1(t, y):\n",
" return -t * y\n",
"\n",
"\n",
"# Solution exacte\n",
"def uex1(t, y0):\n",
" return y0 * np.exp(-(t**2) / 2)\n",
"\n",
"\n",
"plt.figure(1)\n",
"\n",
"## Solutions de l'EDO 1 telles que y(0)=1 et y(0)=2\n",
"\n",
"tt=np.linspace(-3, 3, 100) # vecteur représentant l'intervalle de temps\n",
"y1=uex1(tt, 1) # sol. exacte avec y_0=1\n",
"y2=uex1(tt, 2) # sol. exacte avec y_0=2\n",
"plt.plot(tt,y1,label='y(0)=1')\n",
"plt.plot(tt,y2,label='y(0)=2')\n",
"tt = np.linspace(-3, 3, 100) # vecteur représentant l'intervalle de temps\n",
"y1 = uex1(tt, 1) # sol. exacte avec y_0=1\n",
"y2 = uex1(tt, 2) # sol. exacte avec y_0=2\n",
"plt.plot(tt, y1, label=\"y(0)=1\")\n",
"plt.plot(tt, y2, label=\"y(0)=2\")\n",
"\n",
"##Tracé du champ de vecteurs\n",
"\n",
"plt.title('Solution exacte pour y0 et y1')\n",
"plt.title(\"Solution exacte pour y0 et y1\")\n",
"plt.legend()\n",
"plt.xlabel('t')\n",
"plt.ylabel('y')\n",
"plt.xlabel(\"t\")\n",
"plt.ylabel(\"y\")\n",
"\n",
"t=np.linspace(-5,5,35) # abcisse des points de la grille \n",
"y=np.linspace(0,2.1,23) # ordonnées des points de la grille \n",
"T,Y=np.meshgrid(t,y) # grille de points dans le plan (t,y) \n",
"U=np.ones(T.shape)/np.sqrt(1+f1(T,Y)**2) # matrice avec les composantes horizontales des vecteurs (1), normalisées \n",
"V=f1(T,Y)/np.sqrt(1+f1(T,Y)**2) # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées \n",
"plt.quiver(T,Y,U,V,angles='xy',scale=20,color='cyan')\n",
"plt.axis([-5,5,0,2.1])"
"t = np.linspace(-5, 5, 35) # abcisse des points de la grille\n",
"y = np.linspace(0, 2.1, 23) # ordonnées des points de la grille\n",
"T, Y = np.meshgrid(t, y) # grille de points dans le plan (t,y)\n",
"U = np.ones(T.shape) / np.sqrt(\n",
" 1 + f1(T, Y) ** 2\n",
") # matrice avec les composantes horizontales des vecteurs (1), normalisées\n",
"V = f1(T, Y) / np.sqrt(\n",
" 1 + f1(T, Y) ** 2\n",
") # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées\n",
"plt.quiver(T, Y, U, V, angles=\"xy\", scale=20, color=\"cyan\")\n",
"plt.axis([-5, 5, 0, 2.1])"
]
},
{
@@ -185,40 +192,47 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"# Fonction f définissant l'EDO\n",
"def f2(t,y):\n",
"def f2(t, y):\n",
" return t * y**2\n",
"\n",
"# Solution exacte \n",
"def uex2(t,y0):\n",
"\n",
"# Solution exacte\n",
"def uex2(t, y0):\n",
" return y0 / (1 - y0 * t**2 / 2)\n",
"\n",
"\n",
"plt.figure(1)\n",
"\n",
"## Solutions de l'EDO 1 telles que y(0)=1 et y(0)=2\n",
"\n",
"tt=np.linspace(-3, 3, 100) # vecteur représentant l'intervalle de temps\n",
"y1=uex2(tt, 1) # sol. exacte avec y_0=1\n",
"y2=uex2(tt, 2) # sol. exacte avec y_0=2\n",
"plt.plot(tt,y1,label='y(0)=1')\n",
"plt.plot(tt,y2,label='y(0)=2')\n",
"tt = np.linspace(-3, 3, 100) # vecteur représentant l'intervalle de temps\n",
"y1 = uex2(tt, 1) # sol. exacte avec y_0=1\n",
"y2 = uex2(tt, 2) # sol. exacte avec y_0=2\n",
"plt.plot(tt, y1, label=\"y(0)=1\")\n",
"plt.plot(tt, y2, label=\"y(0)=2\")\n",
"\n",
"##Tracé du champ de vecteurs\n",
"\n",
"plt.title('Solution exacte pour y0 et y1')\n",
"plt.title(\"Solution exacte pour y0 et y1\")\n",
"plt.legend()\n",
"plt.xlabel('t')\n",
"plt.ylabel('y')\n",
"plt.xlabel(\"t\")\n",
"plt.ylabel(\"y\")\n",
"\n",
"xmin, xmax = -2, 2\n",
"ymin, ymax = -4, 4\n",
"\n",
"t=np.linspace(xmin, xmax,35) # abcisse des points de la grille \n",
"y=np.linspace(ymin, ymax) # ordonnées des points de la grille \n",
"T,Y=np.meshgrid(t,y) # grille de points dans le plan (t,y) \n",
"U=np.ones(T.shape)/np.sqrt(1+f2(T,Y)**2) # matrice avec les composantes horizontales des vecteurs (1), normalisées \n",
"V=f1(T,Y)/np.sqrt(1+f2(T,Y)**2) # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées \n",
"plt.quiver(T,Y,U,V,angles='xy',scale=20,color='cyan')\n",
"t = np.linspace(xmin, xmax, 35) # abcisse des points de la grille\n",
"y = np.linspace(ymin, ymax) # ordonnées des points de la grille\n",
"T, Y = np.meshgrid(t, y) # grille de points dans le plan (t,y)\n",
"U = np.ones(T.shape) / np.sqrt(\n",
" 1 + f2(T, Y) ** 2\n",
") # matrice avec les composantes horizontales des vecteurs (1), normalisées\n",
"V = f1(T, Y) / np.sqrt(\n",
" 1 + f2(T, Y) ** 2\n",
") # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées\n",
"plt.quiver(T, Y, U, V, angles=\"xy\", scale=20, color=\"cyan\")\n",
"plt.axis([xmin, xmax, ymin, ymax])"
]
},
@@ -339,33 +353,34 @@
],
"source": [
"# Euler explicite\n",
"def euler_exp(t0,T,y0,h,f):\n",
" t = np.arange(t0, t0+T+h, h)\n",
"def euler_exp(t0, T, y0, h, f):\n",
" t = np.arange(t0, t0 + T + h, h)\n",
" y = np.empty(t.size)\n",
" y[0] = y0\n",
" N = len(t)-1\n",
" N = len(t) - 1\n",
" for n in range(N):\n",
" y[n+1] = y[n] + h * f(t[n], y[n])\n",
" y[n + 1] = y[n] + h * f(t[n], y[n])\n",
" return t, y\n",
"\n",
"t0=0\n",
"T=1\n",
"y0=1\n",
"\n",
"t0 = 0\n",
"T = 1\n",
"y0 = 1\n",
"\n",
"for f, uex in zip([f1, f2], [uex1, uex2]):\n",
" plt.figure()\n",
" t = np.arange(0, 1, 1e-3)\n",
" y = uex(t, y0)\n",
" plt.plot(t, y, label='Solution exacte')\n",
" plt.plot(t, y, label=\"Solution exacte\")\n",
"\n",
" for h in [1/5, 1/10, 1/50]:\n",
" for h in [1 / 5, 1 / 10, 1 / 50]:\n",
" tt, y = euler_exp(t0, T, y0, h, f)\n",
" plt.plot(tt, y, label=f'Solution approchée pour h={h}')\n",
" \n",
" plt.plot(tt, y, label=f\"Solution approchée pour h={h}\")\n",
"\n",
" plt.title(f\"Solutions exacte et approchées de la fonction {f.__name__}\")\n",
" plt.legend()\n",
" plt.xlabel('t')\n",
" plt.ylabel('y')"
" plt.xlabel(\"t\")\n",
" plt.ylabel(\"y\")"
]
},
{
@@ -441,12 +456,14 @@
],
"source": [
"# Modèle de Verhulst\n",
"n=1\n",
"d=0.75\n",
"K=200\n",
"n = 1\n",
"d = 0.75\n",
"K = 200\n",
"\n",
"\n",
"def fV(t, P):\n",
" return (n - d) * P * (1 - P / K)\n",
"\n",
"def fV(t,P):\n",
" return (n-d)*P*(1-P/K)\n",
"\n",
"plt.figure()\n",
"\n",
@@ -455,20 +472,24 @@
"t0, T = 0, 50\n",
"\n",
"\n",
"for P in range(1, K+100, 15):\n",
"for P in range(1, K + 100, 15):\n",
" tt, yy = euler_exp(t0, T, P, h, fV)\n",
" plt.plot(tt, yy, label=f'Solution approchée pour h={h}')\n",
" \n",
"plt.title(f\"Solution approchée du modèle de Verhulst pour la population d'individus\")\n",
"plt.xlabel('t')\n",
"plt.ylabel('y')\n",
" plt.plot(tt, yy, label=f\"Solution approchée pour h={h}\")\n",
"\n",
"t=np.linspace(t0, T, 35) # abcisse des points de la grille \n",
"y=np.linspace(0, K+100, 23) # ordonnées des points de la grille \n",
"T,P=np.meshgrid(t,y) # grille de points dans le plan (t,y) \n",
"U=np.ones(T.shape)/np.sqrt(1+fV(T,P)**2) # matrice avec les composantes horizontales des vecteurs (1), normalisées \n",
"V=fV(T,P)/np.sqrt(1+fV(T,P)**2) # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées \n",
"plt.quiver(T,P,U,V,angles='xy',scale=20,color='cyan')\n",
"plt.title(\"Solution approchée du modèle de Verhulst pour la population d'individus\")\n",
"plt.xlabel(\"t\")\n",
"plt.ylabel(\"y\")\n",
"\n",
"t = np.linspace(t0, T, 35) # abcisse des points de la grille\n",
"y = np.linspace(0, K + 100, 23) # ordonnées des points de la grille\n",
"T, P = np.meshgrid(t, y) # grille de points dans le plan (t,y)\n",
"U = np.ones(T.shape) / np.sqrt(\n",
" 1 + fV(T, P) ** 2\n",
") # matrice avec les composantes horizontales des vecteurs (1), normalisées\n",
"V = fV(T, P) / np.sqrt(\n",
" 1 + fV(T, P) ** 2\n",
") # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées\n",
"plt.quiver(T, P, U, V, angles=\"xy\", scale=20, color=\"cyan\")\n",
"plt.legend(fontsize=4)"
]
},
@@ -527,29 +548,34 @@
}
],
"source": [
"def fS(t,P):\n",
" return (2-np.cos(t))*P-(P**2)/2-1\n",
"def fS(t, P):\n",
" return (2 - np.cos(t)) * P - (P**2) / 2 - 1\n",
"\n",
"P0=5\n",
"t0=0\n",
"T=10\n",
"h=0.1\n",
"\n",
"P0 = 5\n",
"t0 = 0\n",
"T = 10\n",
"h = 0.1\n",
"\n",
"plt.figure()\n",
"tt, yy = euler_exp(t0, T, P0, h, fS)\n",
"plt.plot(tt, yy, label=f'Solution approchée pour h={h}')\n",
" \n",
"plt.title(f\"Solutions approchée du modèle de Verhulst pour une population de saumons\")\n",
"plt.legend()\n",
"plt.xlabel('t')\n",
"plt.ylabel('y')\n",
"plt.plot(tt, yy, label=f\"Solution approchée pour h={h}\")\n",
"\n",
"t=np.linspace(t0, T, 35) # abcisse des points de la grille \n",
"y=np.linspace(0, 6, 23) # ordonnées des points de la grille \n",
"T,P=np.meshgrid(t,y) # grille de points dans le plan (t,y) \n",
"U=np.ones(T.shape)/np.sqrt(1+fS(T,P)**2) # matrice avec les composantes horizontales des vecteurs (1), normalisées \n",
"V=fS(T,P)/np.sqrt(1+fS(T,P)**2) # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées \n",
"plt.quiver(T,P,U,V,angles='xy',scale=20,color='cyan')\n"
"plt.title(\"Solutions approchée du modèle de Verhulst pour une population de saumons\")\n",
"plt.legend()\n",
"plt.xlabel(\"t\")\n",
"plt.ylabel(\"y\")\n",
"\n",
"t = np.linspace(t0, T, 35) # abcisse des points de la grille\n",
"y = np.linspace(0, 6, 23) # ordonnées des points de la grille\n",
"T, P = np.meshgrid(t, y) # grille de points dans le plan (t,y)\n",
"U = np.ones(T.shape) / np.sqrt(\n",
" 1 + fS(T, P) ** 2\n",
") # matrice avec les composantes horizontales des vecteurs (1), normalisées\n",
"V = fS(T, P) / np.sqrt(\n",
" 1 + fS(T, P) ** 2\n",
") # matrice avec les composantes verticales des vecteurs (f(t,y)), normalisées\n",
"plt.quiver(T, P, U, V, angles=\"xy\", scale=20, color=\"cyan\")"
]
},
{

237
L3/Equations Différentielles/TP2_Lokta_Volterra.ipynb
View File

@@ -95,53 +95,63 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"# Fonction F définissant l'EDO\n",
"def F(Y):\n",
" x=Y[0]\n",
" y=Y[1]\n",
" A=np.array([[0,1],[-2,-3]])\n",
" return np.dot(A,Y)\n",
" x = Y[0]\n",
" y = Y[1]\n",
" A = np.array([[0, 1], [-2, -3]])\n",
" return np.dot(A, Y)\n",
"\n",
"# ou \n",
"def F1(x,y):\n",
"\n",
"# ou\n",
"def F1(x, y):\n",
" return y\n",
"\n",
"def F2(x,y):\n",
" return -2*x-3*y\n",
"\n",
"def F2(x, y):\n",
" return -2 * x - 3 * y\n",
"\n",
"\n",
"# Solution exacte de Y'=AY, Y(t_0)=Y_0\n",
"def uex(t,t0,Y0):\n",
" U1=np.array([1,-1])\n",
" U2=np.array([1,-2])\n",
" P=np.ones((2,2))\n",
" P[:,0]=U1\n",
" P[:,1]=U2\n",
" C=np.linalg.solve(P,Y0)\n",
" return np.array([(C[0]*np.exp(-(t-t0))*U1[0]+C[1]*np.exp(-2*(t-t0))*U2[0]),(C[0]*np.exp(-(t-t0))*U1[1]+C[1]*np.exp(-2*(t-t0))*U2[1])])\n",
"def uex(t, t0, Y0):\n",
" U1 = np.array([1, -1])\n",
" U2 = np.array([1, -2])\n",
" P = np.ones((2, 2))\n",
" P[:, 0] = U1\n",
" P[:, 1] = U2\n",
" C = np.linalg.solve(P, Y0)\n",
" return np.array(\n",
" [\n",
" (C[0] * np.exp(-(t - t0)) * U1[0] + C[1] * np.exp(-2 * (t - t0)) * U2[0]),\n",
" (C[0] * np.exp(-(t - t0)) * U1[1] + C[1] * np.exp(-2 * (t - t0)) * U2[1]),\n",
" ]\n",
" )\n",
"\n",
"## Représentation des solutions pour chaque valeur de la donnée initiale \n",
"\n",
"## Représentation des solutions pour chaque valeur de la donnée initiale\n",
"tt = np.linspace(-10, 10, 100)\n",
"t0 = tt[0]\n",
"for x, y in zip([1, -2, 0, 1, 3], [2, -2, -4, -2, 4]):\n",
" sol = uex(tt, t0, [x, y])\n",
" plt.plot(sol[0], sol[1], label=f'$((x0, y0) = ({x}, {y})$')\n",
" plt.plot(sol[0], sol[1], label=f\"$((x0, y0) = ({x}, {y})$\")\n",
" plt.scatter(x, y)\n",
"\n",
"#Tracé du champ de vecteurs\n",
"x=np.linspace(-5,5,26)\n",
"y=x\n",
"xx,yy=np.meshgrid(x,y)\n",
"U=F1(xx,yy)/np.sqrt(F1(xx,yy)**2+F2(xx,yy)**2)\n",
"V=F2(xx,yy)/np.sqrt(F1(xx,yy)**2+F2(xx,yy)**2)\n",
"plt.quiver(xx,yy,U,V,angles='xy', scale=20, color='gray')\n",
"plt.axis([-5.,5,-5,5])\n",
"# Tracé du champ de vecteurs\n",
"x = np.linspace(-5, 5, 26)\n",
"y = x\n",
"xx, yy = np.meshgrid(x, y)\n",
"U = F1(xx, yy) / np.sqrt(F1(xx, yy) ** 2 + F2(xx, yy) ** 2)\n",
"V = F2(xx, yy) / np.sqrt(F1(xx, yy) ** 2 + F2(xx, yy) ** 2)\n",
"plt.quiver(xx, yy, U, V, angles=\"xy\", scale=20, color=\"gray\")\n",
"plt.axis([-5.0, 5, -5, 5])\n",
"\n",
"## Représentation des espaces propres \n",
"plt.plot(tt, -tt, label=f'SEP associé à -1', linewidth=3)\n",
"plt.plot(tt, -2*tt, label=f'SEP associé à -2', linewidth=3)\n",
"## Représentation des espaces propres\n",
"plt.plot(tt, -tt, label=\"SEP associé à -1\", linewidth=3)\n",
"plt.plot(tt, -2 * tt, label=\"SEP associé à -2\", linewidth=3)\n",
"\n",
"plt.legend(fontsize=7)\n",
"plt.title('Représentation des solutions et champs de vecteurs pour le système (S)')"
"plt.title(\"Représentation des solutions et champs de vecteurs pour le système (S)\")"
]
},
{
@@ -210,13 +220,15 @@
"source": [
"a, b, c, d = 0.1, 5e-5, 0.04, 5e-5\n",
"H0, P0 = 2000, 1000\n",
"He, Pe = c/d, a/b\n",
"He, Pe = c / d, a / b\n",
"\n",
"\n",
"def F1(H, P):\n",
" return H * (a - b*P)\n",
" return H * (a - b * P)\n",
"\n",
"\n",
"def F2(H, P):\n",
" return P * (-c + d*H)"
" return P * (-c + d * H)"
]
},
{
@@ -256,17 +268,17 @@
}
],
"source": [
"xx, yy = np.linspace(0, 3000, 20), np.linspace (0, 4500, 30)\n",
"h,p = np.meshgrid(xx, yy)\n",
"n=np.sqrt(F1(h,p)**2+F2(h,p)**2)\n",
"plt.quiver(h, p, F1(h,p)/n, F2(h,p)/n, angles='xy', scale=20, color='gray')\n",
"xx, yy = np.linspace(0, 3000, 20), np.linspace(0, 4500, 30)\n",
"h, p = np.meshgrid(xx, yy)\n",
"n = np.sqrt(F1(h, p) ** 2 + F2(h, p) ** 2)\n",
"plt.quiver(h, p, F1(h, p) / n, F2(h, p) / n, angles=\"xy\", scale=20, color=\"gray\")\n",
"\n",
"plt.vlines(He, 0, 4500, label=f'H=He={He}')\n",
"plt.hlines(Pe, 0, 3000, label=f'P=Pe={Pe}')\n",
"plt.scatter(He, Pe, label=f'(H0, P0) = (He, Pe)')\n",
"plt.scatter(H0, P0, label=f'(H, P)=(H0, P0)=({H0},{P0})', color='red')\n",
"plt.vlines(He, 0, 4500, label=f\"H=He={He}\")\n",
"plt.hlines(Pe, 0, 3000, label=f\"P=Pe={Pe}\")\n",
"plt.scatter(He, Pe, label=\"(H0, P0) = (He, Pe)\")\n",
"plt.scatter(H0, P0, label=f\"(H, P)=(H0, P0)=({H0},{P0})\", color=\"red\")\n",
"\n",
"plt.title('Le modèle de Lotka-Volterra')\n",
"plt.title(\"Le modèle de Lotka-Volterra\")\n",
"plt.legend(fontsize=7)"
]
},
@@ -374,19 +386,20 @@
"outputs": [],
"source": [
"a, b, c, d = 0.1, 5e-5, 0.04, 5e-5\n",
"T=200\n",
"T = 200\n",
"H0, P0 = 2000, 1000\n",
"He, Pe = c/d, a/b\n",
"p=0.02\n",
"He, Pe = c / d, a / b\n",
"p = 0.02\n",
"\n",
"\n",
"def voltEE(T, X0, h):\n",
" t = np.arange(0, T+h, h)\n",
" t = np.arange(0, T + h, h)\n",
" H = 0 * t\n",
" P = 0 * t\n",
" H[0], P[0] = X0\n",
" for n in range(len(t)-1):\n",
" H[n+1] = H[n] + h * F1(H[n], P[n])\n",
" P[n+1] = P[n] + h * F2(H[n], P[n])\n",
" for n in range(len(t) - 1):\n",
" H[n + 1] = H[n] + h * F1(H[n], P[n])\n",
" P[n + 1] = P[n] + h * F2(H[n], P[n])\n",
" return np.array([t, H, P])"
]
},
@@ -438,22 +451,22 @@
],
"source": [
"t, H, P = voltEE(T, [H0, P0], 0.001)\n",
"plt.plot(t, H, label='Population de sardines')\n",
"plt.plot(t, P, label='Population de requins')\n",
"plt.plot(t, H, label=\"Population de sardines\")\n",
"plt.plot(t, P, label=\"Population de requins\")\n",
"plt.legend(fontsize=7)\n",
"\n",
"plt.figure()\n",
"xx, yy = np.linspace(0, 3000, 20), np.linspace (0, 4500, 30)\n",
"h,p = np.meshgrid(xx, yy)\n",
"n=np.sqrt(F1(h,p)**2 + F2(h,p)**2)\n",
"plt.quiver (h, p, F1(h,p)/n, F2(h,p)/n, angles='xy', scale=20, color='gray')\n",
"xx, yy = np.linspace(0, 3000, 20), np.linspace(0, 4500, 30)\n",
"h, p = np.meshgrid(xx, yy)\n",
"n = np.sqrt(F1(h, p) ** 2 + F2(h, p) ** 2)\n",
"plt.quiver(h, p, F1(h, p) / n, F2(h, p) / n, angles=\"xy\", scale=20, color=\"gray\")\n",
"\n",
"plt.vlines(He, 0, 4500, label=f'H=He={He}')\n",
"plt.hlines(Pe, 0, 3000, label=f'P=Pe={Pe}')\n",
"plt.scatter(He, Pe, label=f'(H0, P0) = (He, Pe)')\n",
"plt.scatter(H0, P0, label=f'(H, P)=(H0, P0)=({H0},{P0})', color='red')\n",
"plt.vlines(He, 0, 4500, label=f\"H=He={He}\")\n",
"plt.hlines(Pe, 0, 3000, label=f\"P=Pe={Pe}\")\n",
"plt.scatter(He, Pe, label=\"(H0, P0) = (He, Pe)\")\n",
"plt.scatter(H0, P0, label=f\"(H, P)=(H0, P0)=({H0},{P0})\", color=\"red\")\n",
"\n",
"plt.title('Le modèle de Lotka-Volterra')\n",
"plt.title(\"Le modèle de Lotka-Volterra\")\n",
"plt.legend(fontsize=7)\n",
"plt.plot(H, P)"
]
@@ -467,25 +480,28 @@
"outputs": [],
"source": [
"a, b, c, d = 0.1, 5e-5, 0.04, 5e-5\n",
"T=200\n",
"T = 200\n",
"H0, P0 = 2000, 1000\n",
"He, Pe = c/d, a/b\n",
"p=0.02\n",
"He, Pe = c / d, a / b\n",
"p = 0.02\n",
"\n",
"\n",
"def F1_p(H, P):\n",
" return (a - p) * H - b * H * P\n",
"\n",
"\n",
"def F2_p(H, P):\n",
" return (-c-p) * P + d * H * P\n",
" return (-c - p) * P + d * H * P\n",
"\n",
"\n",
"def voltEE_p(T, X0, h):\n",
" t = np.arange(0, T+h, h)\n",
" t = np.arange(0, T + h, h)\n",
" H = 0 * t\n",
" P = 0 * t\n",
" H[0], P[0] = X0\n",
" for n in range(len(t)-1):\n",
" H[n+1] = H[n] + h * F1_p(H[n], P[n])\n",
" P[n+1] = P[n] + h * F2_p(H[n], P[n])\n",
" for n in range(len(t) - 1):\n",
" H[n + 1] = H[n] + h * F1_p(H[n], P[n])\n",
" P[n + 1] = P[n] + h * F2_p(H[n], P[n])\n",
" return np.array([t, H, P])"
]
},
@@ -535,22 +551,22 @@
],
"source": [
"t, H, P = voltEE_p(T, [H0, P0], 0.001)\n",
"plt.plot(t, H, label='Population de sardines')\n",
"plt.plot(t, P, label='Population de requins')\n",
"plt.plot(t, H, label=\"Population de sardines\")\n",
"plt.plot(t, P, label=\"Population de requins\")\n",
"plt.legend(fontsize=7)\n",
"\n",
"plt.figure()\n",
"xx, yy = np.linspace(0, 3000, 20), np.linspace (0, 4500, 30)\n",
"h,p = np.meshgrid(xx, yy)\n",
"n=np.sqrt(F1(h,p)**2 + F2(h,p)**2)\n",
"plt.quiver (h, p, F1(h,p)/n, F2(h,p)/n, angles='xy', scale=20, color='gray')\n",
"xx, yy = np.linspace(0, 3000, 20), np.linspace(0, 4500, 30)\n",
"h, p = np.meshgrid(xx, yy)\n",
"n = np.sqrt(F1(h, p) ** 2 + F2(h, p) ** 2)\n",
"plt.quiver(h, p, F1(h, p) / n, F2(h, p) / n, angles=\"xy\", scale=20, color=\"gray\")\n",
"\n",
"plt.vlines(He, 0, 4500, label=f'H=He={He}')\n",
"plt.hlines(Pe, 0, 3000, label=f'P=Pe={Pe}')\n",
"plt.scatter(He, Pe, label=f'(H0, P0) = (He, Pe)')\n",
"plt.scatter(H0, P0, label=f'(H, P)=(H0, P0)=({H0},{P0})', color='red')\n",
"plt.vlines(He, 0, 4500, label=f\"H=He={He}\")\n",
"plt.hlines(Pe, 0, 3000, label=f\"P=Pe={Pe}\")\n",
"plt.scatter(He, Pe, label=\"(H0, P0) = (He, Pe)\")\n",
"plt.scatter(H0, P0, label=f\"(H, P)=(H0, P0)=({H0},{P0})\", color=\"red\")\n",
"\n",
"plt.title('Le modèle de Lotka-Volterra')\n",
"plt.title(\"Le modèle de Lotka-Volterra\")\n",
"plt.legend(fontsize=7)\n",
"plt.plot(H, P)"
]
@@ -697,46 +713,49 @@
}
],
"source": [
"gamma=0.5\n",
"a=0.004\n",
"lm=10\n",
"gamma = 0.5\n",
"a = 0.004\n",
"lm = 10\n",
"\n",
"\n",
"def fphi(d):\n",
" return (gamma/d)*np.log(d/a)\n",
" return (gamma / d) * np.log(d / a)\n",
"\n",
"\n",
"def F(X):\n",
" (n,m)=np.shape(X)\n",
" Y=np.zeros((n,m))\n",
" Y[0,:]=X[1,:]\n",
" Xaux=np.zeros(m+2)\n",
" Xaux[-1]=1\n",
" Xaux[1:-1]=X[0,:]\n",
" Y[1,:]=fphi(Xaux[2:]-Xaux[1:-1])-fphi(Xaux[1:-1]-Xaux[0:-2])-lm*X[1,:]\n",
" (n, m) = np.shape(X)\n",
" Y = np.zeros((n, m))\n",
" Y[0, :] = X[1, :]\n",
" Xaux = np.zeros(m + 2)\n",
" Xaux[-1] = 1\n",
" Xaux[1:-1] = X[0, :]\n",
" Y[1, :] = fphi(Xaux[2:] - Xaux[1:-1]) - fphi(Xaux[1:-1] - Xaux[0:-2]) - lm * X[1, :]\n",
" return Y\n",
"\n",
"h=0.0002\n",
"T=15\n",
"N=100\n",
"\n",
"t=np.arange(0,T+h,h)\n",
"Nt=np.size(t)\n",
"X=np.zeros((2,N,Nt))\n",
"R0=-1+2*np.random.rand(N)\n",
"X0=np.arange(1/(N+1),1,1/(N+1))+0.1*R0*(1/(N+1))\n",
"h = 0.0002\n",
"T = 15\n",
"N = 100\n",
"\n",
"X[0,:,0]=X0\n",
"X[1,:,0]=X0\n",
"t = np.arange(0, T + h, h)\n",
"Nt = np.size(t)\n",
"X = np.zeros((2, N, Nt))\n",
"R0 = -1 + 2 * np.random.rand(N)\n",
"X0 = np.arange(1 / (N + 1), 1, 1 / (N + 1)) + 0.1 * R0 * (1 / (N + 1))\n",
"\n",
"X[0, :, 0] = X0\n",
"X[1, :, 0] = X0\n",
"\n",
"plt.figure(1, figsize=(24, 18))\n",
"for n in range(Nt - 1):\n",
" Y = F(X[:, :, n])\n",
" X[:, :, n + 1] = X[:, :, n] + (h / 2) * Y + (h / 2) * F(X[:, :, n] + h * Y)\n",
"\n",
"plt.figure(1,figsize=(24,18))\n",
"for n in range(Nt-1):\n",
" Y=F(X[:,:,n])\n",
" X[:,:,n+1]=X[:,:,n]+(h/2)*Y+(h/2)*F(X[:,:,n]+h*Y)\n",
" \n",
"for i in range(N):\n",
" plt.plot(t,X[0,i,:],'k')\n",
"plt.xlabel('t')\n",
"plt.ylabel('$x_i$')\n",
"plt.title('position x_i des globules rouges au cours du temps')\n"
" plt.plot(t, X[0, i, :], \"k\")\n",
"plt.xlabel(\"t\")\n",
"plt.ylabel(\"$x_i$\")\n",
"plt.title(\"position x_i des globules rouges au cours du temps\")"
]
},
{

242
L3/Equations Différentielles/TP3_Convergence.ipynb
View File

@@ -76,14 +76,15 @@
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"## Question 1\n",
"def mon_schema(t0, T, y0, h, f):\n",
" t = np.arange(t0, t0 + T + h, h)\n",
" y = 0 * t\n",
" y[0] = y0\n",
" for n in range(len(t)-1):\n",
" yn1 = y[n] + h * f(t[n] + h / 2, y[n] + h/2 * f(t[n], y[n]))\n",
" y[n+1] = y[n] + h / 2 * (f(t[n], y[n]) + f(t[n+1], yn1))\n",
" for n in range(len(t) - 1):\n",
" yn1 = y[n] + h * f(t[n] + h / 2, y[n] + h / 2 * f(t[n], y[n]))\n",
" y[n + 1] = y[n] + h / 2 * (f(t[n], y[n]) + f(t[n + 1], yn1))\n",
" return t, y"
]
},
@@ -140,25 +141,27 @@
"## Question 2\n",
"\n",
"# f second membre de l'EDO\n",
"def f(t,y):\n",
"def f(t, y):\n",
" return y + np.sin(t) * np.power(y, 2)\n",
"\n",
"\n",
"# sol. exacte de (P)\n",
"def yex(t):\n",
" return 1 / (1/2 * np.exp(-t) + (np.cos(t)-np.sin(t)) / 2)\n",
" return 1 / (1 / 2 * np.exp(-t) + (np.cos(t) - np.sin(t)) / 2)\n",
"\n",
"\n",
"t0, T = 0, 0.5\n",
"N = 100\n",
"y0 = 1\n",
"t, y_app = mon_schema(t0, T, y0, T/N, f)\n",
"t, y_app = mon_schema(t0, T, y0, T / N, f)\n",
"y_ex = yex(t)\n",
"\n",
"plt.figure(1)\n",
"plt.plot(t, y_ex, label='Solution exacte', lw=3)\n",
"plt.plot(t, y_app, '+ ', label=f'Solution approchée pour N={N}')\n",
"plt.title('Solutions approchées pour mon_schema')\n",
"plt.xlabel(f'$t$')\n",
"plt.ylabel(f'$y$')\n",
"plt.plot(t, y_ex, label=\"Solution exacte\", lw=3)\n",
"plt.plot(t, y_app, \"+ \", label=f\"Solution approchée pour N={N}\")\n",
"plt.title(\"Solutions approchées pour mon_schema\")\n",
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(\"$y$\")\n",
"plt.legend(fontsize=7)"
]
},
@@ -212,15 +215,17 @@
"# Question 3\n",
"plt.figure(2)\n",
"\n",
"for s in range(2,11):\n",
" h=(1/2)/(2**s)\n",
"for s in range(2, 11):\n",
" h = (1 / 2) / (2**s)\n",
" t, y_app = mon_schema(t0, T, y0, h, f)\n",
" plt.plot(t, y_app, label=f'h={h}')\n",
" plt.plot(t, y_app, label=f\"h={h}\")\n",
"\n",
"plt.legend(fontsize=7)\n",
"plt.xlabel(f'$t$')\n",
"plt.ylabel('$y^n$')\n",
"plt.title('Solutions approchées de (P) obtenus avec mon_schema pour différentes valeurs du pas h')"
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(\"$y^n$\")\n",
"plt.title(\n",
" \"Solutions approchées de (P) obtenus avec mon_schema pour différentes valeurs du pas h\"\n",
")"
]
},
{
@@ -252,23 +257,25 @@
}
],
"source": [
"E=[]\n",
"H=[]\n",
"E = []\n",
"H = []\n",
"\n",
"plt.figure(3)\n",
"for s in range(2,11):\n",
" h=(1/2)/(2**s)\n",
"for s in range(2, 11):\n",
" h = (1 / 2) / (2**s)\n",
" t, y_app = mon_schema(t0, T, y0, h, f)\n",
" y_ex = yex(t)\n",
" err = np.max(np.abs(y_app - y_ex))\n",
" plt.plot(t, np.abs(y_app - y_ex), label=f'h={h}')\n",
" plt.plot(t, np.abs(y_app - y_ex), label=f\"h={h}\")\n",
" E.append(err)\n",
" H.append(h)\n",
" \n",
"\n",
"plt.legend()\n",
"plt.xlabel(f'$t$')\n",
"plt.ylabel('$|y(t_n) - y^n|$')\n",
"plt.title('différence en valeur absolue entre sol. exacte et sol. approchée par mon_schema, pour différentes valeurs du pas h')"
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(\"$|y(t_n) - y^n|$\")\n",
"plt.title(\n",
" \"différence en valeur absolue entre sol. exacte et sol. approchée par mon_schema, pour différentes valeurs du pas h\"\n",
")"
]
},
{
@@ -304,13 +311,13 @@
"\n",
"H, E = np.array(H), np.array(E)\n",
"\n",
"plt.plot(T/H, E/H, label=f'$E_h / h$')\n",
"plt.plot(T/H, E/H**2, label=f'$E_h / h^2$')\n",
"plt.plot(T / H, E / H, label=\"$E_h / h$\")\n",
"plt.plot(T / H, E / H**2, label=\"$E_h / h^2$\")\n",
"\n",
"plt.legend()\n",
"plt.xlabel('N')\n",
"plt.ylabel('Erreur globale')\n",
"plt.title('Erreur globale en fonction de N')"
"plt.xlabel(\"N\")\n",
"plt.ylabel(\"Erreur globale\")\n",
"plt.title(\"Erreur globale en fonction de N\")"
]
},
{
@@ -344,14 +351,16 @@
"source": [
"plt.figure(5)\n",
"\n",
"plt.plot(np.log(H), np.log(E), '+-', label='erreur')\n",
"plt.plot(np.log(H), np.log(H), '-', label='droite pente 1')\n",
"plt.plot(np.log(H), 2*np.log(H), '-', label='droite pente 2')\n",
"plt.plot(np.log(H), np.log(E), \"+-\", label=\"erreur\")\n",
"plt.plot(np.log(H), np.log(H), \"-\", label=\"droite pente 1\")\n",
"plt.plot(np.log(H), 2 * np.log(H), \"-\", label=\"droite pente 2\")\n",
"\n",
"plt.legend()\n",
"plt.title('Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)')\n",
"plt.xlabel(f'$log(h)$')\n",
"plt.ylabel(f'$log(E)$')"
"plt.title(\n",
" \"Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)\"\n",
")\n",
"plt.xlabel(\"$log(h)$\")\n",
"plt.ylabel(\"$log(E)$\")"
]
},
{
@@ -444,24 +453,26 @@
" t = np.arange(t0, t0 + T + h, h)\n",
" y = 0 * t\n",
" y[0] = y0\n",
" for n in range(len(t)-1):\n",
" y[n+1] = y[n] + h * f(t[n+1], y[n+1])\n",
" for n in range(len(t) - 1):\n",
" y[n + 1] = y[n] + h * f(t[n + 1], y[n + 1])\n",
" return t, y\n",
"\n",
"\n",
"def crank(t0, T, y0, h, f):\n",
" t = np.arange(t0, t0 + T + h, h)\n",
" y = 0 * t\n",
" y[0] = y0\n",
" for n in range(len(t)-1):\n",
" y[n+1] = y[n] + h/2 * (f(t[n], y[n]) + f(t[n+1], y[n+1]))\n",
" for n in range(len(t) - 1):\n",
" y[n + 1] = y[n] + h / 2 * (f(t[n], y[n]) + f(t[n + 1], y[n + 1]))\n",
" return t, y\n",
"\n",
"\n",
"def adams(t0, T, y0, h, f):\n",
" t = np.arange(t0, t0 + T + h, h)\n",
" y = 0 * t\n",
" y[0] = y0\n",
" for n in range(1, len(t)-1):\n",
" y[n+1] = y[n] + h/2 * (3* f(t[n], y[n]) - f(t[n-1], y[n-1]))\n",
" for n in range(1, len(t) - 1):\n",
" y[n + 1] = y[n] + h / 2 * (3 * f(t[n], y[n]) - f(t[n - 1], y[n - 1]))\n",
" return t, y"
]
},
@@ -474,39 +485,41 @@
"outputs": [],
"source": [
"def euler_explicite(t0, T, y0, h, f):\n",
" N = int(T/h)\n",
" N = int(T / h)\n",
" n = len(y0)\n",
" t = np.linspace(t0, t0 + T, N+1)\n",
" y = np.zeros((N+1, n))\n",
" t = np.linspace(t0, t0 + T, N + 1)\n",
" y = np.zeros((N + 1, n))\n",
" y[0,] = y0\n",
" for n in range(N):\n",
" y[n+1] = y[n] + h * f(t[n], y[n])\n",
" y[n + 1] = y[n] + h * f(t[n], y[n])\n",
" return t, y\n",
"\n",
"\n",
"def heun(t0, T, y0, h, f):\n",
" N = int(T/h)\n",
" N = int(T / h)\n",
" n = len(y0)\n",
" t = np.linspace(t0, t0 + T, N+1)\n",
" y = np.zeros((N+1, n))\n",
" t = np.linspace(t0, t0 + T, N + 1)\n",
" y = np.zeros((N + 1, n))\n",
" y[0,] = y0\n",
" for n in range(N):\n",
" p1 = f(t[n], y[n])\n",
" p2 = f(t[n] + h, y[n] + h * p1)\n",
" y[n+1] = y[n] + h/2 * (p1 + p2)\n",
" y[n + 1] = y[n] + h / 2 * (p1 + p2)\n",
" return t, y\n",
"\n",
"\n",
"def runge(t0, T, y0, h, f):\n",
" N = int(T/h)\n",
" N = int(T / h)\n",
" n = len(y0)\n",
" t = np.linspace(t0, t0 + T, N+1)\n",
" y = np.zeros((N+1, n))\n",
" t = np.linspace(t0, t0 + T, N + 1)\n",
" y = np.zeros((N + 1, n))\n",
" y[0,] = y0\n",
" for n in range(N):\n",
" p1 = f(t[n], y[n])\n",
" p2 = f(t[n] + h/2, y[n] + h/2 * p1)\n",
" p3 = f(t[n] + h/2, y[n] + h/2 * p2)\n",
" p4 = f(t[n] + h, y[n] + h* p3)\n",
" y[n+1] = y[n] + h/6 * (p1 + 2*p2 + 2*p3 + p4)\n",
" p2 = f(t[n] + h / 2, y[n] + h / 2 * p1)\n",
" p3 = f(t[n] + h / 2, y[n] + h / 2 * p2)\n",
" p4 = f(t[n] + h, y[n] + h * p3)\n",
" y[n + 1] = y[n] + h / 6 * (p1 + 2 * p2 + 2 * p3 + p4)\n",
" return t, y"
]
},
@@ -624,13 +637,16 @@
"# Question 1\n",
"a, b, c = 0.1, 2, 1\n",
"\n",
"\n",
"# f second membre de l'EDO\n",
"def f1(t,y):\n",
" return c * y * (1 - y/b)\n",
"def f1(t, y):\n",
" return c * y * (1 - y / b)\n",
"\n",
"\n",
"# sol. exacte de (P1)\n",
"def yex1(t):\n",
" return b / (1 + (b-a)/a * np.exp(-c*t))\n",
" return b / (1 + (b - a) / a * np.exp(-c * t))\n",
"\n",
"\n",
"t0, T = 0, 15\n",
"h = 0.2\n",
@@ -639,25 +655,27 @@
"plt.figure()\n",
"for schema in [euler_explicite, heun, runge]:\n",
" t, y_app = schema(t0, T, y0, h, f1)\n",
" plt.plot(t, y_app.ravel(), '+ ', label=f'Solution approchée par {schema.__name__}')\n",
" plt.plot(t, y_app.ravel(), \"+ \", label=f\"Solution approchée par {schema.__name__}\")\n",
"\n",
"t = np.arange(t0, t0+T+h, h)\n",
"t = np.arange(t0, t0 + T + h, h)\n",
"y_ex = yex1(t)\n",
"plt.plot(t, y_ex, label='Solution exacte', lw=1)\n",
"plt.title(f'Solutions approchées pour le schema {schema.__name__}')\n",
"plt.xlabel(f'$t$')\n",
"plt.ylabel(f'$y$')\n",
"plt.plot(t, y_ex, label=\"Solution exacte\", lw=1)\n",
"plt.title(f\"Solutions approchées pour le schema {schema.__name__}\")\n",
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(\"$y$\")\n",
"plt.legend(fontsize=7)\n",
"\n",
"plt.figure()\n",
"for schema in [euler_explicite, heun, runge]: \n",
"for schema in [euler_explicite, heun, runge]:\n",
" t, y_app = schema(t0, T, y0, h, f1)\n",
" plt.plot(t, np.abs(y_app.ravel() - yex1(t)), label=f'Schema {schema.__name__}')\n",
" \n",
" plt.plot(t, np.abs(y_app.ravel() - yex1(t)), label=f\"Schema {schema.__name__}\")\n",
"\n",
"plt.legend()\n",
"plt.xlabel(f'$t$')\n",
"plt.ylabel('$|y(t_n) - y^n|$')\n",
"plt.title(f'différence en valeur absolue entre sol. exacte et sol. approchée, pour différents schemas')"
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(\"$|y(t_n) - y^n|$\")\n",
"plt.title(\n",
" \"différence en valeur absolue entre sol. exacte et sol. approchée, pour différents schemas\"\n",
")"
]
},
{
@@ -696,24 +714,26 @@
" Z[1] = -Y[0] + np.cos(t)\n",
" return Z\n",
"\n",
"\n",
"def yexF(t):\n",
" return 1/2 * np.sin(t) * t + 5 * np.cos(t) + np.sin(t),\n",
" return (1 / 2 * np.sin(t) * t + 5 * np.cos(t) + np.sin(t),)\n",
"\n",
"\n",
"t0, T = 0, 15\n",
"h = 0.2\n",
"Y0 = np.array([5,1])\n",
"Y0 = np.array([5, 1])\n",
"\n",
"plt.figure()\n",
"for schema in [euler_explicite, heun, runge]:\n",
" t, y_app = schema(t0, T, Y0, h, F)\n",
" plt.plot(t, y_app[:,0], '+ ', label=f'Solution approchée par {schema.__name__}')\n",
" plt.plot(t, y_app[:, 0], \"+ \", label=f\"Solution approchée par {schema.__name__}\")\n",
"\n",
"t = np.arange(t0, t0+T+h, h)\n",
"t = np.arange(t0, t0 + T + h, h)\n",
"y_ex = yexF(t)\n",
"plt.plot(t, y_ex[0], label='Solution exacte', lw=3)\n",
"plt.title('Solutions approchées par differents schemas')\n",
"plt.xlabel(f'$t$')\n",
"plt.ylabel(f'$y$')\n",
"plt.plot(t, y_ex[0], label=\"Solution exacte\", lw=3)\n",
"plt.title(\"Solutions approchées par differents schemas\")\n",
"plt.xlabel(\"$t$\")\n",
"plt.ylabel(\"$y$\")\n",
"plt.legend(fontsize=7)"
]
},
@@ -730,16 +750,19 @@
"import matplotlib.pyplot as plt\n",
"# Question 3\n",
"\n",
"\n",
"# f second membre de l'EDO\n",
"def f3(t,y):\n",
" return (np.cos(t) - y) / (1+t)\n",
"def f3(t, y):\n",
" return (np.cos(t) - y) / (1 + t)\n",
"\n",
"\n",
"# sol. exacte de (P1)\n",
"def yex3(t):\n",
" return (np.sin(t) - 1/4) / (1 + t)\n",
" return (np.sin(t) - 1 / 4) / (1 + t)\n",
"\n",
"\n",
"t0, T = 0, 10\n",
"y0 = np.array([-1/4])"
"y0 = np.array([-1 / 4])"
]
},
{
@@ -793,15 +816,17 @@
"plt.figure()\n",
"for schema in [euler_explicite, heun, runge]:\n",
" plt.figure()\n",
" for s in range(2,11):\n",
" h=1/(2**s)\n",
" for s in range(2, 11):\n",
" h = 1 / (2**s)\n",
" t, y_app = schema(t0, T, y0, h, f3)\n",
" plt.plot(t, y_app, label=f'h={h}')\n",
" plt.plot(t, y_app, label=f\"h={h}\")\n",
"\n",
" plt.legend(fontsize=7)\n",
" plt.xlabel(f'$t$')\n",
" plt.ylabel('$y^n$')\n",
" plt.title(f'Solutions approchées de (P) obtenus avec {schema.__name__} pour différentes valeurs du pas h')"
" plt.xlabel(\"$t$\")\n",
" plt.ylabel(\"$y^n$\")\n",
" plt.title(\n",
" f\"Solutions approchées de (P) obtenus avec {schema.__name__} pour différentes valeurs du pas h\"\n",
" )"
]
},
{
@@ -846,22 +871,24 @@
"for schema in [euler_explicite, heun, runge]:\n",
" H, E = [], []\n",
" plt.figure()\n",
" for s in range(2,11):\n",
" h=1/(2**s)\n",
" for s in range(2, 11):\n",
" h = 1 / (2**s)\n",
" t, y_app = schema(t0, T, y0, h, f3)\n",
" E.append(np.max(np.abs(yex3(t) - y_app.ravel())))\n",
" H.append(h)\n",
" \n",
"\n",
" H, E = np.array(H), np.array(E)\n",
" plt.plot(np.log(H), np.log(E), '+-', label='erreur')\n",
" plt.plot(np.log(H), np.log(H), '-', label='droite pente 1')\n",
" plt.plot(np.log(H), 2*np.log(H), '-', label='droite pente 2')\n",
" plt.plot(np.log(H), 3*np.log(H), '-', label='droite pente 3')\n",
" \n",
" plt.plot(np.log(H), np.log(E), \"+-\", label=\"erreur\")\n",
" plt.plot(np.log(H), np.log(H), \"-\", label=\"droite pente 1\")\n",
" plt.plot(np.log(H), 2 * np.log(H), \"-\", label=\"droite pente 2\")\n",
" plt.plot(np.log(H), 3 * np.log(H), \"-\", label=\"droite pente 3\")\n",
"\n",
" plt.legend()\n",
" plt.title(f'Erreur pour la méthode {schema.__name__} en echelle logarithmique : log(E) en fonction de log(h)')\n",
" plt.xlabel(f'$log(h)$')\n",
" plt.ylabel(f'$log(E)$')"
" plt.title(\n",
" f\"Erreur pour la méthode {schema.__name__} en echelle logarithmique : log(E) en fonction de log(h)\"\n",
" )\n",
" plt.xlabel(\"$log(h)$\")\n",
" plt.ylabel(\"$log(E)$\")"
]
},
{
@@ -943,19 +970,20 @@
"a = 5\n",
"t0, T = 0, 20\n",
"H = [0.1, 0.05, 0.025]\n",
"Y0 = np.array([1/10, 1])\n",
"Y0 = np.array([1 / 10, 1])\n",
"\n",
"\n",
"def P(t, Y):\n",
" return np.array([Y[1], (a - Y[0]**2) * Y[1] - Y[0]])\n",
" return np.array([Y[1], (a - Y[0] ** 2) * Y[1] - Y[0]])\n",
"\n",
"\n",
"for schema in [euler_explicite, heun, runge]:\n",
" plt.figure()\n",
" for h in H:\n",
" t, y_app = schema(t0, T, Y0, h, P)\n",
" plt.plot(t, y_app[:,0], '+ ', label=f'Solution approchée pour h={h}')\n",
" plt.plot(t, y_app[:, 0], \"+ \", label=f\"Solution approchée pour h={h}\")\n",
" plt.legend()\n",
" plt.title(f'Solutions approchees pour differents pas h par {schema.__name__}')\n",
" "
" plt.title(f\"Solutions approchees pour differents pas h par {schema.__name__}\")"
]
},
{

207
L3/Méthodes Numériques/TP1_Equation_de_Poisson.ipynb
View File

@@ -240,22 +240,25 @@
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"def A(M):\n",
" return 2 * np.eye(M) - np.eye(M, k=-1) - np.eye(M, k=1)\n",
"\n",
"\n",
"def solution_approchée(f, M, a=0, b=1, ua=0, ub=0):\n",
" X = np.linspace(a, b + 1 if b == 0 else b, M+2)\n",
" U = np.zeros(M+2)\n",
" X = np.linspace(a, b + 1 if b == 0 else b, M + 2)\n",
" U = np.zeros(M + 2)\n",
" U[0], U[-1] = ua, ub\n",
" \n",
" h = (b-a)/(M+1)\n",
"\n",
" h = (b - a) / (M + 1)\n",
"\n",
" delta = np.zeros(M)\n",
" delta[0], delta[-1] = ua/np.power(h, 2), ub/np.power(h, 2)\n",
" \n",
" delta[0], delta[-1] = ua / np.power(h, 2), ub / np.power(h, 2)\n",
"\n",
" U[1:-1] = np.linalg.solve(A(M), h**2 * (f(X)[1:-1] + delta))\n",
" return X, U"
]
@@ -290,22 +293,25 @@
],
"source": [
"M = 50\n",
"h = 1/(M+1)\n",
"h = 1 / (M + 1)\n",
"\n",
"\n",
"def f1(x):\n",
" return (2 * np.pi)**2 * np.sin(2 * np.pi * x)\n",
" return (2 * np.pi) ** 2 * np.sin(2 * np.pi * x)\n",
"\n",
"\n",
"def u1(x):\n",
" return np.sin(2 * np.pi * x)\n",
"\n",
"\n",
"x, U_app = solution_approchée(f1, M)\n",
"plt.figure()\n",
"plt.scatter(x, U_app, label=f\"$A_M U_h = h² F$\")\n",
"plt.plot(x, u1(x), label='Solution exacte', color='red')\n",
"plt.scatter(x, U_app, label=\"$A_M U_h = h² F$\")\n",
"plt.plot(x, u1(x), label=\"Solution exacte\", color=\"red\")\n",
"plt.legend()\n",
"plt.ylabel('f(x)')\n",
"plt.xlabel('x')\n",
"plt.title(f\"$-u''=f(x)$\")"
"plt.ylabel(\"f(x)\")\n",
"plt.xlabel(\"x\")\n",
"plt.title(\"$-u''=f(x)$\")"
]
},
{
@@ -340,21 +346,23 @@
"H, E = [], []\n",
"for k in range(2, 12):\n",
" M = 2**k\n",
" h = 1/(M+1)\n",
" \n",
" h = 1 / (M + 1)\n",
"\n",
" x, U_app = solution_approchée(f1, M)\n",
" \n",
"\n",
" e = np.abs(u1(x) - U_app)\n",
" plt.plot(x, e, label=f'{h}')\n",
" plt.plot(x, e, label=f\"{h}\")\n",
" E.append(np.max(e))\n",
" H.append(h)\n",
"\n",
"H, E = np.array(H), np.array(E) \n",
"H, E = np.array(H), np.array(E)\n",
"\n",
"plt.xlabel('$h$')\n",
"plt.ylabel('$\\max_{j=0,\\dots,M+1}|u(x_j)-u_j|$')\n",
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"$\\max_{j=0,\\dots,M+1}|u(x_j)-u_j|$\")\n",
"plt.legend(fontsize=7)\n",
"plt.title('Différence en valeur absolue entre la solution exacte et la solution approchée')"
"plt.title(\n",
" \"Différence en valeur absolue entre la solution exacte et la solution approchée\"\n",
")"
]
},
{
@@ -387,12 +395,12 @@
],
"source": [
"plt.figure()\n",
"plt.plot(H, E/H, label=f'$E_h / h$')\n",
"plt.plot(H, E/H**2, label=f'$E_h / h^2$')\n",
"plt.plot(H, E / H, label=\"$E_h / h$\")\n",
"plt.plot(H, E / H**2, label=\"$E_h / h^2$\")\n",
"plt.legend()\n",
"plt.xlabel('$h$')\n",
"plt.ylabel('Erreur globale')\n",
"plt.title('Erreur globale en fonction de $h$')"
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"Erreur globale\")\n",
"plt.title(\"Erreur globale en fonction de $h$\")"
]
},
{
@@ -425,13 +433,15 @@
],
"source": [
"plt.figure()\n",
"plt.plot(np.log(H), np.log(E), '+-', label='erreur')\n",
"plt.plot(np.log(H), np.log(H), '-', label='droite pente 1')\n",
"plt.plot(np.log(H), 2*np.log(H), '-', label='droite pente 2')\n",
"plt.plot(np.log(H), np.log(E), \"+-\", label=\"erreur\")\n",
"plt.plot(np.log(H), np.log(H), \"-\", label=\"droite pente 1\")\n",
"plt.plot(np.log(H), 2 * np.log(H), \"-\", label=\"droite pente 2\")\n",
"plt.legend()\n",
"plt.xlabel(f'$log(h)$')\n",
"plt.ylabel(f'$log(E)$')\n",
"plt.title('Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)')"
"plt.xlabel(\"$log(h)$\")\n",
"plt.ylabel(\"$log(E)$\")\n",
"plt.title(\n",
" \"Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)\"\n",
")"
]
},
{
@@ -509,27 +519,31 @@
],
"source": [
"import numpy as np\n",
"\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"ua, ub = 1/2, 1/3\n",
"ua, ub = 1 / 2, 1 / 3\n",
"a, b = 1, 2\n",
"M = 49\n",
"\n",
"\n",
"def f2(x):\n",
" return -2/np.power((x+1), 3)\n",
" return -2 / np.power((x + 1), 3)\n",
"\n",
"\n",
"def u2(x):\n",
" return 1/(x+1)\n",
" return 1 / (x + 1)\n",
"\n",
"\n",
"x, U_app = solution_approchée(f2, M, a, b, ua, ub)\n",
"plt.figure()\n",
"plt.scatter(x, U_app, label=f\"$A_M U_h = h² F$\")\n",
"plt.plot(x, u2(x), label='Solution exacte', color='red')\n",
"plt.scatter(x, U_app, label=\"$A_M U_h = h² F$\")\n",
"plt.plot(x, u2(x), label=\"Solution exacte\", color=\"red\")\n",
"plt.legend()\n",
"plt.ylabel('f(x)')\n",
"plt.xlabel('x')\n",
"plt.title(f\"$-u''=f(x)$\")"
"plt.ylabel(\"f(x)\")\n",
"plt.xlabel(\"x\")\n",
"plt.title(\"$-u''=f(x)$\")"
]
},
{
@@ -565,19 +579,21 @@
"H, E = [], []\n",
"for k in range(2, 12):\n",
" M = 2**k\n",
" h = (b-a)/(M+1)\n",
" \n",
" h = (b - a) / (M + 1)\n",
"\n",
" x, U_app = solution_approchée(f2, M, a, b, ua, ub)\n",
" \n",
"\n",
" e = np.abs(u2(x) - U_app)\n",
" plt.plot(x, e, label=f'{h}')\n",
" plt.plot(x, e, label=f\"{h}\")\n",
" E.append(np.max(e))\n",
" H.append(h)\n",
" \n",
"plt.xlabel('$h$')\n",
"plt.ylabel('$\\max_{j=0,\\dots,M+1}|u(x_j)-u_j|$')\n",
"\n",
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"$\\max_{j=0,\\dots,M+1}|u(x_j)-u_j|$\")\n",
"plt.legend(fontsize=7)\n",
"plt.title('Différence en valeur absolue entre la solution exacte et la solution approchée')"
"plt.title(\n",
" \"Différence en valeur absolue entre la solution exacte et la solution approchée\"\n",
")"
]
},
{
@@ -609,14 +625,14 @@
}
],
"source": [
"H, E = np.array(H), np.array(E) \n",
"H, E = np.array(H), np.array(E)\n",
"plt.figure()\n",
"plt.plot(H, E/H, label=f'$E_h / h$')\n",
"plt.plot(H, E/H**2, label=f'$E_h / h^2$')\n",
"plt.plot(H, E / H, label=\"$E_h / h$\")\n",
"plt.plot(H, E / H**2, label=\"$E_h / h^2$\")\n",
"plt.legend()\n",
"plt.xlabel('$h$')\n",
"plt.ylabel('Erreur globale')\n",
"plt.title('Erreur globale en fonction de $h$')"
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"Erreur globale\")\n",
"plt.title(\"Erreur globale en fonction de $h$\")"
]
},
{
@@ -649,13 +665,15 @@
],
"source": [
"plt.figure()\n",
"plt.plot(np.log(H), np.log(E), '+-', label='erreur')\n",
"plt.plot(np.log(H), np.log(H), '-', label='droite pente 1')\n",
"plt.plot(np.log(H), 2*np.log(H), '-', label='droite pente 2')\n",
"plt.plot(np.log(H), np.log(E), \"+-\", label=\"erreur\")\n",
"plt.plot(np.log(H), np.log(H), \"-\", label=\"droite pente 1\")\n",
"plt.plot(np.log(H), 2 * np.log(H), \"-\", label=\"droite pente 2\")\n",
"plt.legend()\n",
"plt.xlabel(f'$log(h)$')\n",
"plt.ylabel(f'$log(E)$')\n",
"plt.title('Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)')"
"plt.xlabel(\"$log(h)$\")\n",
"plt.ylabel(\"$log(E)$\")\n",
"plt.title(\n",
" \"Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)\"\n",
")"
]
},
{
@@ -732,22 +750,25 @@
],
"source": [
"def An(M, h):\n",
" mat = 2*np.eye(M) - np.eye(M, k=1) - np.eye(M, k=-1)\n",
" mat[0,0] = 1\n",
" mat = 2 * np.eye(M) - np.eye(M, k=1) - np.eye(M, k=-1)\n",
" mat[0, 0] = 1\n",
" mat[-1, -1] = 1\n",
" return mat + h**2 * np.eye(M)\n",
"\n",
"\n",
"def f3(x):\n",
" return (np.power(2 * np.pi, 2) + 1 ) * np.cos(2*np.pi * x)\n",
" return (np.power(2 * np.pi, 2) + 1) * np.cos(2 * np.pi * x)\n",
"\n",
"\n",
"def u3(x):\n",
" return np.cos(2*np.pi * x)\n",
" return np.cos(2 * np.pi * x)\n",
"\n",
"\n",
"def solution_neumann(f, M, a, b):\n",
" X = np.linspace(a, b, M+2)\n",
" U = np.zeros(M+2)\n",
" h = 1/(M+1)\n",
" \n",
" X = np.linspace(a, b, M + 2)\n",
" U = np.zeros(M + 2)\n",
" h = 1 / (M + 1)\n",
"\n",
" U[1:-1] = np.linalg.solve(An(M, h), h**2 * f3(X[1:-1]))\n",
" U[0], U[-1] = U[1], U[-2]\n",
" return X, U\n",
@@ -759,12 +780,14 @@
"for M in [49, 99, 499]:\n",
" x, U_app = solution_neumann(f3, M, a, b)\n",
"\n",
" plt.scatter(x, U_app, marker='+', s=3, label=\"$h^2({A_N}_h + I_M)U = h^2F$ pour M={M}\")\n",
"plt.plot(x, u3(x), label='Solution exacte', color='red')\n",
" plt.scatter(\n",
" x, U_app, marker=\"+\", s=3, label=\"$h^2({A_N}_h + I_M)U = h^2F$ pour M={M}\"\n",
" )\n",
"plt.plot(x, u3(x), label=\"Solution exacte\", color=\"red\")\n",
"plt.legend(fontsize=8)\n",
"plt.ylabel('f(x)')\n",
"plt.xlabel('x')\n",
"plt.title(f\"$-u''(x) + u(x)=f(x)$\")"
"plt.ylabel(\"f(x)\")\n",
"plt.xlabel(\"x\")\n",
"plt.title(\"$-u''(x) + u(x)=f(x)$\")"
]
},
{
@@ -800,19 +823,21 @@
"H, E = [], []\n",
"for k in range(2, 12):\n",
" M = 2**k\n",
" h = (b-a)/(M+1)\n",
" \n",
" h = (b - a) / (M + 1)\n",
"\n",
" x, U_app = solution_neumann(f3, M, a, b)\n",
" \n",
"\n",
" e = np.abs(u3(x) - U_app)\n",
" plt.plot(x, e, label=f'{h}')\n",
" plt.plot(x, e, label=f\"{h}\")\n",
" E.append(np.max(e))\n",
" H.append(h)\n",
" \n",
"plt.xlabel('$h$')\n",
"plt.ylabel('$\\max_{j=0,\\dots,M+1}|u(x_j)-u_j|$')\n",
"\n",
"plt.xlabel(\"$h$\")\n",
"plt.ylabel(\"$\\max_{j=0,\\dots,M+1}|u(x_j)-u_j|$\")\n",
"plt.legend(fontsize=7)\n",
"plt.title('Différence en valeur absolue entre la solution exacte et la solution approchée')"
"plt.title(\n",
" \"Différence en valeur absolue entre la solution exacte et la solution approchée\"\n",
")"
]
},
{
@@ -845,13 +870,15 @@
],
"source": [
"plt.figure()\n",
"plt.plot(np.log(H), np.log(E), '+-', label='erreur')\n",
"plt.plot(np.log(H), np.log(H), '-', label='droite pente 1')\n",
"plt.plot(np.log(H), 2*np.log(H), '-', label='droite pente 2')\n",
"plt.plot(np.log(H), np.log(E), \"+-\", label=\"erreur\")\n",
"plt.plot(np.log(H), np.log(H), \"-\", label=\"droite pente 1\")\n",
"plt.plot(np.log(H), 2 * np.log(H), \"-\", label=\"droite pente 2\")\n",
"plt.legend()\n",
"plt.xlabel(f'$log(h)$')\n",
"plt.ylabel(f'$log(E)$')\n",
"plt.title('Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)')"
"plt.xlabel(\"$log(h)$\")\n",
"plt.ylabel(\"$log(E)$\")\n",
"plt.title(\n",
" \"Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)\"\n",
")"
]
},
{
@@ -959,7 +986,11 @@
"outputs": [],
"source": [
"def v(x):\n",
" return 4*np.exp(-500*np.square(x-.8)) + np.exp(-50*np.square(x-.2))+.5*np.random.random(x.shape)"
" return (\n",
" 4 * np.exp(-500 * np.square(x - 0.8))\n",
" + np.exp(-50 * np.square(x - 0.2))\n",
" + 0.5 * np.random.random(x.shape)\n",
" )"
]
}
],

223
L3/Projet Numérique/Segregation.ipynb
View File

File diff suppressed because one or more lines are too long

54
L3/Statistiques/TP1.ipynb
View File

@@ -48,10 +48,10 @@
"nb_repl = 100000\n",
"sample = np.random.exponential(lam, nb_repl)\n",
"intervalle = np.linspace(0, 5, 100)\n",
"plt.hist(sample, bins=intervalle, density=True, label='Echantillon')\n",
"plt.hist(sample, bins=intervalle, density=True, label=\"Echantillon\")\n",
"\n",
"densite = stats.expon().pdf\n",
"plt.plot(intervalle, densite(intervalle), label='Fonction densité')\n",
"plt.plot(intervalle, densite(intervalle), label=\"Fonction densité\")\n",
"plt.legend()"
]
},
@@ -86,15 +86,15 @@
],
"source": [
"np_repl = 10000\n",
"sampleX, sampleY = 2*np.random.rand(np_repl)-1, 2*np.random.rand(np_repl)-1\n",
"sampleX, sampleY = 2 * np.random.rand(np_repl) - 1, 2 * np.random.rand(np_repl) - 1\n",
"\n",
"intervalle = np.linspace(-2, 2, 100)\n",
"plt.hist(sampleX + sampleY, bins=intervalle, density=True, label='Echantillon')\n",
"plt.hist(sampleX + sampleY, bins=intervalle, density=True, label=\"Echantillon\")\n",
"\n",
"densite = stats.uniform(-1, 2).pdf\n",
"plt.plot(intervalle, densite(intervalle), label='Fonction densité')\n",
"plt.plot(intervalle, densite(intervalle), label=\"Fonction densité\")\n",
"\n",
"plt.plot(intervalle, 1/4 * np.maximum(2 - np.abs(intervalle), 0), label='Fonction')\n",
"plt.plot(intervalle, 1 / 4 * np.maximum(2 - np.abs(intervalle), 0), label=\"Fonction\")\n",
"plt.legend()"
]
},
@@ -130,10 +130,22 @@
"source": [
"lam = 1\n",
"nb_repl = 10000\n",
"sampleX, sampleY, sampleZ = np.random.exponential(lam, nb_repl), np.random.exponential(lam, nb_repl), np.random.exponential(lam, nb_repl)\n",
"sampleX, sampleY, sampleZ = (\n",
" np.random.exponential(lam, nb_repl),\n",
" np.random.exponential(lam, nb_repl),\n",
" np.random.exponential(lam, nb_repl),\n",
")\n",
"intervalle = np.linspace(-5, 5, 100)\n",
"plt.hist(sampleX - sampleY/2, bins=intervalle, density=True, alpha=.7, label='Echantillon de X - Y/2')\n",
"plt.hist(sampleZ/2, bins=intervalle, density=True, alpha=.7, label='Echantillon de Z')\n",
"plt.hist(\n",
" sampleX - sampleY / 2,\n",
" bins=intervalle,\n",
" density=True,\n",
" alpha=0.7,\n",
" label=\"Echantillon de X - Y/2\",\n",
")\n",
"plt.hist(\n",
" sampleZ / 2, bins=intervalle, density=True, alpha=0.7, label=\"Echantillon de Z\"\n",
")\n",
"plt.legend()"
]
},
@@ -157,15 +169,17 @@
}
],
"source": [
"p = 1/2\n",
"p = 1 / 2\n",
"nb_repl = 1000\n",
"\n",
"plt.figure(figsize=(18, 5))\n",
"for i, k in enumerate([10, 100, 1000]):\n",
" plt.subplot(1, 3, i+1)\n",
" plt.subplot(1, 3, i + 1)\n",
" sample = np.random.binomial(k, p, nb_repl)\n",
" intervalle = np.linspace(np.min(sample), np.max(sample), 100)\n",
" plt.hist(sample, bins=intervalle, density=True, label=f'Echantillon de X pour n={k}')\n",
" plt.hist(\n",
" sample, bins=intervalle, density=True, label=f\"Echantillon de X pour n={k}\"\n",
" )\n",
" plt.legend()"
]
},
@@ -179,7 +193,7 @@
"outputs": [],
"source": [
"def sample_uniforme(N):\n",
" return 2*np.random.rand(N) - 1"
" return 2 * np.random.rand(N) - 1"
]
},
{
@@ -218,7 +232,7 @@
"\n",
"for _ in range(nb_lgn):\n",
" liste_Sn.append(np.mean(sample_uniforme(nb_repl)))\n",
" \n",
"\n",
"nb_bins = 100\n",
"intervalles = np.linspace(np.min(liste_Sn), np.max(liste_Sn), nb_bins)\n",
"plt.hist(liste_Sn, density=True, bins=intervalles)\n",
@@ -256,12 +270,14 @@
"liste_Sn = []\n",
"\n",
"for _ in range(nb_lgn):\n",
" liste_Sn.append(np.mean(np.sqrt(3*nb_repl) * np.tan(np.pi/2 * sample_uniforme(nb_repl))))\n",
" liste_Sn.append(\n",
" np.mean(np.sqrt(3 * nb_repl) * np.tan(np.pi / 2 * sample_uniforme(nb_repl)))\n",
" )\n",
"\n",
"#nb_bins = 100\n",
"#intervalles = np.linspace(np.min(liste_Sn), np.max(liste_Sn), nb_bins)\n",
"#plt.hist(liste_Sn, density=True, bins=intervalles)\n",
"#plt.show()\n",
"# nb_bins = 100\n",
"# intervalles = np.linspace(np.min(liste_Sn), np.max(liste_Sn), nb_bins)\n",
"# plt.hist(liste_Sn, density=True, bins=intervalles)\n",
"# plt.show()\n",
"\n",
"plt.figure()\n",
"densite = stats.norm(scale=1).pdf\n",

47
L3/Statistiques/TP2.ipynb
View File

@@ -27,7 +27,12 @@
"source": [
"def f(L, z):\n",
" x, y = L[0], L[1]\n",
" return np.power(z*x, 2) + np.power(y/z, 2) - np.cos(2 * np.pi * x) - np.cos(2 * np.pi * y)"
" return (\n",
" np.power(z * x, 2)\n",
" + np.power(y / z, 2)\n",
" - np.cos(2 * np.pi * x)\n",
" - np.cos(2 * np.pi * y)\n",
" )"
]
},
{
@@ -123,17 +128,17 @@
"source": [
"plt.figure(figsize=(18, 5))\n",
"for i, nb_repl in enumerate([100, 1000, 10000]):\n",
" plt.subplot(1, 3, i+1)\n",
" plt.subplot(1, 3, i + 1)\n",
" sample_X1 = np.random.normal(0, 1, nb_repl)\n",
" sample_X2 = np.random.normal(3, np.sqrt(5), nb_repl)\n",
" sample_e = np.random.normal(0, np.sqrt(1/4), nb_repl)\n",
" sample_e = np.random.normal(0, np.sqrt(1 / 4), nb_repl)\n",
" Y = 5 * sample_X1 - 4 * sample_X2 + 2 + sample_e\n",
"\n",
" intervalle = np.linspace(np.min(Y), np.max(Y), 100)\n",
" plt.hist(Y, bins=intervalle, density=True, label='Echantillon de Y')\n",
" plt.hist(Y, bins=intervalle, density=True, label=\"Echantillon de Y\")\n",
"\n",
" densite = stats.norm(-10, np.sqrt(105.25)).pdf\n",
" plt.plot(intervalle, densite(intervalle), label='Fonction densité')\n",
" plt.plot(intervalle, densite(intervalle), label=\"Fonction densité\")\n",
"\n",
" plt.title(f\"Graphique de la somme de gaussiennes pour N={nb_repl}\")\n",
" plt.legend()"
@@ -151,8 +156,9 @@
"def theta_hat(Y):\n",
" return np.mean(Y)\n",
"\n",
"\n",
"def sigma_hat(Y):\n",
" return 1/nb_repl * np.sum(np.power(Y - theta_hat(Y), 2))"
" return 1 / nb_repl * np.sum(np.power(Y - theta_hat(Y), 2))"
]
},
{
@@ -166,7 +172,11 @@
"source": [
"def log_likehood_gauss(X, Y):\n",
" theta, sigma_2 = X[0], X[1]\n",
" return 1/2*np.log(2*np.pi) + 1/2*np.log(sigma_2) + 1/(2*nb_repl*sigma_2) * np.sum(np.power(Y - theta, 2))"
" return (\n",
" 1 / 2 * np.log(2 * np.pi)\n",
" + 1 / 2 * np.log(sigma_2)\n",
" + 1 / (2 * nb_repl * sigma_2) * np.sum(np.power(Y - theta, 2))\n",
" )"
]
},
{
@@ -191,9 +201,9 @@
"nb_repl = 5000\n",
"sample_X1 = np.random.normal(0, 1, nb_repl)\n",
"sample_X2 = np.random.normal(3, np.sqrt(5), nb_repl)\n",
"sample_e = np.random.normal(0, np.sqrt(1/4), nb_repl)\n",
"sample_e = np.random.normal(0, np.sqrt(1 / 4), nb_repl)\n",
"Y = 5 * sample_X1 - 4 * sample_X2 + 2 + sample_e\n",
" \n",
"\n",
"mk = {\"method\": \"BFGS\", \"args\": Y}\n",
"res = opt.basinhopping(log_likehood_gauss, x0=(-1, 98.75), minimizer_kwargs=mk)\n",
"print(res.x)\n",
@@ -211,12 +221,12 @@
"outputs": [],
"source": [
"def simule(a, b, n):\n",
" X = np.random.gamma(a, 1/b, n)\n",
" X = np.random.gamma(a, 1 / b, n)\n",
" intervalle = np.linspace(0, np.max(X), 100)\n",
" plt.hist(X, bins=intervalle, density=True, label='Echantillon de X')\n",
" plt.hist(X, bins=intervalle, density=True, label=\"Echantillon de X\")\n",
"\n",
" densite = stats.gamma.pdf(intervalle, a, 0, 1/b)\n",
" plt.plot(intervalle, densite, label='Fonction densité Gamma(2, 1)')\n",
" densite = stats.gamma.pdf(intervalle, a, 0, 1 / b)\n",
" plt.plot(intervalle, densite, label=\"Fonction densité Gamma(2, 1)\")\n",
" plt.legend()"
]
},
@@ -254,8 +264,13 @@
"source": [
"def log_likehood_gamma(X, sample):\n",
" a, b = X[0], X[1]\n",
" n = len(sample) \n",
" return -n*a*np.log(b) + n * np.log(sp.gamma(a)) - (a-1) * np.sum(np.log(sample)) + b * np.sum(sample)"
" n = len(sample)\n",
" return (\n",
" -n * a * np.log(b)\n",
" + n * np.log(sp.gamma(a))\n",
" - (a - 1) * np.sum(np.log(sample))\n",
" + b * np.sum(sample)\n",
" )"
]
},
{
@@ -296,7 +311,7 @@
"nb_repl = 1000\n",
"a, b = 2, 1\n",
"\n",
"sample = np.random.gamma(a, 1/b, nb_repl)\n",
"sample = np.random.gamma(a, 1 / b, nb_repl)\n",
"mk = {\"method\": \"BFGS\", \"args\": sample}\n",
"res = opt.basinhopping(log_likehood_gamma, x0=(1, 1), minimizer_kwargs=mk)\n",
"print(res.x)"