2022-03-15 07:04:08 +00:00
{
"cells": [
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 1,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
2022-04-03 20:18:44 +00:00
"import matplotlib.pyplot as plt"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 2,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"def maillage_carre(n: int):\n",
" \"\"\"\n",
" Une discrétisation possible d'une EDP elliptique sur le domaine ]0,1[ x ]0,1[.\n",
" Le carre [0,1]x[0,1] est maille uniquement avec des triangles.\n",
" Les conditions limites sont de type Dirichlet uniquement -> `neumann=[]`.\n",
"\n",
" Args:\n",
2022-04-03 20:18:44 +00:00
" n: nombre de points par cote du care => Npts points de discretisation au total\n",
2022-03-15 07:04:08 +00:00
"\n",
" Returns:\n",
2022-04-03 20:18:44 +00:00
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" elements3: matrice a trois colonnes. Chaque ligne contient les indices des sommets d'un element triangle, dans le sens antihoraire.\n",
" dirichlet: vecteur colonne des indices des sommets de la frontiere de Dirichlet.\n",
" neumann: matrice a deux colonnes. Chaque ligne contient les indices des deux sommets d'une arete de la frontiere de Neumann. (neumann est vide sur cet exemple)\n",
2022-03-15 07:04:08 +00:00
" \"\"\"\n",
"\n",
" h = 1 / (n - 1)\n",
" n_pts = n * n\n",
" n_elm = 2 * (n - 1) * (n - 1)\n",
" coordinates = np.zeros((n_pts, 2))\n",
" elements3 = np.zeros((n_elm, 3), dtype=int)\n",
" neumann = []\n",
" dirichlet = np.zeros((4 * n - 4, 1), dtype=int)\n",
"\n",
" # Coordonnees et connectivites :\n",
" e = -1\n",
" p = -1\n",
" x = np.zeros((n + 1, 1))\n",
" x[n, 0] = 1.0\n",
"\n",
" for l in range(n + 1):\n",
" x[l, 0] = l * h\n",
"\n",
" for j in range(n):\n",
" for i in range(n):\n",
" p = p + 1\n",
" coordinates[p, 0] = x[i, 0]\n",
" coordinates[p, 1] = x[j, 0]\n",
" if (i != n - 1) & (j != n - 1):\n",
" p1 = p\n",
" p2 = p1 + 1\n",
" p3 = p1 + n\n",
" p4 = p2 + n\n",
" e = e + 1\n",
" elements3[e, 0] = p1\n",
" elements3[e, 1] = p2\n",
" elements3[e, 2] = p3\n",
" e = e + 1\n",
" elements3[e, 0] = p4\n",
" elements3[e, 1] = p3\n",
" elements3[e, 2] = p2\n",
"\n",
" # Liste des sommets de la frontiere de Dirichlet:\n",
" p = -1\n",
" for j in range(n):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
" for j in range(n * 2 - 1, n * (n - 1), n):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
" for j in range(n * n - 1, n * n - n - 1, -1):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
" for j in range(n * n - 2 * n, n - 1, -n):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
2023-06-21 18:29:07 +00:00
" return coordinates, elements3, dirichlet, neumann"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 14,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def show(coordinates, u, title) -> None:\n",
2022-03-15 07:04:08 +00:00
" \"\"\"Fonction d'affichage de la solution u sur le maillage defini par elements3, coordinates.\n",
"\n",
" Args:\n",
2022-04-03 20:18:44 +00:00
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" u: vecteur colonne de longueur egale au nombre de lignes de coordinates contenant les valeurs de la solution a afficher aux points de la discretisation.\n",
2022-04-03 20:42:26 +00:00
" title: le titre de la figure\n",
2022-03-15 07:04:08 +00:00
"\n",
" Returns:\n",
" None, plots a figure\n",
" \"\"\"\n",
"\n",
2023-06-21 18:29:07 +00:00
" ax = plt.figure(figsize=(20, 20)).add_subplot(projection=\"3d\")\n",
2022-03-15 07:04:08 +00:00
" ax.plot_trisurf(\n",
" coordinates[:, 0], coordinates[:, 1], u, linewidth=0.2, antialiased=True\n",
" )\n",
2022-04-03 20:42:26 +00:00
" plt.title(title)\n",
2022-03-15 07:04:08 +00:00
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2022-04-03 20:18:44 +00:00
"## Partie I : maillage triangulaire et conditions de Dirichlet\n",
"\n",
"$$\n",
"\\left\\{\n",
"\\begin{array}{rll}\n",
"\n",
"\\displaystyle -\\delta u (x, y) &= f(x, y) &\\text{sur } \\Omega \\\\\n",
"\\displaystyle u (x, y) &= u_d(x, y) &\\text{sur } \\partial\\Omega_d \\\\\n",
"\\displaystyle \\frac{\\partial u (x, y)}{\\partial n} &= g(x, y) &\\text{sur } \\partial\\Omega_n\n",
"\n",
"\\end{array}\n",
"\\right.\n",
"$$\n"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 4,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def f(x, y) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\n",
2022-03-15 07:04:08 +00:00
"\n",
"\n",
2022-04-03 20:42:26 +00:00
"def u_ex(x, y) -> np.ndarray:\n",
2022-03-15 07:04:08 +00:00
" return np.sin(np.pi * x) * np.sin(np.pi * y)\n",
"\n",
"\n",
2022-04-03 20:42:26 +00:00
"def u_d(x, y) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" return np.zeros(x.shape[0])"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 5,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [
{
2022-04-03 20:18:44 +00:00
"name": "stdout",
"output_type": "stream",
"text": [
2023-06-21 18:29:07 +00:00
"coords =\n",
"[[0. 0. ]\n",
2022-04-03 20:18:44 +00:00
" [0.5 0. ]\n",
" [1. 0. ]\n",
" [0. 0.5]\n",
" [0.5 0.5]\n",
" [1. 0.5]\n",
" [0. 1. ]\n",
" [0.5 1. ]\n",
" [1. 1. ]]\n",
"\n",
2023-06-21 18:29:07 +00:00
"elems3 =\n",
"[[0 1 3]\n",
2022-04-03 20:18:44 +00:00
" [4 3 1]\n",
" [1 2 4]\n",
" [5 4 2]\n",
" [3 4 6]\n",
" [7 6 4]\n",
" [4 5 7]\n",
" [8 7 5]]\n",
"\n",
2023-06-21 18:29:07 +00:00
"dirichlet =\n",
"[[0]\n",
2022-04-03 20:18:44 +00:00
" [1]\n",
" [2]\n",
" [5]\n",
" [8]\n",
" [7]\n",
" [6]\n",
" [3]]\n",
"\n",
2023-06-21 18:29:07 +00:00
"neumman =\n",
"[]\n"
2022-04-03 20:18:44 +00:00
]
2022-03-15 07:04:08 +00:00
}
],
"source": [
2022-04-03 20:18:44 +00:00
"# affichage d'un petit maillage\n",
"coords, elems3, dirichlet, neumann = maillage_carre(3)\n",
"print(\n",
2023-06-21 18:29:07 +00:00
" f\"coords =\\n{coords}\",\n",
" f\"elems3 =\\n{elems3}\",\n",
" f\"dirichlet =\\n{dirichlet}\",\n",
" f\"neumman =\\n{neumann}\",\n",
2022-04-03 20:18:44 +00:00
" sep=\"\\n\\n\"\n",
")"
]
},
2022-03-15 07:04:08 +00:00
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 6,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def calcul_alpha(x, y) -> float:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Calcul du coefficient alpha.\n",
"\n",
" Args:\n",
" x (np.array): les coordonnées x du triangle.\n",
" y (np.array): les coordonnées y du triangle.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Returns:\n",
2022-04-03 20:42:26 +00:00
" alpha: le coefficient alpha.\n",
2022-04-03 20:18:44 +00:00
" \"\"\"\n",
2022-03-15 07:04:08 +00:00
" mat_alpha = np.array(\n",
" [\n",
" [x[1] - x[0], x[2] - x[0]],\n",
" [y[1] - y[0], y[2] - y[0]]\n",
" ]\n",
" )\n",
"\n",
2022-04-03 20:18:44 +00:00
" return np.linalg.det(mat_alpha)"
]
},
2022-04-03 20:42:26 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
2023-06-21 18:29:07 +00:00
"Pour calculer les coefficients de la matrice de raideur, on calcul simplement en chaque combinaison de sommets:\n",
"\n",
"$$\n",
"[M^A_T]_{ij} = \\displaystyle \\int_T \\nabla \\eta_i (x, y)^\\top \\nabla \\eta_j (x, y) \\ dx \\ dy\n",
"$$"
2022-04-03 20:42:26 +00:00
]
},
2022-04-03 20:18:44 +00:00
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 7,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1. , -0.5, -0.5],\n",
" [-0.5, 0.5, 0. ],\n",
" [-0.5, 0. , 0.5]])"
]
},
2023-06-21 18:29:07 +00:00
"execution_count": 7,
2022-04-03 20:18:44 +00:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2022-04-03 20:42:26 +00:00
"def raideur(triangle) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Construction de la matrice de raideur ́elementaire relative à un ́élément triangle.\n",
"\n",
" Args:\n",
" triangle: les coordonnées x et y des trois points formant le triangle.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Returns:\n",
" M: La matrice de raideur ́elementaire.\n",
" \"\"\"\n",
" M = np.zeros((3, 3))\n",
" x = triangle[:, 0]\n",
" y = triangle[:, 1]\n",
"\n",
" alpha = calcul_alpha(x, y)\n",
"\n",
" # calcul de la matrice M\n",
2022-03-15 07:04:08 +00:00
" for i in range(3):\n",
" grad_eta_i = np.array(\n",
" [\n",
" y[(i+1)%3] - y[(i+2)%3],\n",
" x[(i+2)%3] - x[(i+1)%3]\n",
" ]\n",
" )\n",
" for j in range(3):\n",
" grad_eta_j = np.array(\n",
" [\n",
" y[(j+1)%3] - y[(j+2)%3],\n",
" x[(j+2)%3] - x[(j+1)%3]\n",
" ]\n",
" )\n",
"\n",
" M[i, j] = np.dot(grad_eta_i, grad_eta_j)\n",
"\n",
2022-04-03 20:18:44 +00:00
" return M / alpha / 2\n",
"\n",
"# on affiche la première matrice de raideur pour vérifier\n",
"raideur(coords[elems3[0]])"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 8,
2022-03-15 07:04:08 +00:00
"metadata": {},
2022-04-03 20:18:44 +00:00
"outputs": [],
2022-03-15 07:04:08 +00:00
"source": [
2022-04-03 20:42:26 +00:00
"def assemblage(coordinates, elements3) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Assemblage de la matrice A dans le cas d'un maillage constitué uniquement d'éléments triangles.\n",
"\n",
" Args:\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" elements3: matrice a trois colonnes. Chaque ligne contient les indices des sommets d'un element triangle, dans le sens antihoraire.\n",
"\n",
" Returns:\n",
" A: matrice nécéssaire à la résolution de la formulation variationnelle du problème.\n",
" \"\"\"\n",
" Ns = len(coordinates)\n",
2022-03-15 07:04:08 +00:00
" A = np.zeros((Ns, Ns))\n",
2022-04-03 20:18:44 +00:00
"\n",
" for triangle in elements3:\n",
" M = raideur(coordinates[triangle])\n",
2022-03-15 07:04:08 +00:00
" for i, a in enumerate(triangle):\n",
" for j, b in enumerate(triangle):\n",
" A[a, b] += M[i, j]\n",
2022-04-03 20:18:44 +00:00
" \n",
" return A"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 9,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def second_membre(coordinates, elements3) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Calcul le second membre.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Args:\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" elements3: matrice a trois colonnes. Chaque ligne contient les indices des sommets d'un element triangle, dans le sens antihoraire.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Returns:\n",
" b: vecteur b nécéssaire à la résolution de la formulation variationnelle du problème, sans les conditions de Dirichlet.\n",
" \"\"\"\n",
" Ns = len(coordinates)\n",
2022-03-15 07:04:08 +00:00
" b = np.zeros(Ns)\n",
2022-04-03 20:18:44 +00:00
" for triangle in elements3:\n",
" coords_triangle = coordinates[triangle]\n",
2022-03-15 07:04:08 +00:00
" centre = np.mean(coords_triangle, 0)\n",
" x = coords_triangle[:, 0]\n",
" y = coords_triangle[:, 1]\n",
"\n",
2022-04-03 20:18:44 +00:00
" alpha = calcul_alpha(x, y)\n",
"\n",
" # approximation pour la quadrature du second membre\n",
2022-03-15 07:04:08 +00:00
" b[triangle] += alpha / 6 * f(centre[0], centre[1])\n",
"\n",
2022-04-03 20:18:44 +00:00
" return b"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 10,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def calcul_Ud(coords, dirichlet) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Calcul le vecteur Ud nécéssaire à l'application des conditions de Dirichlet.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Args:\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" dirichlet: vecteur colonne des indices des sommets de la frontiere de Dirichlet.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Returns:\n",
" Ud: vecteur pour appliquer les conditions de Dirichlet.\n",
" \"\"\"\n",
2022-03-15 07:04:08 +00:00
" Ns = len(coords)\n",
" U = np.zeros(Ns)\n",
2022-04-03 20:18:44 +00:00
"\n",
2022-03-15 07:04:08 +00:00
" U[dirichlet.T] = u_d(coords[dirichlet, 0], coords[dirichlet, 1])\n",
"\n",
2022-04-03 20:18:44 +00:00
" return U"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 11,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
"def tildage(A, b, coordinates, dirichlet):\n",
" \"\"\"Permet de retirer les parties de A et b soumis au conditions de Dirichlet, nécéssaire avant la résolution numérique.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Args:\n",
" A: La matrice A de la résolution numérique.\n",
" b: Le vecteur b de la résolution numérique.\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" dirichlet: vecteur colonne des indices des sommets de la frontiere de Dirichlet.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Returns:\n",
" A: La matrice A de la résolution numérique tildée.\n",
" b: Le vecteur b de la résolution numérique tildé.\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation non soumis ausx conditions de Dirichlet. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" \"\"\"\n",
2022-03-15 07:04:08 +00:00
" A_tild = np.delete(A, dirichlet, 0)\n",
" A_tild = np.delete(A_tild, dirichlet, 1)\n",
2022-04-03 20:18:44 +00:00
" \n",
2022-03-15 07:04:08 +00:00
" b_tild = np.delete(b, dirichlet, 0)\n",
2022-04-03 20:18:44 +00:00
" \n",
" coords_tild = np.delete(coordinates, dirichlet, 0)\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" return A_tild, b_tild, coords_tild"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 12,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def untildage(x, dirichlet, U_d) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Opération inverse de la fonction tildage, place dans le vecteur x aux coordonnées de dirichlet les valeurs des conditions\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Args:\n",
" x: le vecteur solution trouvé après résolution.\n",
" dirichlet: vecteur colonne des indices des sommets de la frontiere de Dirichlet.\n",
" Ud: vecteur pour appliquer les conditions de Dirichlet.\n",
2022-03-15 07:04:08 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Returns:\n",
" x: le vecteur solution complet, avec les conditions aux bords.\n",
" \"\"\"\n",
2022-03-15 07:04:08 +00:00
" x_untild = np.zeros(U_d.shape[0])\n",
2022-03-16 14:44:07 +00:00
" not_dirichlet = np.setdiff1d(range(U_d.shape[0]), dirichlet)\n",
"\n",
2022-03-15 07:04:08 +00:00
" x_untild[dirichlet] = U_d[dirichlet]\n",
" x_untild[not_dirichlet] = x\n",
"\n",
2022-04-03 20:18:44 +00:00
" return x_untild"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 15,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [
{
"data": {
2023-06-21 18:29:07 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABE0AAARdCAYAAACkSby6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9eZQkWXoWeD/3mpmv4eERkXtGrkVVV3VXdXdVVmZJgGAAgYSaUTNsghlJoNHyoUZw0CBp0IaQdAQIBiGBgOF80ELfqEECCQmNkNRCahiJEa1eKLpyXyL3tbIywyMjfLPl3vv9YW4W7h4ekRGZvpi5Pb9z8nRXbG7m5m5m9/H7vlcYY0BERERERERERL3kpDeAiIiIiIiIiCiJGJoQEREREREREQ3A0ISIiIiIiIiIaACGJkREREREREREAzA0ISIiIiIiIiIagKEJEREREREREdEADE2IiIgyTAjxB4QQd57j979bCPEvhrlNz0MI8XVCiP93CH/nhhDiDw/4+ncIIX5KCMF7KCIiogywJ70BRERElA5CiD8A4BPGmEPR14wxf3tiGzRmQoivAPAmgK82xuhJbw8RERGNHkMTIiIiom0wxvwqgF+d9HYQERHR+HBqKRER0ZQQQvx1IcRdIcSaEOKSEOJLO1/PCyF+TAhxr/Pvx4QQ+U3+hhFCvNj13z8phPghIUQZYWBwUAhR7/w7KIT4fiHEJ7p+/qNCiHNCiBUhxP8jhHh/1/duCCG+XQhxWgjxRAjxb4QQhS3255uEEBc6+3NeCHGi8/XvFEJc7fr6n9jib7wqhPh1IcSyEOJdIcR3d+9X189tWqYkhJBdj/lYCPFvhRALXd//YiHEf+3s8zudGTlEREQ0BRiaEBERTQEhxMsA/jKAU8aYCoAvB3Cj8+3vAfDFAF4H8GEAbwH43p38fWNMA8BXALhnjJnp/LvXtw3vA/DTAL4VwB4AvwLgl4QQua4f+yoAfxTAcQAfAvB1m+zPnwHw/QD+PIBZAB8F8Ljz7asAfh+AKoAfAPAJIcSBAX+jAuA3AHwSwEEALwL41E72u+OvAPifAPwPnb9TA/BPOo+xCOCXAfwQgAUA3w7g3wkh9jzD4xAREVHCMDQhIiKaDgpAHsAHhBCOMeaGMeZq53tfDeAHjTEPjTHvIQwavnYE2/BnAfyyMebXjTE+gL8PoAjg93T9zD8yxtwzxiwD+CWEQc4g3wjg7xljPmdCS8aYmwBgjPnZzt/Qxph/A+AKwiCo3/8I4IEx5keMMW1jzJox5jPPsF/fDOB7jDF3jDEuwjDnTwshbABfA+BXjDG/0tmeXwfweQAfeYbHISIiooRhaEJERDQFjDFLCGd4fD+Ah0KInxFCHOx8+yCAm10/frPztWHreZxOs9TbABa7fuZB1/9vApjZ5G8dRjijZAMhxJ8XQnyhUw6zAuA1ALt38jd26CiAX+h6vAsIQ6p9ne/9meh7ne9/CYANM1+IiIgofRiaEBERTQljzL82xnwJwoG8AfB3O9+61/la5Ejna4M0AZS6/nt/90M8ZRN6HkcIIRAGF3efuvEb3Qbwu/q/KIQ4CuCfIyxF2mWMmQNwFoDY5G+8sMnfb2Dz/Rz0d77CGDPX9a9gjLnb+d5P9X2vbIz54aftIBERESUfQxMiIqIpIIR4WQjxhzoNXtsAWgCiZXF/GsD3CiH2CCF2A/g+AJ/Y5E99AcD/IoSwhBB/FGEfj8i7AHYJIaqb/O6/BfDHhBBfKoRwAHwbABfAf32GXfoXAL5dCPGmCL3YCUzKCMOb9zr7/b8inGkyyH8AcEAI8a2dZrgVIcQXde3nR4QQC0KI/Qhn6WzmnwH4W53HR+d5/OOd730CwFcKIb6885wVOk1lD23614iIiCg1GJoQERFNhzyAHwbwCGEJzF4A39X53g8h7LNxGsAZAG93vjbIXwXwlQBWEPZC+ffRN4wxFxEGMNc6pSg9JT7GmEsIe3z8eGc7vhLAVxpjvJ3ujDHmZwH8LQD/GsBaZzsWjDHnAfwIgE8jDHE+COC3N/kbawD+SGc7HiDsffIHO9/+KQDvIGyW+x8B/JstNucfAvi/AfxHIcQagN8B8EWdx7gN4I8D+G6EQc5tAN8B3mMRERFNBWHM02baEhERERERERFlDz8FISIiIiIiIiIagKEJEREREREREdEADE2IiIiIiIiIiAZgaEJERERERERENABDEyIiIiIiIiKiAeynfJ9L6xARERERERHRNBObfYMzTYiIiIiIiIiIBmBoQkREREREREQ0AEMTIiIiIiIiIqIBGJoQEREREREREQ3A0ISIiIiIiIiIaACGJkREREREREREAzA0ISIiIiIiIiIagKEJEREREREREdEADE2IiIiIiIiIiAZgaEJERERERERENABDEyIiIiIiIiKiARiaEBERERERERENwNCEiIiIiIiIiGgAhiZERERERERERAMwNCEiIiIiIiIiGoChCRERERERERHRAAxNiIiIiIiIiIgGYGhCRERERERERDQAQxMiIiIiIiIiogEYmhARERERERERDcDQhIiIiIiIiIhoAIYmREREREREREQDMDQhIiIiIiIiIhqAoQkRERERERER0QAMTYiIiIiIiIiIBmBoQkREREREREQ0AEMTIiIiIiIiIqIBGJoQEREREREREQ3A0ISIiIiIiIiIaACGJkREREREREREAzA0ISIiIiIiIiIagKEJEREREREREdEADE2IiIiIiIiIiAZgaEJERERERERENABDEyIiIiIiIiKiARiaEBERERERERENwNCEiIiIiIiIiGgAhiZERERERERERAMwNCEiIiIiIiIiGoChCRERERERERHRAAxNiIiIiIiIiIgGYGhCRERERERERDQAQxMiIiIiIiIiogEYmhARERERERERDcDQhIiIiIiIiIhoAIYmREREREREREQDMDQhIiIiIiIiIhqAoQkRERERERER0QAMTYiIiIiIiIiIBmBoQkREREREREQ0AEMTIiIiIiIiIqIBGJoQEREREREREQ3A0ISIiIiIiIiIaACGJkREREREREREAzA0ISIiIiIiIiIagKEJEREREREREdEADE2IiIiIiIiIiAZgaEJERERERERENABDEyIiIiIiIiKiARiaEBERERERERENwNCEiIiIiIiIiGgAhiZERERERERERAMwNCEiIiIiIiIiGoChCRERERERERHRAAxNiIiIiIiIiIgGYGhCRERERERERDQAQxMiIiIiIiIiogEYmhARERERERERDcDQhIiIiIiIiIhoAIYmREREREREREQDMDQhIiIiIiIiIhqAoQkRERERERER0QAMTYiIiIiIiIiIBmBoQkREREREREQ0AEMTIiIiIiIiIqIBGJoQEREREREREQ3A0ISIiIiIiIiIaACGJkREREREREREAzA0ISIiIiIiIiIagKEJEREREREREdEADE2IiIiIiIiIiAZgaEJERERERERENIA96Q0gIiIaJa01Wq0WbNuGbduQUkIIMenNIiIiIqIUEMaYrb6/5TeJiIiSyhiDIAgQBAFc142/LoSIAxSGKEREREQEYNObQYYmREQ0dYwx8DwPWmsIIeB5XhyMGGOgtY5/NgpRHMeBZVkMUYiIiIiyh6EJERFlg1IKvu/DGBOHH92hSb9BIYrjOPFMFCEEQxQiIiKi6cbQhIiIplt3OU530BHNOtlu8BGFKNH1UUoJx3HimSgMUYiIiIimDkMTIiKaXlpr+L4fl+N0hxo7DU26RdfI7pkoUsoN5TxERERElGoMTYiIaPoYY+JyHAADZ4E8T2gy6PEAhihEREREU4ahCRERTRdjDHzfh1Jqy5KZYYYmg/42wBCFiIiIKOUYmhAR0fTQWsPzvLjZ61aBSBSujIMxJv4XiXqi2LbNEIWIiIgomRiaEBFR+nWX42y3Ies4Q5NBj71ViBKtzkNEREREE8XQhIiI0i0qsxnU7PVpvzep0KRfFKBcvXoVCwsLmJubg2VZPTNRGKIQERERjd2mN2D2OLeCiIjoWeykHCfJom2P/kkpobVGu92Of4YhChE
2022-04-03 20:18:44 +00:00
"text/plain": [
2023-06-21 18:29:07 +00:00
"<Figure size 1440x1440 with 1 Axes>"
2022-04-03 20:18:44 +00:00
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
2023-06-21 18:29:07 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABE0AAARdCAYAAACkSby6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9eZgk214Xen/XiogcKjMrq6vH6rlh77PHs/c5e3e3oN5zlaOocD3KfXFARI8MFxB94b3oI6Ai4HG6iqByr/go4ICKoiIigkwiXjyeCc+unquru6vn3r17Z1ZX5RTDWuv9IzKiM7Oyqqu6KzMjMr6f59nP3ruGzIiMrIhY31y/3xLGGBARERERERERUT856Q0gIiIiIiIiIkoihiZEREREREREREMwNCEiIiIiIiIiGoKhCRERERERERHREAxNiIiIiIiIiIiGYGhCRERERERERDQEQxMiIqIpJYR4UQixKIQ4tQuP9atCiK/fje3a7efcjf0UQnyXEOIfPevvExER0XRiaEJERDSFhBBVAP8QwFcaY27s8mN/XAjx/+7g508KIYwQwt7iZ75HCPHjz7Atu7Kfxpi/aowZayhEREREybfpzQsRERGlixDCNsYEAGCMeQzgt012i0ZvN/az93UjIiIi6sWZJkRERCkmhFgRQvw5IcQigKYQwhZCfJEQ4r8LIVaFEO8IIX5bz89/XAhxXQixLoS4IYT46u7X+2Z6bDY7RAjxCoAfBvDFQoiGEGK1+/UvF0L8TyHEmhDithDie3p+7de6/17t/s4XDzzm7wbwXQD+UPf77/R8+4QQ4te72/sLQoh9Pb+31X4eFkL8ByFETQixLIT4hp7vfY8Q4t8IIX5cCLEG4OND9v9rhBA3hRDvCyH+fPd1/h3d7/1jIcQnen72twkh7gw8978VQrzXfY3/v0MPHhERESUeQxMiIqL0+yoAXw5gDsBBAD8L4BMA5gH8GQD/VgixXwhRAvB3AfweY0wFwG8G8PmdPJEx5hKAbwLwSWNM2Rgz1/1WE8Af627DlwP4ZiHE7+9+7yPdf891f+eTA4/58wD+KoB/1f3+mz3f/iMA/gSAAwBy3f2BEOLIZvvZ/b2fAHAHwGEAXwngrwohvqTncX8fgH/T3d5/3rs9QohXAfx9AF/T/f29AI5u5/URQkgAPwPgHQBHAHwUwLcJIX7Xdn6fiIiIkoWhCRERUfr9XWPMbWNMG8AfBfCfjDH/yRijjTG/COCzAL6s+7MawOtCiKIx5r4x5sJubIAx5leNMee6z7kI4F8C+F934aF/zBiz1N23fw3gQ92vb7qfQohjAH4LgD9njOkYYz4P4B8hDHUinzTG/Pvu77YHnvMrAfxHY8yvGWNcAH8R4eu2HWcA7DfGfJ8xxjPGXEfYc+UP73jPiYiIaOIYmhAREaXf7Z7/PgHgD3RLVla75TO/FcCCMaYJ4A8hnClyXwjxs0KIl3djA4QQv0kI8V+6JSmPu8+x72m/tw0Pev67BaDc/e9N9xPh7JCaMWa953dvIpz5Eel9zQYd7v1+93V7f5vbewLA4YHt+i6EM4CIiIgoZdgIloiIKP1Mz3/fBvDPjDHfMPQHjfnPAP6zEKKIsLTlHwL4XxCW18z0/OihbT5f5F8A+CGEpT8dIcQP4kloMuznt/OYW9l0P7szTeaFEJWe4OQ4gLvbfL77AF7pebwZhCU6ka1eq9sAbhhjXtzWXhAREVGicaYJERHRdPlxAL9XCPG7hBCWEKLQbVR6VAhxUAjx+7q9TVwADTwpO/k8gI8IIY53l/H9zi2e410AR4UQuZ6vVRDO7ugIIc4i7EUSea/7PF/wlMc82e0J8lz7aYy5DeC/A/hr3a+/AeDrur+zHf8GwP8mhPit3X38PvTfM30eYRnQvBDiEIBv6/nepwGsd5vzFrvb9roQ4sw2n5uIiIgShKEJERHRFOkGBr8PYUnIewhnPvxZhNd8CeD/BHAPQA1hz5Fv7v7eLwL4VwAWAXwOwH/c4ml+BcAFAA+EEI+6X/uTAL5PCLEO4LsR9h+JtqkF4K8A+PVuycoXDXnMn+z++30hxG88534CYXPck919/SkAf8kY80tPe9zuY18A8C0IZ8/cB1BH2FQ28s8QNnpdAfALCF+36HcVgP8NYe+VGwAeIeynUt3OcxMREVGyCGN2OhuWiIiIKFuEECsAvn67wQsRERFNB840ISIiIiIiIiIagqEJEREREREREdEQLM8hIiIiIiIiIhqCM02IiIiIiIiIiIZgaEJERERERERENIT9lO+zdoeIiIiIiIiIppnY7BucaUJERERERERENARDEyIiIiIiIiKiIRiaEBERERERERENwdCEiIiIiIiIiGgIhiZEREREREREREMwNCEiIiIiIiIiGoKhCRERERERERHREAxNiIiIiIiIiIiGYGhCRERERERERDQEQxMiIiIiIiIioiEYmhARERERERERDcHQhIiIiIiIiIhoCIYmRERERERERERDMDQhIiIiIiIiIhqCoQkRERERERER0RAMTYiIiIiIiIiIhmBoQkREREREREQ0BEMTIiIiIiIiIqIhGJoQEREREREREQ3B0ISIiIiIiIiIaAiGJkREREREREREQzA0ISIiIiIiIiIagqEJEREREREREdEQDE2IiIiIiIiIiIZgaEJERERERERENARDEyIiIiIiIiKiIRiaEBERERERERENwdCEiIiIiIiIiGgIhiZEREREREREREMwNCEiIiIiIiIiGoKhCRERERERERHREAxNiIiIiIiIiIiGYGhCRERERERERDQEQxMiIiIiIiIioiEYmhARERERERERDcHQhIiIiIiIiIhoCIYmRERERERERERDMDQhIiIiIiIiIhqCoQkRERERERER0RAMTYiIiIiIiIiIhmBoQkREREREREQ0BEMTIiIiIiIiIqIhGJoQEREREREREQ3B0ISIiIiIiIiIaAiGJkREREREREREQzA0ISIiIiIiIiIagqEJEREREREREdEQDE2IiIiIiIiIiIZgaEJERERERERENARDEyIiIiIiIiKiIRiaEBERERERERENwdCEiIiIiIiIiGgIhiZEREREREREREMwNCEiIiIiIiIiGoKhCRERERERERHREAxNiIiIiIiIiIiGYGhCRERERERERDQEQxMiIiIiIiIioiEYmhARERERERERDcHQhIiIiIiIiIhoCIYmRERERERERERDMDQhIiIiIiIiIhqCoQkRERERERER0RAMTYiIiIiIiIiIhmBoQkREREREREQ0BEMTIiIiIiIiIqIhGJoQEREREREREQ3B0ISIiIiIiIiIaAiGJkREREREREREQzA0ISIiIiIiIiIagqEJEREREREREdEQDE2IiIiIiIiIiIZgaEJERERERERENARDEyIiIiIiIiKiIRiaEBERERERERENwdCEiIiIiIiIiGgIhiZEREREREREREMwNCEiIiIiIiIiGoKhCRERERERERHREAxNiIiIiIiIiIiGYGhCRERERERERDQEQxMiIiIiIiIioiHsSW8AERHRKGmt0W63Yds2bNuGlBJCiElvFhERERGlgDDGbPX9Lb9JRESUVMYYBEGAIAjgum78dSFEHKAwRCEiIiIiAJveDDI0ISKiqWOMged50FpDCAHP8+JgxBgDrXX8s1GI4jgOLMtiiEJERESUPQxNiIgoG5RS8H0fxpg4/OgNTQYNC1Ecx4lnogghGKIQERERTTeGJkRENN16y3F6g45o1sl2g48oRImuj1JKOI4Tz0RhiEJEREQ0dRiaEBHR9NJaw/f9uBynN9TYaWjSK7pG9s5EkVJuKOchIiIiolRjaEJERNPHGBOX4wAYOgvkeUKTYc8HMEQhIiIimjIMTYiIaLoYY+D7PpRSW5bM7GZoMuyxAYYoRERERCnH0ISIiKaH1hqe58XNXrcKRKJwZRyMMfE/kagnim3bDFGIiIiIkomhCRERpV9vOc52G7KOMzQZ9txbhSjR6jxERERENFEMTYiIKN2iMpthzV6f9nuTCk0GRQHKtWvXMD8/j7m5OViW1Tc
2022-04-03 20:18:44 +00:00
"text/plain": [
2023-06-21 18:29:07 +00:00
"<Figure size 1440x1440 with 1 Axes>"
2022-04-03 20:18:44 +00:00
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
2022-03-15 07:04:08 +00:00
"n = 50\n",
"coords, elems3, dirichlet, neumann = maillage_carre(n)\n",
"\n",
2022-04-03 20:18:44 +00:00
"# calcul du premier membre de l'équation\n",
2022-03-15 07:04:08 +00:00
"A = assemblage(coords, elems3)\n",
2022-04-03 20:18:44 +00:00
"\n",
"# calcul du second membre de l'équation\n",
"b = second_membre(coords, elems3)\n",
"\n",
"# calcul du vecteur des conditions de dirichlet\n",
2022-03-15 07:04:08 +00:00
"U_d = calcul_Ud(coords, dirichlet)\n",
2022-04-03 20:18:44 +00:00
"\n",
"# on modifie b pour vérifier les conditions \n",
2022-03-15 07:04:08 +00:00
"b -= np.dot(A, U_d)\n",
"\n",
2022-04-03 20:18:44 +00:00
"# on enlève les conditions aux bords avant résolution\n",
2022-03-15 07:04:08 +00:00
"A_tild, b_tild, coords_tild = tildage(A, b, coords, dirichlet)\n",
"\n",
2022-04-03 20:18:44 +00:00
"# on résoud le système\n",
2022-03-15 07:04:08 +00:00
"x = np.linalg.solve(A_tild, b_tild)\n",
2022-04-03 20:18:44 +00:00
"\n",
"# on remet les conditions aux bords\n",
2022-03-15 07:04:08 +00:00
"x_untild = untildage(x, dirichlet, U_d)\n",
"\n",
2022-04-03 20:18:44 +00:00
"# on affiche le résultat\n",
2022-04-03 20:42:26 +00:00
"show(coords, x_untild, \"solution calculée\")\n",
2022-04-03 20:18:44 +00:00
"\n",
"# on compare avec le résultat théorique exacte\n",
2022-04-03 20:42:26 +00:00
"show(coords, u_ex(coords[:, 0], coords[:, 1]), \"résultat théorique\")"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2022-04-03 20:18:44 +00:00
"## Partie II : maillage mixte et ajoût des conditions de Neumann\n"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 14,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def f(x, y) -> int:\n",
2022-04-03 20:18:44 +00:00
" return 1\n",
"\n",
"\n",
2022-04-03 20:42:26 +00:00
"def u_d(x, y) -> int:\n",
2022-04-03 20:18:44 +00:00
" return 1\n",
"\n",
"\n",
2022-04-03 20:42:26 +00:00
"def g(x) -> int:\n",
2022-04-03 20:18:44 +00:00
" return 1"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 15,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"# Création d'un maillage mixte\n",
"\n",
2022-03-15 07:04:08 +00:00
"e3 = np.array(\n",
" [[1, 2, 12], [2, 3, 12], [3, 4, 14], [4, 5, 14], [2, 15, 3], [3, 15, 4]]\n",
").astype(int)\n",
"\n",
"e4 = np.array(\n",
" [\n",
" [0, 1, 12, 11],\n",
" [11, 12, 13, 10],\n",
" [12, 3, 14, 13],\n",
" [10, 13, 8, 9],\n",
" [13, 14, 7, 8],\n",
" [14, 5, 6, 7],\n",
" ]\n",
").astype(int)\n",
"\n",
"dds = np.array([2, 15, 4, 6, 7, 8, 9, 10, 11, 0]).astype(int)\n",
"\n",
"nns = np.array([[4, 5], [5, 6], [0, 1], [1, 2]]).astype(int)\n",
"\n",
"ccs = np.array(\n",
" [\n",
2022-04-03 20:18:44 +00:00
" [0, 0],\n",
2022-03-15 07:04:08 +00:00
" [1 / 3, 0],\n",
2022-04-03 20:18:44 +00:00
" [16 / 30, 0],\n",
2022-03-15 07:04:08 +00:00
" [2 / 3, 1 / 3],\n",
2022-04-03 20:18:44 +00:00
" [1, 14 / 30],\n",
2022-03-15 07:04:08 +00:00
" [1, 2 / 3],\n",
2022-04-03 20:18:44 +00:00
" [1, 1],\n",
" [2 / 3, 1],\n",
" [1 / 3, 1],\n",
" [0, 1],\n",
" [0, 2 / 3],\n",
" [0, 1 / 3],\n",
2022-03-15 07:04:08 +00:00
" [1 / 3, 1 / 3],\n",
" [1 / 3, 2 / 3],\n",
" [2 / 3, 2 / 3],\n",
2022-04-03 20:18:44 +00:00
" [1, 0],\n",
2022-03-15 07:04:08 +00:00
" ]\n",
")"
]
},
2023-06-21 18:29:07 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Pour calcul plus facilement l'intégrale de notre carré, on change de base de telle manière à ce que la calcul soit trivial.\n",
"\n",
"On cherche la fonction affine $\\phi_Q$ qui transforme notre quadrangle en carré de côté unitaire.\n",
"\n",
"$$\n",
"\\begin{cases}\n",
" x & = (x_2 - x_1) \\xi + (x_4 - x_1) \\zeta + x_1 \\\\\n",
" y & = (y_2 - y_1) \\xi + (y_4 - y_1) \\zeta + y_1\n",
"\\end{cases}\n",
"$$\n",
"\n",
"$$\n",
"\\begin{cases}\n",
" \\phi_1(\\xi, \\zeta) & = (1 - \\xi) (1 - \\zeta) \\\\\n",
" \\phi_2(\\xi, \\zeta) & = \\xi (1 - \\zeta) \\\\\n",
" \\phi_3(\\xi, \\zeta) & = \\xi \\zeta \\\\\n",
" \\phi_4(\\xi, \\zeta) & = (1 - \\xi) \\zeta\n",
"\\end{cases}\n",
"$$\n",
"\n",
"On pose alors:\n",
"\n",
"$$\n",
"J_k =\n",
"\\begin{bmatrix}\n",
" x_2 - x_1 & x_4 - x_1 \\\\\n",
" y_2 - y_1 & y_4 - y_1\n",
"\\end{bmatrix}\n",
"=\n",
"\\begin{bmatrix}\n",
" a & b \\\\\n",
" b & c\n",
"\\end{bmatrix}\n",
"\n",
"\\quad et \\quad \n",
"\n",
"[M^A_T]_{ij} = \\displaystyle \\int_{[0,1]} \\int_{[0,1]} \\phi_i(\\xi, \\zeta) \\ |J_k| \\ \\phi_j(\\xi, \\zeta) \\ d\\xi \\ d\\zeta\n",
"$$\n",
"\n",
"Donc par calcul de toutes ces combinaisons, on obtient:\n",
"\n",
"$$\n",
"M = \n",
"\\begin{bmatrix}\n",
" 2 a + 3 b + 2 c & -2 a + c & -a - 3 b - c & a - 2 c \\\\\n",
" -2 a + c & 2 a - 3 b + 2 c & a - 2 c & -a + 3 b - c \\\\\n",
" -a - 3 b - c & a - 2 c & 2 a + 3 b + 2 c & -2 a + c \\\\\n",
" a - 2 c & -a + 3 b - c & -2 a + c & 2 a - 3 b + 2 c\n",
"\\end{bmatrix}\n",
"$$"
]
},
2022-03-15 07:04:08 +00:00
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 16,
2022-03-15 07:04:08 +00:00
"metadata": {},
2022-03-16 14:44:07 +00:00
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.66666667, -0.16666667, -0.33333333, -0.16666667],\n",
" [-0.16666667, 0.66666667, -0.16666667, -0.33333333],\n",
" [-0.33333333, -0.16666667, 0.66666667, -0.16666667],\n",
" [-0.16666667, -0.33333333, -0.16666667, 0.66666667]])"
]
},
2023-06-21 18:29:07 +00:00
"execution_count": 16,
2022-03-16 14:44:07 +00:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
2022-04-03 20:42:26 +00:00
"def raideur_quadrangle(quadrangle) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Construction de la matrice de raideur ́elementaire relative à un ́élément quadrangle.\n",
"\n",
" Args:\n",
" quadrangle: les coordonnées x et y des quatres points formant le quadrangle.\n",
"\n",
" Returns:\n",
" M: La matrice de raideur ́elementaire.\n",
" \"\"\"\n",
2022-03-16 14:44:07 +00:00
" x = quadrangle[:, 0]\n",
" y = quadrangle[:, 1]\n",
"\n",
2022-04-03 20:18:44 +00:00
" # calcul de la jacobienne et de son déterminant\n",
2022-03-16 14:44:07 +00:00
" J_kk = np.array([[x[1] - x[0], x[3] - x[0]], [y[1] - y[0], y[3] - y[0]]])\n",
" det_J_kk = np.linalg.det(J_kk)\n",
"\n",
" # on récupère les coefficients\n",
" coeffs = np.linalg.inv(np.matmul(J_kk.T, J_kk))\n",
" a = coeffs[0, 0]\n",
" b = coeffs[0, 1]\n",
" c = coeffs[1, 1]\n",
"\n",
" # on calcul M (on a calculé toutes les intégrales au préalable)\n",
" M = np.array(\n",
" [\n",
" [2 * a + 3 * b + 2 * c, -2 * a + c, -a - 3 * b - c, a - 2 * c],\n",
" [-2 * a + c, 2 * a - 3 * b + 2 * c, a - 2 * c, -a + 3 * b - c],\n",
" [-a - 3 * b - c, a - 2 * c, 2 * a + 3 * b + 2 * c, -2 * a + c],\n",
" [a - 2 * c, -a + 3 * b - c, -2 * a + c, 2 * a - 3 * b + 2 * c],\n",
" ]\n",
" )\n",
"\n",
" return det_J_kk / 6 * M\n",
"\n",
2022-04-03 20:18:44 +00:00
"# on affiche la première matrice de raideur pour vérifier\n",
2022-03-16 14:44:07 +00:00
"raideur_quadrangle(ccs[e4[0]])"
]
2022-03-15 07:04:08 +00:00
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 17,
2022-03-15 07:04:08 +00:00
"metadata": {},
2022-04-03 20:18:44 +00:00
"outputs": [],
2022-03-16 14:44:07 +00:00
"source": [
2022-04-03 20:42:26 +00:00
"def assemblage_quadrangle(coordinates, elements4) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Assemblage de la matrice A dans le cas d'un maillage constitué uniquement d'éléments quadrangles.\n",
"\n",
" Args:\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" elements4: matrice a quatre colonnes. Chaque ligne contient les indices des sommets d'un element quadrangle, dans le sens antihoraire.\n",
"\n",
" Returns:\n",
" A: matrice nécéssaire à la résolution de la formulation variationnelle du problème.\n",
" \"\"\"\n",
" Ns = len(coordinates)\n",
2022-03-16 14:44:07 +00:00
" A = np.zeros((Ns, Ns))\n",
2022-04-03 20:18:44 +00:00
"\n",
" for quadrangle in elements4:\n",
" M = raideur_quadrangle(coordinates[quadrangle])\n",
2022-03-16 14:44:07 +00:00
" for i, a in enumerate(quadrangle):\n",
" for j, b in enumerate(quadrangle):\n",
" A[a, b] += M[i, j]\n",
2022-04-03 20:18:44 +00:00
" \n",
" return A"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 18,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def second_membre_quadrangle(coordinates, elements4) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Calcul le second membre.\n",
2022-03-16 14:44:07 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Args:\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" elements4: matrice a quatre colonnes. Chaque ligne contient les indices des sommets d'un element quadrangle, dans le sens antihoraire.\n",
"\n",
" Returns:\n",
" b: vecteur b nécéssaire à la résolution de la formulation variationnelle du problème, sans les conditions de Dirichlet.\n",
" \"\"\"\n",
" Ns = len(coordinates)\n",
2022-03-16 14:44:07 +00:00
" b = np.zeros(Ns)\n",
2022-04-03 20:18:44 +00:00
" for quadrangle in elements4:\n",
" coords_quadrangle = coordinates[quadrangle]\n",
" centre = np.mean(coords_quadrangle, 0)\n",
" x = coords_quadrangle[:, 0]\n",
" y = coords_quadrangle[:, 1]\n",
2022-03-16 14:44:07 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" alpha = calcul_alpha(x, y)\n",
2022-03-16 14:44:07 +00:00
"\n",
" b[quadrangle] += alpha / 4 * f(centre[0], centre[1])\n",
"\n",
2022-04-03 20:18:44 +00:00
" return b"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 19,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def condition_neumann(coordinates, neumann) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" \"\"\"Calcul le vecteur nécéssaire à l'application des conditions de Neumann.\n",
2022-03-16 14:44:07 +00:00
"\n",
2022-04-03 20:18:44 +00:00
" Args:\n",
" coordinates: matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" neumann: vecteur colonne des indices des sommets de la frontiere de Neumann.\n",
"\n",
" Returns:\n",
" Ud: vecteur pour appliquer les conditions de Neumann.\n",
" \"\"\"\n",
" Ns = len(coordinates)\n",
2022-03-16 14:44:07 +00:00
" coeffs = np.zeros(Ns)\n",
2022-04-03 20:18:44 +00:00
" for i, j in neumann:\n",
" point1 = coordinates[i]\n",
" point2 = coordinates[j]\n",
2022-03-16 14:44:07 +00:00
" \n",
" valeur = np.linalg.norm(point1 - point2) / 2 * g((point1 + point2) / 2)\n",
" coeffs[i] += valeur\n",
" coeffs[j] += valeur\n",
"\n",
2022-04-03 20:18:44 +00:00
" return coeffs"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 20,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [
{
"data": {
2022-04-03 20:42:26 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAECCAYAAADXWsr9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABj/ElEQVR4nO2dd3hc1bW33z0zGvXei1VsWe6yZVu2gYApocQUOySAQ0InBEJJAmkk+XIJuSHcJNw0knAvkJtGhwAOxYTQIeCqbtmS1Xvvmj77+0OaYTSacmY0I8to3ucx2NKZc/bMnN/Za6+19lpCSkmIECE+2ahO9ABChAgRfEJCDxFiERASeogQi4CQ0EOEWASEhB4ixCIgJPQQIRYBIaFPI4RYLoSoFEIUBOBcbwshbgzEuBYSQohvCSH+KoRwe98IIXKFEONCCPV8ji2EZ0JCB4QQ8cDDwOellE0BPve1Qoj3fTg+XwghhRCaQI5jrgghPgNsAq6VUlrdHSelbJVSxkgpLfM3uhDeWFA303wihNBIKc0AUsoR4MwTO6KFjZTyVeBVT8c4fqYhFhaLakYXQjQLIb4jhKgEJoQQGiHENiHEv4UQw0KICiHEmQ7HXyuEaBRCjAkhmoQQX5z++T1CiL85HOdyFhZCrAIeAk6ZNmeHp39+oRCiTAgxKoRoE0Lc4/Cyd6f/Pzz9mlNcvI97hBBPCyH+Mj22GiHEZoffSyFEocO//ySE+M/pv58phGgXQnxbCNErhOgSQuwSQuwQQtQJIQaFEN9zeK1KCPFdIUSDEGJg+rpJTu/7BiFEK/Cm82chhCgQQrwzPc7XhRAP2j4721hcfEef9nbtEL6xqIQ+zReAC4EEIB14GfhPIAn4JvCcECJVCBEN/Ab4jJQyFjgVKPflQlLKWuBm4MNpczZh+lcTwNXTY7gQuEUIsWv6d2dM/z9h+jUfujn9JcCT0+fYAzzow9AygAggG/ghU8uWLzFlmp8O/D8HX8XtwC5gO5AFDAG/czrfdmAVcL6Laz0OHAJSgB8D1/gwTiXXDqGAxSj030gp26SUOqZu7leklK9IKa1SyteBg8CO6WOtwFohRKSUsktKWROIAUgp35ZSVk1fsxJ4gqmb2Rfenx63BfgrsN6H15qAn0gpTUw9LFKAX0spx6bf4xGH890MfF9K2S6lNAD3AJ93sl7ukVJOTH+mdoQQuUAp8P+klAYp5bvAP3wYp5Jrh1DAYhR6m8Pf84DLps324WnT+lNAppRyAriCqZutSwjxshBiZSAGIITYKoR4SwjRJ4QYmb5Gio+n6Xb4+yQQ4YMABhycZTZx9jj8XgfETP89D3je4fOpBSxMWUM2HD9TR7KAoenP0kaLwjEqvXYIBSxGoTtu12sD/iqlTHD4Ey2lvB9ASvmalPJcIBM4ypSJC1Omd5TDeTIUXs/G40yZ20uklPFMreOFh+N9ZdKH8Xmjjanli+NnFCGl7HA4xt2Yu4DE6WWQjVyHv8/4HKdDcqk+XjuEAhaj0B35G3CxEOJ8IYRaCBEx7SDKEUKkCyF2Tt+kBmCcKVMeptbqZ0zHjOOBuz1cowfIEUJoHX4WCwxKKfVCiC3AlQ6/65u+ztI5vK9y4Mrp93QBvi8LHHkI+IkQIg9g2n+xU8kLpZQtTC2FfiSE0AohPgVc7HBIHVOWyIVCiDDgB0B4IK4dYiaLWuhSyjZgJ/A9pgTWBnyLqc9FBdwJdAKDTInllunXvQ48BVQy5Wh6ycNl3gRqgG4hRP/0z74K3CuEGGPKGfa0w5gmgZ8AH0ybrNv8eGtfY0pQw8AXgRf8OIeNXzNlffxzerwfAVt9eP2V08cPAv8B/MX2i+mw5leBR4AOpmZ4Ry/8XK8dYhoRKjwRYj6ZDiUWSim/dKLHsphY1DN6iBCLhZDQQ4RYBIRM9xAhFgGhGT1EiEVASOghQiwCvGVShez6ECGCj/B+yNwIzeghQiwCQkIPEWIREBJ6iBCLgJDQQ4RYBISEHiLEIiAk9BAhFgEhoYcIsQgICT1EiEVASOghQiwCQkIPEWIREBJ6iBCLgJDQQ4RYBISEHiLEIiAk9BAhFgEhoYcIsQgICf0EIKXEaDRiNpsJlfIKMR+EeljNM1arFaPRiF6vt/9MrVYTFhaGRqNBrVYjRNDrEIRYZHgrDhmabgKElBKz2YzZbEYIgclksv9cSonVarUL3GAwEBsbi1arDQl/cRD0Lzg0o88DNlPdUcw2hBAIIVCpVPZjGxoayM/PJypqqi1ZaMYPMVdCQg8yZrOZ9vZ2LBYL2dnZCCHss7grwdqEr1arUavV9tlep9PZj9doNPY/IeGHUEJI6EHC0VS3Wq12k91XXM34FosFs9lsP0aj0dhnfJVKFRJ+iFmEhB4ErFYrJpPJbqrbZnGleDredj4bzsIXQsyY8UPCDwEhoQcUm+hsjjbbLOxOuO7Md19wJXyz2Wwfg80a0Gg0aLXakPAXKSGhBwgpJSaTCYvFMkt8zkL3Nsv7agE4v9ZZ+O3tU52IMzMzQzP+IiUk9ABgi43bZmhXnvUTlRjjOB6bc89kMs2Y8W1rfLVaHRL+J5SQ0OeAc2zcZqo7E8g1+lyxefRtuBK+zbGn0WhcPrhCnHyEhO4nzrFxT2LwR7iBFrq78bkSvtFoxGAwAFN+hrCwMPuMHxL+yUlI6H5gc7i5M9Wd8WdGP1F4Er6jY8/R1A+x8Al9Sz5gM3PLy8sxGAyK17Mn2nSfy7lswrc57gCMRiMHDhxgeHiY0dFRJicn7eHEEAuT0IyuEMfYuM3xppQT6YwLJI5OPb1eb5/5jUYjRqMRIDTjL1BCQveCc2zcZr4GU+gny4PBMV0XPt6g4yx8xzz9kPBPDCGhe8BdbFwI4ZOZerII11fcbdCxYRO+wWCY5dwLCX9+CQndDZ5i48GeoT8pDwZvwpdSzjDzbeG8EIEnJHQnHE11d7FxlUoVmtH9wJXwrVarvQhHZ2cnOTk5aLXa0M68ABMSugNKY+OBmKGDlQJ7MuH8GXd3d5OVlRWqvhMEQkKfxlsaqyPBdsYtZpxj+M578UPC949FL3SlaayOBNsZF3owTOFqL35I+P6xqIVutVrp7e1FSkliYqLim+RkFO4nQQBKhB+qvuOaRSl0R4fb2NgYUkqSkpIUv36uzrjh4WGqq6tRqVQkJiaSmJhIfHy83WxdSJlxCxl31XdsjlQbWq2W8PDwRb0zb9EJ3dlUV6vV9mQYpfg7o0spaWlpoaenh+LiYtRqNSMjI/T399PQ0IBGoyExMdEeegrhG66E39raSlhYGGlpaTO25C62vfiLSuiuSjz5OjuDf844i8VCeXk5ERERlJaWYrFYsFqtpKSkkJKSAkyVeR4aGqK3t5ehoSG6urpISkoiMTGRqKioRXNTBgrbA9Zmyi/msluLQuieYuO+ihZ8d8aNj48zMjLCunXryMjIAKZ2wDkTHh5ORkYGer2eyMhIYmNjGRoaorGxEZ1OR0xMDImJiSQlJREeHu7TmBcrVqt1Rkkvb2W3PqnC/8QL3VOJJ/BdtLbXKHk4SClpa2ujvb2dmJgYu8iVXiMqKoqoqCiys7ORUjI+Ps7g4CC1tbWYTCbi4uJISkoiISGBsLAwn97DYsFR6M64Ev4ntfrOJ1roSmLj/pru3l5jNpupqalBrVazadMmysvLFZ/fXYJNbGwssbGx5OXlYbVaGRkZYWhoiNbWVnvkwNmxd7IQLJ+EJ6E740v1nZNN+J9IofsSG/dH6N5m9LGxMaqqqsjPzycrKwuLxRLwG9nRYw9TD5bh4eFZjr2kpKSTwrEXiIq4rvBF6M54KsLR2dlJeno6UVFRJ0XZrU+c0H0p8QT+C93VGhugvb2dtrY2iouLiYmJsR/vLDZP1/UnvKbRaFw69jo6OhgYGECr1WKxWBasY28hCt0ZR+EPDg6Snp4+o/qObcZfiHvxP1FCt5lZtg9dyY0TKK+7xWLhyJEjSCkpLS1Fo/n4o3Ul3GDPsjbHXkZGht20F0LQ2NjI5OQksbGxC8qxF0hBzsd5LRaLvW0WfPx92vbi33bbbfzgBz9g5cqVAb+2P3w
2022-04-03 20:18:44 +00:00
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# calcul de premier membre de l'équation\n",
2022-03-16 14:44:07 +00:00
"A3 = assemblage(ccs, e3)\n",
"A4 = assemblage_quadrangle(ccs, e4)\n",
"A = A3 + A4\n",
"\n",
2022-04-03 20:18:44 +00:00
"# calcul du second membre de l'équation\n",
"b3 = second_membre(ccs, e3)\n",
"b4 = second_membre_quadrangle(ccs, e4)\n",
"b = b3 + b4\n",
2022-03-16 14:44:07 +00:00
"\n",
2022-04-03 20:18:44 +00:00
"# calcul de Ud pour les conditions de Dirichlet\n",
2022-03-16 14:44:07 +00:00
"U_d = calcul_Ud(ccs, dds)\n",
2022-04-03 20:18:44 +00:00
"\n",
"# modifiction de b pour vérifier Dirichlet\n",
2022-03-16 14:44:07 +00:00
"b -= np.dot(A, U_d)\n",
"\n",
2022-04-03 20:18:44 +00:00
"# modification de b pour vérifier Neumann\n",
"b += condition_neumann(ccs, nns)\n",
2022-03-16 14:44:07 +00:00
"\n",
2022-04-03 20:18:44 +00:00
"# on enlève les conditions aux bords (Dirichlet) avant résolution\n",
2022-03-16 14:44:07 +00:00
"A_tild, b_tild, ccs_tild = tildage(A, b, ccs, dds)\n",
"\n",
2022-04-03 20:18:44 +00:00
"# on résoud le système\n",
2022-03-16 14:44:07 +00:00
"x = np.linalg.solve(A_tild, b_tild)\n",
2022-04-03 20:18:44 +00:00
"\n",
"# on remet les conditions aux bords (Dirichlet)\n",
2022-03-16 14:44:07 +00:00
"x_untild = untildage(x, dds, U_d)\n",
"\n",
2022-04-03 20:18:44 +00:00
"# on affiche le résultat\n",
2022-04-03 20:42:26 +00:00
"show(ccs, x_untild, \"résultat numérique\")"
2022-04-03 20:18:44 +00:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compléments d'analyse du système\n"
2022-03-16 14:44:07 +00:00
]
2022-03-24 07:10:59 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2022-04-03 20:18:44 +00:00
"### Analyse de l’
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 21,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:42:26 +00:00
"def f(x, y) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\n",
"\n",
"\n",
2022-04-03 20:42:26 +00:00
"def u_ex(x, y) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" return np.sin(np.pi * x) * np.sin(np.pi * y)\n",
"\n",
"\n",
2022-04-03 20:42:26 +00:00
"def u_d(x, y) -> np.ndarray:\n",
2022-04-03 20:18:44 +00:00
" return np.zeros(x.shape[0])"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 22,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"98/98\r"
]
}
],
"source": [
"erreurs = []\n",
"hs = []\n",
"range_n = range(3, 100, 5)\n",
"\n",
"for n in range_n:\n",
" print(f\"{n}/{max(range_n)}\", end=\"\\r\")\n",
" coords, elems3, dirichlet, neumann = maillage_carre(n)\n",
"\n",
" A = assemblage(coords, elems3)\n",
" b = second_membre(coords, elems3)\n",
" U_d = calcul_Ud(coords, dirichlet)\n",
" b -= np.dot(A, U_d)\n",
"\n",
" A_tild, b_tild, coords_tild = tildage(A, b, coords, dirichlet)\n",
"\n",
" x = np.linalg.solve(A_tild, b_tild)\n",
" x_untild = untildage(x, dirichlet, U_d)\n",
" x_ex = u_ex(coords[:, 0], coords[:, 1])\n",
"\n",
" v = x_untild - x_ex\n",
" h = np.sqrt(1/len(v))\n",
" hs.append(h)\n",
" erreur = h * np.linalg.norm(v)\n",
" erreurs.append(erreur)"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 23,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"regression linéaire: \n",
"2.078 x - 0.0363\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy7klEQVR4nO3dd3xUZfr//9edBAgpdAgEAoQOSWgJRUBJBGkqCCqo2N1FXetH0V2X3Y/rR3d/rsvadf2yrmKPUqRJLwHpJEAgQHqjBULPhNSZ+/dHBpZoyiTM5MxJrufjMQ+nnDnnzQ3myjn3OddRWmuEEEKIKzyMDiCEEMK9SGEQQghRjhQGIYQQ5UhhEEIIUY4UBiGEEOV4GR2gJtq0aaO7du3qknXn5+fj6+vrknU7m5myguR1JTNlBXPlNVNWqDpvXFzcGa11W4dXprU2zSM8PFy7yqZNm1y2bmczU1atJa8rmSmr1ubKa6asWledF4jVNfhZa8ihJKXUP5RSiUqpA0qpH5VSLYzIIYQQ4teMmmNYB4RqrfsDycArBuUQQgjxC4YUBq31Wq11qf3lTqCTETmEEEL8mtIGt8RQSi0Hvtdaf13J57OAWQABAQHh0dHRv/wcX19fPD09ryuH1hql1HWto67UdVar1Up+fj61/bdisVjw8/NzcirXMVNeM2UFc+U1U1aoOm9UVFSc1jrC4ZXVZEKiJg9gPZBQwWPKNcvMAX7EXqCqe1Q0+Zyenq5zc3O1zWZzaIKmMpcuXbqu79elusxqs9l0bm6uTk9Pr/U66tMknrsxU1atzZXXTFm1du7ks8tOV9Vaj63qc6XUw8BtwBh78FopLCyka9eupvlt32yUUrRu3Zrc3Fyjowgh6ogh1zEopSYALwOjtdaXnbC+6w8lKiXjK0TDYtRZSR8C/sA6pdR+pdQnBuUQQgi3duFyMX9Zdoi8wpI626ZRZyX10FoHaa0H2h9PGJHDLPLy8vjXv/5V68lfIYQ5bUnOZfy7W/h6Zxa70s/V2XalV5KbKy4u5qmnnmL06NGVHtJZtmwZb775Zh0nE0K4yuXiUv53aQIPfrabZt6NWPLUSMb2C6iz7ZuqV1JD1LhxY7788stKPy8tLWXy5MlMnjy5DlMJIVxlX/Z5Xvghnowz+Tw2KpiXxvfGu/hCnWaoV4XhteWHOHziUq2+a7VaK7wWol9gM169PaTK72ZmZjJx4kRGjRrF9u3b6dixI0uXLmXixInMnTuXiIgIzpw5Q0REBJmZmcyfP58lS5aQn59PSkoKs2fPpri4mK+++oomTZqwcuVKWrVqRVpaGk899RS5ubn4+Pjw73//mz59+vDEE0/g7+/Pvn37GDlyJP379yc2NpYPP/yQjIwM7rvvPiwWC1OmTOHdd9/FYrEQExPD3LlzWbFiBQBPP/00ERERPPzww8TFxfHCCy9gsVho06YN8+fPp0OHDrUaRyFE7ZRYbby/IYWPNqXSoXlTvv3tMEZ08YeMjXA8DkLvhDY96ySLHEpykpSUFJ566ikOHTpEixYtWLRoUZXLJyQksHjxYvbs2cOcOXPw8fFh37593HDDDVf3EGbNmsUHH3xAXFwcc+fO5Xe/+93V7x87dozt27fz9ttvl1vvc889x5NPPsnBgwcd+uFeUlLCM888w8KFC4mLi+PRRx9lzpw5tRgBIURtpZzKY+rH2/hgYypTB3Vi1fM3MqL5BdjzaVlRCBwMLTrXWZ56tcdQ3W/2VcnLy8Pf37/W3w8ODmbgwIEAhIeHk5mZWeXyUVFR+Pv74+/vT/Pmzbn99tsBCAsL48CBA1gsFrZv387dd9999TtFRUVXn999990V7uFs27btalF64IEH+P3vf19ljqSkJBISErjllluAsj0n2VsQom7YbJrPt2fy99WJ+DXx4pP7w5nQyx9SfoLTR8C3DQy6H5rXbdegelUYjNSkSZOrzz09PSkoKMDLywubzQaUXYhX2fIeHh5XX3t4eFBaWorNZqNFixbs37+/wu1V1Se+oknqa7Ncm0drTUhICDt27KjmTyiEcKbjFwqY/UM8O9LPMrZvO/6/qWG0zU+G3d+BtQSCb4Sg4eBZ9z+m5VCSC3Xt2pW4uDgAFi5cWKPvNmvWjODgYBYsWACU/QCPj4+v9nsjR47kSj+pb7755ur7Xbp04fDhwxQVFXHhwgU2bNgAQO/evcnNzb1aGEpKSjh06FCNsgohHKe1ZlHcMSa8s4UDxy7w9zvD+Pfd3Wmb/iMk/gQ+bSDiUeg6ypCiAFIYXGr27Nn861//YtCgQZw5c6bG3//mm2/4z3/+w4ABAwgJCWHp0qXVfue9997jo48+IiwsjOPHj199PygoiOnTpxMaGsr06dMZNGgQUHbW08KFC/n973/PgAEDGDhwINu3b69xViFE9c5ainjy6728uCCePh38WfXsKGYEnEDF/gfyTkCv8WWHjnzbGBu0Jo2VjH5U1ETv8OHDlTaOqon62kTP19fXKdu8nnGuT83I3I2ZsmptrrzOzrr+cI4Of32d7vnHlfqTmFRdeuG41rs/1Xrj37Q+sEDrgovXtX5TNNETQggBlqJS3lhxmOg9R+nT3p+vHh5I38J42LcHGvlA6DRo29vomOVIYajnLBaL0RGEaLD2ZJ7jhR/2c/x8AU9Gduf5wV40SV8ABRcgcCB0i4JG3kbH/BUpDEII4WRFpVbeXpfMvC3pBLX0YeFjAxlcshcOJYBPKxg0s06vS6gpKQxCCOFER05e4n++309iTh73DunEn4eCT9YCKC2CLiOgy0jDzjZylHunE0IIk7DaNPO2pPP2uiSaN23Ml/f14ia9F1LToVkg9J4Ifu2MjukQKQxCCHGdss9e5sUF+9mTeZ5JIe14c1gxzXKWAAp63lLW0sLDPFcHmCepcNiJEye46667jI4hRL2ntSZ6dzYT39tCYk4eH08J5KO+CTQ7vgVadIWhv4VOEaYqCiB7DE535Txgj1r+QygtLcXL6/r+WgIDA2t8pbUQomZO5xXyyqKDbEg8zY3dmvPO8HzanFsLNm/oNwXa9QWT3ha3fhWGlPVgOVWrr3oVFoJ3BaeN+QVAz7FVfjczM5Px48czbNgw4uLimD59OitWrKCoqIipU6fy2muvAfD666/z9ddf07ZtW4KCgggPD2f27NlERkYycOBAtm7dyr333ktkZGSFbbDff/99PvnkEzw8PAgNDSU6OprNmzfz3HPPAWU9krZs2cLZs2e57bbbSEhIoLCwkCeffJLY2Fi8vLx4++23iYqKYv78+SxbtozLly+TlpbG1KlTeeutt2o1dkI0NKsTTvLK4oNcLrby1pjm3OV3EI+zF6BDf+h+MzRqanTE62JIYVBKvQ5MAWzAaeBhrfUJI7I4S0pKCl988QWXLl1i4cKF7N69G601kydPZsuWLTRt2pRFixYRHx9PSUkJgwcPJjw8/Or3i4uLiY2NpaSkhNGjR7N06VLatm3L999/z5w5c/jss8948803ycjIoLi4GKvVCsDcuXP56KOPGDlyJBaLBe9fFLePPvoIpRQHDx4kMTGRcePGkZycDMD+/fvZt28fTZo0oXfv3jzzzDMEBQXV3aAJYTKXCkv4y9JDLN53nPBAbz4cdokOl/eDagkD7oFWwUZHdAqj9hj+obX+M4BS6lngf4Hrv+9zNb/ZV6U0Lw+uo+12ly5dGD58OLNnz2bt2rVXexFZLBZSUlLIy8tjypQpeHt74+3tfbXN9hUzZswAqm6D3b9/f2bOnMn48eO59957gbKmeS+88AIzZ85k2rRpdOpUvj3v1q1beeaZZwDo06cPXbp0uVoYxowZQ/PmzQHo168fWVlZUhiEqMT21DPMXhDPqbxCXrvBg/tbJ+BZUAidh9sb3jUyOqLTGFIYtNbX3mbNFzD9Xe6vtMHWWvPKK6/w+OOPl/v83Xffdfj7lbXB/umnn9iyZQuLFi3i7bff5uDBg/zhD3/g1ltvZeXKlYwcOZI1a9b8aq+hMr9sFV5aWurQ94RoSApLrLy1OonPtmUQ1hq+HXeGrvoEeLeH3jPAv73REZ3OsDkGpdRfgQe
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"log_hs = np.log(hs)\n",
"log_erreurs = np.log(erreurs)\n",
"\n",
"# affichage des erreurs\n",
"plt.plot(log_hs, log_erreurs, label=\"numérique\")\n",
"\n",
"# regression linéaire\n",
"coeffs = np.polyfit(log_hs, log_erreurs, 1)\n",
"reg = np.poly1d(coeffs)\n",
"\n",
"# affichage de la regression linéaire\n",
"plt.plot(log_hs, reg(log_hs), alpha=0.5, label=\"regression\")\n",
"plt.xlabel(\"log(h)\")\n",
"plt.ylabel(\"log(erreur)\")\n",
"plt.grid()\n",
"plt.legend()\n",
"\n",
"print(f\"regression linéaire: {reg}\")"
]
},
2022-04-03 20:42:26 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On observe alors que l'on obtient un ordre de discrétisation d'environ 2."
]
},
2022-04-03 20:18:44 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Résolution du système linéaire par une méthode directe\n"
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 24,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"98/98\r"
]
}
],
"source": [
"zeros_A = []\n",
"zeros_L = []\n",
"range_n = range(3, 100, 5)\n",
"\n",
"for n in range_n:\n",
" print(f\"{n}/{max(range_n)}\", end=\"\\r\")\n",
" coords, elems3, dirichlet, neumann = maillage_carre(n)\n",
"\n",
" A = assemblage(coords, elems3)\n",
"\n",
" A_tild = np.delete(A, dirichlet, 0)\n",
" A_tild = np.delete(A_tild, dirichlet, 1)\n",
"\n",
" L = np.linalg.cholesky(A_tild)\n",
"\n",
2023-06-21 18:29:07 +00:00
" zeros_A.append(len(np.where(A != 0)[0]))\n",
" zeros_L.append(len(np.where(L != 0)[0]))"
2022-04-03 20:18:44 +00:00
]
},
{
"cell_type": "code",
2023-06-21 18:29:07 +00:00
"execution_count": 25,
2022-04-03 20:18:44 +00:00
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
2023-06-21 18:29:07 +00:00
"<matplotlib.legend.Legend at 0x7f51b25237f0>"
2022-04-03 20:18:44 +00:00
]
},
2023-06-21 18:29:07 +00:00
"execution_count": 25,
2022-04-03 20:18:44 +00:00
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
2023-06-21 18:29:07 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEGCAYAAABYV4NmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAp6ElEQVR4nO3deZgcZbXH8e/p7lky2TcCJIFECLusEYIsRhAIa0QQgiwR0OgVBBX0Il5EVBSXi4AKGgUB5SZg2AKEJSzDToCwQ4KEhGwEErJPZu3uc/+omqQzzEwmmampnu7f53n66aq3tvN2Teqk3rcWc3dEREQ6WiLuAEREpDApwYiISCSUYEREJBJKMCIiEgklGBERiUQq7gDyxYABA3zYsGEbla1bt47u3bvHE1CeKPbfQPVX/VX/1us/c+bMT9x9YHPTlGBCw4YN4+WXX96orLKyktGjR8cTUJ4o9t9A9Vf9Vf/Rrc5jZvNbmqYmMhERiYQSjIiIREIJRkREIqE+mFaYGfPmzaO2tjbuUFpVXl7OkCFDKCkpiTsUEZH1lGBa0b17d3r27MmwYcMws7jDaZa7s3z5chYtWsTw4cPjDkdEZD01kbUimUzSv3//vE0uEJxl9e/fP+/PskSk+CjBbEI+J5dGXSFGESk+SjAiIsVq7cfw8E+gekUkq1eC6QLuuecezIzZs2fHHYqIFJInr4IZf4GalZGsXgmmC5g0aRIHH3wwkyZNijsUESkUy/4DM2+BkedA/x0i2YQSTJ6rqqrimWee4cYbb2Ty5MlxhyMiheKxK6CkAg79UWSb0GXKbXTFfW/zzodrOnSdu23bi8uP373Vee69917GjBnDTjvtRP/+/Zk5cyb77bdfh8YhIkVm/vMw+3447H+gR7PPqewQOoPJc5MmTWLcuHEAjBs3Ts1kItI+7jD9Mui5DYw6L9JN6QymjTZ1phGFFStW8Pjjj/Pmm29iZmQyGcyM3/3ud7o0WUS2zKypsOglOOGPUFoR6aZ0BpPHpkyZwplnnsn8+fP54IMPWLhwIcOHD+fpp5+OOzQR6YoyDfDoFTBwV9jra5FvTgkmj02aNIkTTzxxo7KTTjpJzWQismVm3gwr3ocjroBk9A1YaiLLY0888cSnyi644IIYIhGRLq9uLVReBcMOgRFHdsomdQYjIlIMnr0Oqj8Jzl46qQ9XCUZEpNCtWQLP/wn2OAkGd95tDkowIiKFrvLXQQf/YZd16maVYERECtnS2fDqP+Fz34B+nfvOKCUYEZFC9ujPoLQHHPrDTt+0EoyISKH64Bn4z4Nw8Pehe/9O37wSTJ7r0aNH3CGISFfkDo9cBr0Gw6j/iiWESBOMmX3fzN42s7fMbJKZlZvZcDObYWZzzOx2MysN5y0Lx+eE04flrOfHYfm7ZnZUTvmYsGyOmV2SU97sNkREisbbd8OHr8AXfwIl3WIJIbIEY2aDgQuAke6+B5AExgG/Af7g7jsCK4Fzw0XOBVaG5X8I58PMdguX2x0YA1xvZkkzSwJ/Bo4GdgNOC+ellW2IiBS+dH3wOP6tdoe9xsUWRtR38qeAbmbWAFQAS4DDgMaH4NwC/Ay4ARgbDgNMAf5kwRMdxwKT3b0OmGdmc4D9w/nmuPtcADObDIw1s1mtbGPLPXgJfPRmu1bxKVt/Fo6+qmPXKSLy8k2w8gM4/U5IJGMLI7IE4+6Lzez3wAKgBngEmAmscvd0ONsiYHA4PBhYGC6bNrPVQP+w/IWcVecus7BJ+QHhMi1tYyNmNgGYADBo0CAqKys3mt6rVy/Wrl0LQFlDPYlMuukq2iXbUE9duP7WrG3DPLW1tZ+KvyNUVVVFst6uQvVX/bta/ZPpdYx64ZdU9dmT1xclYXHlFq+rvfWPLMGYWV+Cs4/hwCrg3wRNXHnD3ScCEwFGjhzpo0eP3mj6q6++Ss+ePYORE66OJIa2dA6tj6EV5eXl7LPPPu0PqInKykqa/i7FRPVX/btc/R+9AtJr6XvKHxm97d7tWlV76x9lJ/+XgHnuvszdG4C7gIOAPmbWmNiGAIvD4cXAUIBwem9geW55k2VaKl/eyjZERArX6sXwwvXw2VOgncmlI0SZYBYAo8ysIuxLORx4B3gCODmcZzxwbzg8NRwnnP64u3tYPi68ymw4MAJ4EXgJGBFeMVZKcCHA1HCZlrbR5VRXVzNkyJD1n6uvjuZMSkQKwBO/As8Gr0LOA1H2wcwwsynAK0AaeJWgOeoBYLKZ/TIsuzFc5Ebgn2En/gqChIG7v21mdxAkpzRwnrtnAMzsfOBhgivUbnL3t8N1/XcL2+hystls3CGISFfw8dvw2m1w4HnQd/u4owEivorM3S8HLm9SPJcNV4HlzlsLfLWF9VwJXNlM+TRgWjPlzW5DRKRgTb8cynvBIRfFHcl6upNfRKSrm1sJc6bDIRdDRb+4o1lPCWYTgi6d/NYVYhSRiGSzMP2n0Hso7D8h7mg2ogTTikwmw/Lly/P6AO7uLF++nPLy8rhDEZE4vHUnLHk9eNdLSX4dB6K+k79LW7duHWvXrmXZsmVxh9Kq8vJyhgwZEncYItLZ0nXw+M+Dp4J8ttku7FgpwbTC3Rk+vHNf0CMi0mYvXA+rFsCZd0Mi/xqk8i8iERHZtFUL4Mnfwi7HwQ6HxR1Ns5RgRES6ood+HHyP+XW8cbRCCUZEpKt59yGYfT984UfQZ7u4o2mREoyISFdSXw0P/hAG7Ayjzos7mlapk19EpCt5+n+D/pfx90Mqv1/WqzMYEZGu4pP34NlrYc9xMPyQuKPZJCUYEZGuwB0euAhKKuDIX8QdTZuoiUxEpCt4606Y9yQc83vosVXc0bSJzmBERPJd7Wp4+FLYZm8YeU7c0bSZzmBERPLdE7+CqqVw2mRIJOOOps10BiMiks+WvA4vToTPnQuD9407ms2iBCMikq+yWbj/B1DRP29eg7w51EQmIpKvXr0VFr8MJ/4VuvWNO5rNpjMYEZF8tO6T4DXI2x8Ee54adzRbRAlGRCQfTb8c6qvg2P8Fs7ij2SJKMCIi+Wb+8/Dav+DA82GrXeOOZospwYiI5JNMAzzwA+g1JHhachemTn4RkXwy46+w9B049TYo7R53NO2iMxgRkXyxejFU/hpGHAW7HBt3NO2mBCMiki8e/jFk03D0b7psx34uJRgRkXzw3qPwzr1wyMXQb3jc0XQIJRgRkbg11MK0i6H/jnDQBXFH02HUyS8iErdnr4GV8+DMeyBVFnc0HUZnMCIicVr+Pjx9NexxEuzwxbij6VBKMCIicXGHaT+EZCkceWXc0XQ4JRgRkbi8dSe8/xgc9hPotU3c0XQ4JRgRkTis/Sjo2N92X/jcN+OOJhJKMCIinc0dpn4XGmqCR/EnC/N6q8KslYhIPnvlFnjvERjzGxi4U9zRREZnMCIinWnFPHj4JzDsENh/QtzRREoJRkSks2QzcM93wBLw5RsgUdiH4EhrZ2Z9zGyKmc02s1lmdqCZ9TOz6Wb2XvjdN5zXzOw6M5tjZm+Y2b456xkfzv+emY3PKd/PzN4Ml7nOLHh4T0vbEBGJ1QvXw4LnYMxV0Gdo3NFELur0eS3wkLvvAuwFzAIuAR5z9xHAY+E4wNHAiPAzAbgBgmQBXA4cAOwPXJ6TMG4Avpmz3JiwvKVtiIjEY+kseOznsPOxsPfX4o6mU0SWYMysN3AocCOAu9e7+ypgLHBLONstwJfD4bHArR54AehjZtsARwHT3X2Fu68EpgNjwmm93P0Fd3fg1ibram4bIiKdL10Pd02Asp5w/LUF8aTktojyKrLhwDLgH2a2FzATuBAY5O5Lwnk+AgaFw4OBhTnLLwrLWitf1Ew5rWxjI2Y2geBsiUGDBlFZWbnR9Kqqqk+VFZti/w1Uf9W/I+o/bN5tDPvoDd7a/RI+efnt9gfWSdpb/ygTTArYF/iuu88ws2tp0lTl7m5mHmEMrW7D3ScCEwFGjhzpo0eP3mh6ZWUlTcuKTbH/Bqq/6t/u+i+aCU/eCXuOY4+v/LhD4uos7a1/lH0wi4BF7j4jHJ9CkHA+Dpu3CL+XhtMXA7m9XkPCstbKhzR
2022-04-03 20:18:44 +00:00
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(range_n, zeros_A, label=\"A\")\n",
"plt.plot(range_n, zeros_L, label=\"L\")\n",
"plt.xlabel(\"n\")\n",
2023-06-21 18:29:07 +00:00
"plt.ylabel(\"#!zéros\")\n",
2022-04-03 20:18:44 +00:00
"plt.grid()\n",
"plt.legend()"
2022-03-24 07:10:59 +00:00
]
2022-04-03 20:42:26 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2023-06-21 18:29:07 +00:00
"On observe que la matrice L possède bien plus de valeurs différentes de zeros que la matrice A, il n'est donc pas bénéfique de l'utiliser pour gagner en espace mémoire (surtout pour n très grand). Le plus efficace serait de stocker A via une matrice creuse."
2022-04-03 20:42:26 +00:00
]
2022-03-15 07:04:08 +00:00
}
],
"metadata": {
"kernelspec": {
2023-06-21 18:29:07 +00:00
"display_name": "Python 3.10.8 64-bit",
2022-03-15 07:04:08 +00:00
"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",
2023-06-21 18:29:07 +00:00
"version": "3.10.8"
},
"vscode": {
"interpreter": {
"hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
}
2022-03-15 07:04:08 +00:00
}
},
"nbformat": 4,
"nbformat_minor": 4
}
