2022-03-15 07:04:08 +00:00
{
"cells": [
{
"cell_type": "code",
2022-04-03 20:18:44 +00:00
"execution_count": 54,
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",
2022-04-03 20:18:44 +00:00
"execution_count": 55,
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",
" return coordinates, elements3, dirichlet, neumann\n"
]
},
{
"cell_type": "code",
2022-04-03 20:18:44 +00:00
"execution_count": 56,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"def show(coordinates, u) -> None:\n",
" \"\"\"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-03-15 07:04:08 +00:00
"\n",
" Returns:\n",
" None, plots a figure\n",
" \"\"\"\n",
"\n",
" ax = plt.figure().add_subplot(projection=\"3d\")\n",
" ax.plot_trisurf(\n",
" coordinates[:, 0], coordinates[:, 1], u, linewidth=0.2, antialiased=True\n",
" )\n",
" 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",
2022-04-03 20:18:44 +00:00
"execution_count": 57,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"def f(x, y):\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",
"def u_ex(x, y):\n",
" return np.sin(np.pi * x) * np.sin(np.pi * y)\n",
"\n",
"\n",
"def u_d(x, y):\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",
2022-04-03 20:18:44 +00:00
"execution_count": 58,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [
{
2022-04-03 20:18:44 +00:00
"name": "stdout",
"output_type": "stream",
"text": [
"coords = [[0. 0. ]\n",
" [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",
"elems3 = [[0 1 3]\n",
" [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",
"dirichlet = [[0]\n",
" [1]\n",
" [2]\n",
" [5]\n",
" [8]\n",
" [7]\n",
" [6]\n",
" [3]]\n",
"\n",
"neumman = []\n"
]
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",
" f\"coords = {coords}\",\n",
" f\"elems3 = {elems3}\",\n",
" f\"dirichlet = {dirichlet}\",\n",
" f\"neumman = {neumann}\",\n",
" sep=\"\\n\\n\"\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[M^A_T]_{ij} = \\displaystyle \\int_T \\nabla \\eta_i (x, y) \\eta_j (x, y) \\ dx \\ dy\n",
"$$\n"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2022-04-03 20:18:44 +00:00
"execution_count": 59,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
2022-04-03 20:18:44 +00:00
"def calcul_alpha(x, y):\n",
" \"\"\"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",
" float: le coefficient alpha.\n",
" \"\"\"\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)"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1. , -0.5, -0.5],\n",
" [-0.5, 0.5, 0. ],\n",
" [-0.5, 0. , 0.5]])"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def raideur(triangle):\n",
" \"\"\"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",
2022-04-03 20:18:44 +00:00
"execution_count": 61,
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:18:44 +00:00
"def assemblage(coordinates, elements3):\n",
" \"\"\"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",
"execution_count": 62,
"metadata": {},
"outputs": [],
"source": [
"def second_membre(coordinates, elements3):\n",
" \"\"\"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",
"execution_count": 63,
"metadata": {},
"outputs": [],
"source": [
"def calcul_Ud(coords, dirichlet):\n",
" \"\"\"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",
"execution_count": 64,
"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",
"execution_count": 65,
"metadata": {},
"outputs": [],
"source": [
"def untildage(x, dirichlet, U_d):\n",
" \"\"\"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",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAADyCAYAAACRQVPgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAADDpElEQVR4nOz9eZRdWVreCf/2PuMdYx6kCIWk1JRSKmdJFLC+xssYD7hNN+tze+iv23hhu2k34Ak3YEyzABswZrINNhgKcBnbYBowFBjKmIIqY6gsKVNSKjWPISkiFIrxxp3OtIfvjxM3CEVKiiEl5VDxrJUrpdA959w4Zz/nffc7PK+w1rKNbWzjowf5fn+BbWxjG08H2+TexjY+otgm9za28RHFNrm3sY2PKLbJvY1tfESxTe5tbOMjCnedf9/Ok21jG08f4mmcdNtyb2MbH1Fsk3sb2/iIYpvc29jGRxTb5N7GNj6i2Cb3NrbxEcU2ubexjY8otsm9jW18RLFN7m1s4yOKbXJvYxsfUWyTexvb+Ihim9zb2MZHFNvk3sY2PqLYJvc2tvERxTa5t7GNjyi2yb2NbXxEsU3u9wHWWtI0RSnFtrT0Np4W1hNr2MYThjGGNE2J43jlZ47j4HkeruviOA5CPJXe/W18gUGsYzm2zcoTgrUWpRRKKYQQZFm28nNrLcaYFVInSUKlUsH3/W2yf2HgqTzgbcv9DNBxw1cTuAMhBEIIpJQrn71x4wZ79uyhWCwC25Z9G1vDNrmfMpRSTExMoLVmZGQEIcSKtX4YSTtkdxwHx3FWrHoURSufd1135b9tsm/jUdgm91PCajfcGLPijm8WD7PsWmuUUiufcV13xbJLKbfJvg1gm9xPBcYYsixbccM71nqjeNznO+frYC3ZhRAPWPZtsn/hYpvcTxAdonWCZR1r+yiyPso13wweRnal1Mp36Fh913XxfX+b7F9A2Cb3E4K1lizL0Fq/i3Bryb2eNd+spV977FqyT0xMALBjx45ty/4FhG1yPwF0ctcdS/ywiPj7Vayy+vt0AnRZlj1g2Tt7dsdxtsn+EcI2ud8D1uauO274WjzJPfd7RScS38HDyN4Jzrmu+9CX1TY+HNgm9xaxNnf9OAJshaxPmtyP+n4PI3uapiRJAuRxA8/zViz7Ntk/PNgm9xbQCZo9yg1fi61Y7vcLjyP76uDcajd+Gx9MbD+ZTaDjwp49e5YkSTa8P32/3fL3cq4O2TvBN4A0TTl16hS1Wo16vU673V5J/W3jg4Nty71BrM5dd4JnG8X7GVB7klgdmIvjeMXCp2lKmqYA25b9A4Rtcq+Dtbnrjmv6NMn9YXkZrC6VhT9qgllL9tV18dtkf3bYJvdj8KjctRBiUy7oh4Wsm8WjmmA66JA9SZJ3Bei2yf70sU3uR+BxueunbYk/Ki+D9churX3Ahe+k3rbxZLBN7jVY7YY/Knctpdy23FvAw8hujFkRrpiammJ0dBTf97c73p4Atsm9ChvNXb9XS2yMYWJiAs/z6Onpwff993T+DyvW3uPp6Wl27ty5rVLzhLBN7mWsV0K6Gu8loBZFEefOnaO7u5s0TZmcnMQYQ1dXFz09PXR3d7/XX+VDjbU59rW97Ntk3zi+4Mm90RLS1dhqQG12dparV69y5MgRyuXyioegtaZWq7G4uMj4+PjKnlQIQVdX1xds0OlhvezbZN84vqDJbYxhZmYGay09PT0bXhhbcZvr9TpKKY4fP47v+yupNcgXaF9fH319fQBcvXoVx3GYmZnh+vXrK+57T08PlUply6IPH3ZshOzbKjV/hC9Icq8OmjUaDay19Pb2bvj4zQTU4jjm3LlzCCF4/fXXN7TYXNelWq3S39+/co7FxUUmJiZoNpsUCoUVsheLxXXP+VHdvz9KpaYTDO3A932CIPiC63j7giP3WjfccZwHrOhGsFHL3XHD9+3bx+Tk5JYXVhiG7Nixgx07dmCtpd1us7i4yM2bN4miiHK5TE9PD729vQRBsKVrfBTwMLLfuXMHz/MYHBx8oL31C6GX/QuK3A+TP9psWgvWD6gZY7h+/Tr1ep1jx44BrAgmbATrCTmUSiVKpRKjo6NYa2k0GiwuLnLp0iWyLKNara5Yds/zNvW7fZTQuY8dN/0LTZLqC4Lcj8tdbzbyDY8PqHXc8N7e3hU3/GnWogshqFarVKtVdu/ejTGGpaUlFhcXuXv3LtZaHMehXC6jtX4gGv2FAGPMA3JX60lSfZTI/pEn9+Pkj2Dzke/OMQ8j39zcHFeuXOH5559fCY497vNPA1LKFasNubTy9evXabVanD59GsdxVlz4SqXykY/Eryb3WjyM7B8llZqPNLk3krveqlu++hhrLdevX6dWq3Hs2LF37Xvfz/JT13UplUp0dXWxY8cO0jRlcXGRe/fuceXKFYIgWHkZlMvl923xPq2X3+PIvRabUan5MJD9I0nuzeSut0Lu1eRLkmSlKOXYsWMPfdgfpIoz3/cZGhpiaGgIyItqFhcXuXPnDs1mk1KptEL2QqHwzBbvk1CCfRg2Q+61eJxwxdTUFENDQxSLxQ+sJNVHjtybkT+CrZNba838/DyXL1/m0KFDK2mrR31+Lbkfd91n+TIoFAoUCgV27tyJtZZWq8Xi4iLXr18njmMqlcoK2Z9mJP6DSO61WE32hYUFhoaGHlCp6Vj2D0ov+0eK3B0XqnOjN7JYtkru2dlZZmZmeP311wnDcN3PryXrB8WSr4YQgnK5TLlcZteuXRhjViLxFy9eRCm1Uib7pL//kyThszhvJzi5upcd/ki44hu+4Rv49m//dp5//vknfu2N4iNB7o4bfvnyZbq7uxkcHNzwsZsld5IkjI+P43kex44deyoL52lY7q1YRSklXV1ddHV1sWfPHrTWK5H4drvNW2+9RXd3Nz09PXR1db2nSPyHwXKvPe/q33d1OSzklr1QKDzx624GH3pyr85dP42c9WosLCxw6dKllYKIzQRqPgpwHIfe3l56enpYWFjgpZdeolarMTc3x40bN3Bd94FI/GZ+b2vth8pyw+Ofa6vVolwuP5XrbhQfWnI/bHSP4zhbcrHXO8Zay82bN5mfn+f111+nXq+ztLS05e+++rzPIgD3pL2Azvf2PI+BgQEGBgaA3KtZXFxkcnKSRqNBGIYr+/VSqfRYMjxsvPGTwNMi93r3tN1ub5N7K3hU7vpJpLXWIk1Tzp07R6VSWXHDG43GttLnQxAEAcPDwwwPD2OtXYnEj4+Pr1iy1ZH41fiwueXrIcuyd/XpP2t86Mj9uNy14zhorTd1vseRuxNIOnjw4Ip16hzzXqyhtZbJyUmyLKOvr+9dzR8fpNTZw7ARIgohKBaLFItFRkZGsNbSbDZZXFzk6tWrJEnyQJns07SwTysu8kHHh4bcG8ldSyk33QTyMHJba7l16xazs7O89tpr77Iy74V8SinOnz+P67oUi8UHmj96e3s31Z32fmErVlYIQaVSoVKpMDY2hjGGer2+4sZ3ntvc3Bzd3d247gd3aa737J+WF7JZfHDv4CpsNHf9JNzyNE155513KJVKHD9+/IloqHXQbDY5d+4ce/bsYXBwEKXUA80fCwsLnD9/njiOKRQKhGFId3f3B64e/El4FVJKuru76e7uZu/evSuFNLVajfHxcYQQK1b9gyZYsZEa/Q8CwT/w5N7M6J6tBtQ66LjhBw4ceGw6bSuWO8syzp07x4svvkilUlnpTOqcr9P8sWfPHiYnJ2k2mywsLHDz5s33FIV+WnjS30EIQaFQYP/+/UB+vxYXF98lWNHb2/u+lsnC+uRWSn0gXsgfWHJvRf5ISrnpPXfnWrdu3WJmZuahbvjDrrPRl4gxhitXrpBlGV/6pV+64m4+bnFKKSkWi+zatQvIo9ALCwtMTEzQaDRWSkR7e3vfl1zq07BKa8/Z6cHuvGQ7ghV37959QLCicw+eJdm11o9dj61Wi1Kp9My+z6PwgSS31ppms4nv+5uq192Ku5ymKVEUEcfxI93wtdio5Y7jmLfffpvBwUEKhcKG95Frzx8EwQNiDZ0S0U5gqlM19qz6t58Gudc
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAADyCAYAAACRQVPgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAADENElEQVR4nOz9d5heWX7Xi37WWju9sXJQlVTKsdXqoFbPjDn3mnuI9nNteLgc7oXLwVwbzsUH+8A55gAGXxwwDjgQbGODx8YmGYPHOHucxvYAnh6pW2rlHEqqKlUOb9pxrXX/2PXWlKolVWhJ091T3+fR0+rSu/d+a+/13b/f+oXvT1hr2cY2tvHRg/xif4FtbGMbzwfb5N7GNj6i2Cb3NrbxEcU2ubexjY8otsm9jW18RLFN7m1s4yMKZ51/386TbWMbzx/ieZx023JvYxsfUWyTexvb+Ihim9zb2MZHFNvk3sY2PqLYJvc2tvERxTa5t7GNjyi2yb2NbXxEsU3ubWzjI4ptcm9jGx9RbJN7G9v4iGKb3NvYxkcU2+TexjY+otgm9za28RHFNrm3sY2PKLbJvY1tfESxTe4vAqy1JElClmVsS0tv43lhPbGGbTxjGGNIkoQoilZ+ppTCdV0cx0EphRDPpXd/G19iEOtYjm2z8oxgrSXLMrIsQwhBmqYrP7fWYoxZIXUcx1QqFTzP2yb7lwaeywPettwvAG03fDWB2xBCIIRASrny2du3b7Nnzx6KxSKwbdm3sTVsk/s5I8syxsbG0FozPDyMEGLFWj+OpG2yK6VQSq1Y9TAMVz7vOM7Kn22yb+NJ2Cb3c8JqN9wYs+KObxaPs+xaa7IsW/mM4zgrll1KuU32bQDb5H4uMMaQpumKG9621hvF0z7fPl8ba8kuhHjEsm+T/UsX2+R+hmgTrR0sa1vbJ5H1Sa75ZvA4smdZtvId2lbfcRw8z9sm+5cQtsn9jGCtJU1TtNbvIdxacq9nzTdr6dceu5bsY2NjAOzYsWPbsn8JYZvczwDt3HXbEj8uIv7FKlZZ/X3aAbo0TR+x7O09u1Jqm+wfIWyT+31gbe667YavxbPcc79ftCPxbTyO7O3gnOM4j31ZbePDgW1ybxFrc9dPI8BWyPqsyf2k7/c4sidJQhzHQB43cF13xbJvk/3Dg21ybwHtoNmT3PC12Irl/mLhaWRfHZxb7cZv44OJ7SezCbRd2HfffZc4jje8P/1iu+Xv51xtsreDbwBJknDmzBkWFxep1Wq0Wq2V1N82PjjYttwbxOrcdTt4tlF8MQNqzxKrA3NRFK1Y+CRJSJIEYNuyf4CwTe51sDZ33XZNnye5Pywvg9WlsvCFJpi1ZF9dF79N9heHbXI/BU/KXQshNuWCfljIulk8qQmmjTbZ4zh+T4Bum+zPH9vkfgKelrt+3pb4o/IyWI/s1tpHXPh26m0bzwbb5F6D1W74k3LXUspty70FPI7sxpgV4YqJiQl27tyJ53nbHW/PANvkXoWN5q7fryW21vLw4UMcx6GrqwvHcZ76+Y8q1t7jyclJhoaGtlVqnhG2yb2M9UpIV+P9BNTiOObChQuUSiUARkdHEULQ1dVFV1cXHR0d7+8X+ZBjbY59bS/7Ntk3ji95cm+0hHQ1thpQm5ub49q1axw+fJiOjo4VDyFNUxYXF5mZmeHWrVukaUqSJDiOQ6VS+ZJdwI/rZd8m+8bxJU1uYwzT09NYa+nq6trwwtiK21yv17l9+zYnT54kCIKV1BqA67r09fXR19cHwPXr11FKMTY2RqPRIAgCuru76erqolgsbln04cOOjZB9W6XmC/iSJPfqoFm9XsdaS3d394aP30xALY5jLl26hLWWN954Y0Oegeu6VKtVent7sdYShiELCwvcuXOHMAwplUorZA+CYN3zfVT3709SqWkHQ9vwPA/f97/kOt6+5Mi91g1XSj1iRTeCjVru+fl5rl69yt69e3n48OGWcrpCCIrFIsVikeHhYay1NBoNFhYWuHbtGkmS0NHRsbJnd11309f4qOBxZL9//z6u69Lf3/9Ie+uXQi/7lxS5Hyd/tNm0FqwfULPWcufOHebm5jh58iRCCCYmJjZ8/vWEHCqVCpVKhZGREYwx1Go15ufnGRsbwxhDZ2cnXV1ddHZ2PhKg+lJD+z623fQvNUmqLwlyPy13vdnINzw9oJYkCRcuXKBSqay44c+zFl1KSWdnJ52dnUCutrq0tMT8/Dx3795FSomUkkqlgjHmS64ibPXvvBFJqo8S2T/y5H6a/BFsPvLdPuZx5FtYWODKlSscOnRoJTj2tM8/DziOQ09PDz09PUD+srl16xZLS0u8/fbb+L6/4sKXy+UP9eLdCJ72Qnsc2T9KKjUfaXJvJHe9Vbd89THWWu7evcvMzAyvv/46hULhkc9/MctPPc+jUqnQ1dXFjh07iKKI+fl57t+/T6PRoFgsrgTnCoXCF23xPq+X32a8lc2o1HwYyP6RJPdmctdbIfdq8iVJwsWLFymXy5w6deqx1/ogVZwFQcDQ0BBDQ0NYa2m1WiwsLHDr1i2iKKJcLq+Q3ff9F/a9noUS7OPwfrYiTxOumJiYYGBggGKx+IGVpPrIkXsz8kewdXJrrVfc8IMHD9Lf3//Uz68l99Ou+6JeBkIISqUSpVKJnTt3Yq2lXq8zPz/PlStXyLJsJRLf2dn5XCPxH0Ryr8Vqss/PzzMwMPCISk3bsn9Qetk/UuRuu1DtG72RxbJVcs/NzTE1NfVYN/xxn19L1g+KJV8NIQTVapVqtcqePXvQWq9E4u/fv79S7NPV1fXMv//zCvY9r/NqrVdGPgGPeHJJkvAN3/ANfMu3fAtHjhx55tfeKD4S5G674deuXaOzs/OpVnQtNkvuNE0ZHR1FSsmbb775XBbO87DcW7GKSqkVMkMeiV9cXGR2dpZWq8W5c+dW/r1Sqbyve/FhsNxrz7vaZV9dDgu5ZV/vpf+88aEn9+rc9fPIWa/G4uIily9fpre3d1Nu1wdpH/Z+4DgOvb299PT0sLS0xEsvvcT8/DwTExPU63WCIFghe6lU2tTvba39UFluePpzbTablMvl53LdjeJDS+7Hje5RSm3JxV7vGGsto6OjTE5O8tprr9FoNFhaWtryd1993sctkGdtuZ+1F9D+3p7nMTg4yODgIABhGDI/P8+9e/dWFneb7OtZsceNN34WeF7kXu+etlqtbXJvBU/KXT+LtNZapGnKpUuXCIJgxQ1vNpvbSp+PQaFQYHh4eKVMttlsMj8/z40bN4jjmEqlshKJ9zzvkWM/bG75ekjT9D2/44vGh47cT8tdK6XQWm/qfE8j99LSEpcvX2bfvn0r1ql9zPuxhm2xhjRN6enped958ReNjRBRCEG5XKZcLq+UybYj8ePj42itHymTfZ4W9nnFRT7o+NCQeyO5aynlpptAHkfudsPBw4cPefXVVykWi4/8+/shn9aay5cvA1AsFlesWrVaXbFqH3RsxcpKKeno6KCjo4O9e/eitWZxcZGFhQXu3bu38gwWFhbo6Oj4oqeRnob1nv3z8kI2iw8FuTeau34Wbnmaply+fBnP8zh16tRjGy+2ch3IgywXLlxg165dDA4OkmXZe5o/Hjx4QJIkFAoFisXiB3KhPwuvQin1SJlsO902PT3NrVu3cBxn5WX3QROsaKfBnoYPAsE/8OTezOierQbU2mi74Xv37mXHjh1PPWazCzxNU86fP8/x48epVqsrnUnw3uaPtkjD9PQ0N2/eXKkH7+7u3nQU+nnhWX8HKSXFYpFDhw4BeR98u9Ot0WhQKBRWgnNbFax4VliP3FmWfSC68T6w5N6K/JGUctN77va17t+/z/j4OK+88sqKvtnTrrPRl4gxhps3b5KmKV/2ZV+2UuW13kuqVCqxa9cugBWxhtVR6O7ubrq7u19oiWgbz8MqrT2n7/vs2LGDHTt2rAhWzM/PrwhWrI7Eb0Sw4llCa/3U9dhsNtddQy8CH0hya61pNBp4nrepet2tuMtZlhGGIfV6nTfffHNDb9yNWu62GGJ3dzeFQmHD5Ztrz18oFCgUCiv14I1G45ES0dW
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"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-03-15 07:04:08 +00:00
"show(coords, x_untild)\n",
2022-04-03 20:18:44 +00:00
"\n",
"# on compare avec le résultat théorique exacte\n",
2022-03-15 07:04:08 +00:00
"show(coords, u_ex(coords[:, 0], coords[:, 1]))"
]
},
{
"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",
2022-04-03 20:18:44 +00:00
"execution_count": 67,
"metadata": {},
"outputs": [],
"source": [
"def f(x, y):\n",
" return 1\n",
"\n",
"\n",
"def u_d(x, y):\n",
" return 1\n",
"\n",
"\n",
"def g(x):\n",
" return 1"
]
},
{
"cell_type": "code",
"execution_count": 68,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"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",
")"
]
},
{
"cell_type": "code",
2022-04-03 20:18:44 +00:00
"execution_count": 69,
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]])"
]
},
2022-04-03 20:18:44 +00:00
"execution_count": 69,
2022-03-16 14:44:07 +00:00
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def raideur_quadrangle(quadrangle):\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",
2022-04-03 20:18:44 +00:00
"execution_count": 70,
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:18:44 +00:00
"def assemblage_quadrangle(coordinates, elements4):\n",
" \"\"\"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",
"execution_count": 71,
"metadata": {},
"outputs": [],
"source": [
"def second_membre_quadrangle(coordinates, elements4):\n",
" \"\"\"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",
"execution_count": 72,
"metadata": {},
"outputs": [],
"source": [
"def condition_neumann(coordinates, neumann):\n",
" \"\"\"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",
"execution_count": 73,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAADyCAYAAABkv9hQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABfUUlEQVR4nO29d3xb9b3//zySLHnvbccrjhNnOHESJwFKwiiQBkhSWiCllz0KJdAWektp++ul3G8pty23ixZuIW1pyyYFUkYoexWyvEdsx3svecnW1vn94UiRZY0jWbKdWM/HI5DYR+cc2Xqdz/vznoIoigQJEuTMRjbfNxAkSJDAExR6kCCLgKDQgwRZBASFHiTIIiAo9CBBFgFBoQcJsghQePh+MPYWJEjgEQJ9geCKHiTIIiAo9CBBFgFBoQcJsggICj1IkEVAUOhBgiwCgkIPEmQREBR6kCCLgKDQgwRZBASFHiTIIiAo9CBBFgFBoQcJsggICj1IkEVAUOhBgiwCgkIPEmQREBR6kCCLgKDQ5wFRFDEYDJhMJoLttoPMBZ4aTwTxMxaLBYPBgE6ns31NLpcTEhKCQqFALpcjCAHvQxBkkSF4WFGCy42fEEURk8mEyWRCEASMRqPt66IoYrFYbALX6/VERUWhVCqDwl8cBPwXHFzR5wCrqW4vZiuCICAIAjKZzHZsU1MTOTk5hIeHA8EVP8jsCQo9wJhMJjo7OzGbzWRkZCAIgm0VdyZYq/Dlcjlyudy22mu1WtvxCoXC9ico/CBSCAo9QNib6haLxWaye4uzFd9sNmMymWzHKBQK24ovk8mCwg8yg6DQA4DFYsFoNNpMdesqLhV3x1vPZ8VR+IIgTFvxg8IPAkGh+xWr6KyONusq7Eq4rsx3b3AmfJPJZLsHqzWgUChQKpVB4S9SgkL3E6IoYjQaMZvNM8TnKHRPq7y3FoDjax2F39nZCUBaWlpwxV+kBIXuB6yxcesK7cyzPl+JMfb3Y3XuGY3GaSu+dY8vl8uDwj9DCQp9FjjGxq2muiP+3KPPFqtH34oz4VsdewqFwumDK8jpR1DoPuIYG3cnBl+E62+hu7o/Z8I3GAzo9Xpgys8QEhJiW/GDwj89CQrdB6wON1emuiO+rOjzhTvh2zv27E39IAuf4G/JC6xmbnl5OXq9XvJ+dr5N99mcyyp8q+MOwGAwcOTIEUZGRhgbG2NyctIWTgyyMAmu6BKxj41bHW9SmU9nnD+xd+rpdDrbym8wGDAYDADBFX+BEhS6Bxxj41bzNZBCP10eDPbpunCqQMdR+PZ5+kHhzw9BobvBVWxcEASvzNTTRbje4qpAx4pV+Hq9foZzLyj8uSUodBe4i40HeoU+Ux4MnoQviuI0M98azgvif4JCd8DeVHcVG5fJZMEV3QecCd9isdiacHR3d5OZmYlSqQxW5vmZoNDtkBob98cKHagU2NMJx59xb28v6enpwe47ASAo9JN4SmO1J9DOuMWMYwzfsRY/KHzfWPRCl5rGak+gnXHBB8MUzmrxg8L3jUUtdIvFQn9/P6IoEhcXJ/lDcjoK90wQgBThB7vvOGdRCt3e4TY+Po4oisTHx0t+/WydcSMjI1RXVyOTyYiLiyMuLo6YmBib2bqQMuMWMq6671gdqVaUSiUqlWpRV+YtOqE7mupyudyWDCMVX1d0URRpa2ujr6+PoqIi5HI5o6OjDA4O0tTUhEKhIC4uzhZ6CuIdzoTf3t5OSEgIycnJ00pyF1st/qISurMWT96uzuCbM85sNlNeXk5oaCglJSWYzWYsFguJiYkkJiYCU22eh4eH6e/vZ3h4mJ6eHuLj44mLiyM8PHzRfCj9hfUBazXlF3PbrUUhdHexcW9FC9474zQaDaOjo6xZs4bU1FRgqgLOEZVKRWpqKjqdjrCwMKKiohgeHqa5uRmtVktkZCRxcXHEx8ejUqm8uufFisVimdbSy1PbrTNV+Ge80N21eALvRWt9jZSHgyiKdHR00NnZSWRkpE3kUq8RHh5OeHg4GRkZiKKIRqNBrVZTV1eH0WgkOjqa+Ph4YmNjCQkJ8eo9LBbshe6IM+Gfqd13zmihS4mN+2q6e3qNyWSipqYGuVzOhg0bKC8vl3x+Vwk2UVFRREVFkZ2djcViYXR0lOHhYdrb222RA0fH3ulCoHwS7oTuiDfdd0434Z+RQvcmNu6L0D2t6OPj41RVVZGTk0N6ejpms9nvH2R7jz1MPVhGRkZmOPbi4+NPC8eePzriOsMboTvirglHd3c3KSkphIeHnxZtt844oXvT4gl8F7qzPTZAZ2cnHR0dFBUVERkZaTveUWzurutLeE2hUDh17HV1dTE0NIRSqcRsNi9Yx95CFLoj9sJXq9WkpKRM675jXfEXYi3+GSV0q5ll/aFL+eD4y+tuNpupra1FFEVKSkpQKE79aJ0JN9CrrNWxl5qaajPtBUGgubmZyclJoqKiFpRjz5+CnIvzms1m29gsOPX7tNbi7927lx/96EesWLHC79f2hTNC6FZT/fjx48TGxpKcnCz5tf4w3TUaDVVVVSxZssQ2X202+DthRhAElEolqampC9axdzqs6I7ntTfr7VNyYWrFDwsL8/t1feW0F7p9bHwuYuLW11iv093dTWtrK6tXryY6Otrp8QvNTJ6tYy8Q1ogoiqfVig7uf68TExO2rdtC4LQVurPxR3K53KfV2dc9ek1NDUajkU2bNk0z1aXibqLqXKbAeuPYi4qKCsjq62yktL/OGwihe/qZTk5OBoU+W1zFxgMVKnPEYDDQ3d1NXl4eWVlZC27Fni3uHHvj4+OoVCoMBgMTExN+c+ydbqa7J4xGI0qlcs6v64rTTujuYuNyudylN9wV3gq9r6+P5uZm4uPjyc7O9upaUlkI1W722Dv2RFFkYmKCysrKGRl7cXFxhIaG+nSNQK68gTjv6fZwP22ELiU2LpPJvC5QkSp0i8VCfX09Wq2WFStWMDQ05NV1zhQEQSA0NJTQ0FDWrFkzzbF3/Phxnx17gVrRA4Gnh/BCfC+nhdClxsYDZbprtVoqKytJTk5mxYoVjI6OBnRYgS9+g7nE/oPur4y9QK28gcAaWnPHQhP7ghe6N+OPfHXGuWNgYICGhgZWrlxpc1YF2rReSB8QV7h72Hrj2LOeJ1DOuEDgSegmk2nBpSAvWKH70uJJJpN5vUd3hcVi4cSJE4yNjVFSUjLNseKL5eCIpw/1QtqjO+LNauXJsRcaGkp8fPyCE4Y7zGaz28/jxMQEERERc3hHnlmQQjebzWg0GpRKpVf5w/4QIIBOp6OyspKEhAQ2bNjgtOIt0Cv6mSJ0Rxwde1qtluHhYbq7u5mcnMRkMs3asRdoHJNlHFloMXRYYEK3xsb1ej1lZWVs2bLFqw+UP4Q+NDTE8ePHWbFiBQkJCU6PWehCnAv8YWbbl+KqVCpGR0dJTk62OfYMBgMxMTE24XubsReo35En0z24orvBPjZuzVbz9sPkS3jN/vpNTU0MDw+zceNGt/nf/rIcXBGIB4k/97+BzIxz5djr6OhAFEViY2OJj4+XVIob6Dx3V2g0muCK7gxveqq7w1cBWiwWjh07RkxMDBs3bvR4/dkK0dqQYmJigoSEhIDXj/tbmIHKjHMUpVTHXlxcHFFRUTNeH0ihuzuvRqMJruj2SBl/5A2+CF2tVjM5Ocny5ctJSkoK2HWsmEwmqqurCQkJISYmhsHBQU6cOIFSqSQ+Pt7WjXYhbw0CIXQp53Tl2Ovu7rY59qwe/fDw8Hlb0Rda+ivMo9C9rRuXgjemuyiKtLS0MDAwQHh4uO3DIwVfV3SNRkNlZSU5OTmkpKRgNBptlXY6nQ61Wk1rayujo6OEhITYPNILKZUSAvMQ8kWUrhx71lLcsLAwDAYDOp3Or449Kc644IqO/0x1R6SutAaDgaqqKiIiIigpKeHQoUNerVK+CN1oNFJRUcGaNWuIjo6e8UAKDQ0lPT2d9PR0+vv7UavV6HQ6qqurbQ0jpO5N54L5WNHd4azH3sDAAK2trX5x7NljNpvdFjEteq+7t7Fxb3/5UoQ+MjJCTU0N+fn5pKSkTHu
"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",
"show(ccs, x_untild)"
]
},
{
"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",
"execution_count": 74,
"metadata": {},
"outputs": [],
"source": [
"def f(x, y):\n",
" return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\n",
"\n",
"\n",
"def u_ex(x, y):\n",
" return np.sin(np.pi * x) * np.sin(np.pi * y)\n",
"\n",
"\n",
"def u_d(x, y):\n",
" return np.zeros(x.shape[0])"
]
},
{
"cell_type": "code",
"execution_count": 75,
"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",
"execution_count": 76,
"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}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Résolution du système linéaire par une méthode directe\n"
]
},
{
"cell_type": "code",
"execution_count": 77,
"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",
" zeros_A.append(len(np.where(A == 0)[0]))\n",
" zeros_L.append(len(np.where(L == 0)[0]))"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7faf9ed53a30>"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAERCAYAAACepNcKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtWElEQVR4nO3deXhU1f3H8ffJThYChE0IkrCJCIKC7CqIWnfcxbrgikurttpWra3W2tr2Z6tSW9u6ohVBxQVQlLJWWWXf97AkkLAkJGTf5vz+uEPJSgJkciczn9fzzJOZe+7c+z1zw5eTc8+cY6y1iIhIYApxOwAREfEdJXkRkQCmJC8iEsCU5EVEApiSvIhIAFOSFxEJYH6X5I0x7xhjDhhj1tdj31eMMau9j63GmOxGCFFEpMkw/jZO3hhzAZAHvG+t7X0C73sEOMdae4/PghMRaWL8riVvrf0WyKq4zRjT1RjzjTFmhTHmO2NMzxreeiswqVGCFBFpIsLcDqCe3gAetNZuM8YMAl4HLjpaaIzpDCQDc12KT0TEL/l9kjfGxAJDgU+MMUc3R1bZbQwwxVpb3pixiYj4O79P8jhdStnW2n7H2WcM8KPGCUdEpOnwuz75qqy1R4CdxpibAIyj79Fyb/98S2CxSyGKiPgtv0vyxphJOAn7DGNMmjHmXuA24F5jzBpgAzC6wlvGAJOtvw0TEhHxA343hFJERBqO37XkRUSk4fjVjdfWrVvbpKSkStvy8/OJiYlxJyA/oPqr/qp/8NYf6v4MVqxYccha26a2cr9K8klJSSxfvrzStvnz5zNixAh3AvIDqr/qr/qPcDsMV9X1GRhjdh/v/equEREJYEryIiIBTEleRCSA+VWffE2MMezcuZOioiK3Q6lVVFQUiYmJhIeHux2KiEglfp/kY2JiiIuLIykpiQpz1/gNay2ZmZmkpaWRnJzsdjgiIpX4fXdNaGgoCQkJfpngwflLIyEhwa//0hCR4OX3SR7w2wR/lL/HJyLBq0kkeRGRQDVvywEmLNxJSZnHJ8dXkq+nL774AmMMmzdvdjsUEQkgr83ZxoRFuwgL8U2PgJJ8PU2aNInhw4czaZJWGBSRhrE2LZuVe7K5c0gSIUry7snLy2PBggW8/fbbTJ482e1wRCRATFi0i5iIUG4ckOizc/j9EMqKnp++gY37jjToMXt1aM5zV5913H2mTp3KZZddRo8ePUhISGDFihX079+/QeMQkeByMLeYL9ekc+vATjSP8t13bNSSr4dJkyYxZswYAMaMGaMuGxE5ZZO+30NJuYc7hyb59DxNqiVfV4vbF7Kyspg7dy7r1q3DGEN5eTnGGF566SUNnRSRk1JS5uGDJbu5oEcburaJ9em51JKvw5QpU7jjjjvYvXs3u3btIjU1leTkZL777ju3QxORJuqbDRkcyC3mbh+34kFJvk6TJk3iuuuuq7TthhtuUJeNiJy0CQt3kpQQzYU9al3ro8E0qe4aN8ybN6/atkcffdSFSEQkEKxJdYZNPnd1L58Nm6xILXkRkUb03tFhk/19N2yyIiV5EZFGcjC3mC/XpnNj/0TifDhssiIleRGRRtJYwyYrUpIXEWkER4dNXtgIwyYrUpIXEWkEX69P50BuMXc1YiselORFRBrFe4t2NdqwyYqU5OshNrbx/rQSkcBzdNjk2KG+m22yNkryIiI+dtxhkwVZUJzns3MryYuI+NDB3GKmr91X+7DJeS/Ca+dCqW/WiW5a33j9+inIWNewx2zfBy7/Y8MeU0TEa9L3eygttzUPm8w/BKs+gD43QniUT86vlryIiI/UOWzy+zehrBCGPuKzGJpWS14tbhFpQo4Om/zTjUnVC0vy4fs34IwroM0ZPotBLXkRER+ZsGgXya1juLB7DcMmV02EwiwY9phPY1CSr4eCggISExP/93j55ZfdDklE/Nya1GxW7cnmziGdqw+bLC+Dxa9Bp0Fw+mCfxuHT7hpjzE+B+wALrAPuttb65hayD3k8HrdDEJEm5rjDJjd+Adl74DLfd0H7rCVvjOkIPAoMsNb2BkKBMb46n4iIvzg6bPKmAZ2qD5u0FhaOh4Tu0ONyn8fi6+6aMKCZMSYMiAb2+fh8IiKu+9+wySGdqxfu/C9krIVhj0KI73vMjbXWdwc35jHg90Ah8B9r7W017DMOGAfQrl27/pMnT65UHhcXR/fu3f160WxrLTt27CAnJ6fBj52XlxfU0yqo/qp/U6t/mcfys/8W0ikuhCcGVB/7fvaa54jJ382SwW9iQ+qeU76uz2DkyJErrLUDat3BWuuTB9ASmAu0AcKBL4Dbj/ee/v3726qWLVtmDx48aD0eT7Uyf+DxeOzBgwdtSkqKT44/b948nxy3qVD957kdgquaYv2/WJVmOz/5pZ27eX/1wn1rrH2uubXf/qXex6vrMwCW2+PkVV/eeL0Y2GmtPQhgjPkMGAp8cCIHyc/PJzc3l4MHD/ogxIYRFRVFYmLjLOUlIv7tuMMmF/0VImJhwD2NFo8vk/weYLAxJhqnu2YUsPxED2KtJTk5uaFjExFpcEeHTf6mpkW6D++G9Z/B4IegWYtGi8lnvf7W2qXAFGAlzvDJEOANX51PRMRtR4dN3lDTsMklr4MxMPjhRo3Jp+PkrbXPAc/58hwiIv7gQG4R09fu47ZBnasPmyzIgpXvQ5+bIb5jo8alb7yKiDSASUtTax82uewtKC3w6URktVGSFxE5RSVlHiYu3c2IM9rQpepsk6WFsPRf0P1SaNer0WNTkhcROUVHZ5scW9Oc8as/hIJDPp+IrDZK8iIip6jWYZOeclj0GnTsD52HuRKbkryIyCk4OmxybE2zTW6aDod3Oq14l761ryQvInIK3lqwk9jIsOrDJo9ORNaqC/S8yp3gUJIXETlpW/fn8uXafdwxpIZhk7sWwL6VzoiakFB3AkRJXkTkpI2fs43o8FDGnd+leuHC8RDTBvre2viBVaAkLyJyErZk5DJjXTp3DUuiZUxE5cL9G2D7LBj4AIQ3cydALyV5EZGT8Nc524iJCOP+mlrxi16D8Gg4797GD6wKJXkRkRO0OeMIX61L5+5hSbSIrtKKz0mDdZ/AuWMhupU7AVagJC8icoLGz95GXGQY9w2voRW/5B/OyJohjTsRWW2U5EVETsDGfUf4en0Gdw9PJj66yoiawsOwYgL0vgFanO5KfFUpyYuInIDxc7YSFxXGvcNrWOdi+TtQkues3+onlORFROppw74cZm7Yzz3DkolvVqUVX1oES/4JXUdB+z7uBFgDJXkRkXp6dfY24qLCuKemVvzayZB/wK9a8aAkLyJSL+v35jBr437uG96leive43GGTZ7WF5IvdCfAWijJi4jUw6uzt9I8Koy7hydVL9wyAzK3uzoRWW2U5EVE6rA2LZvZmw5w//ldaF51jhprYeGr0KIznDnalfiOR0leRKQOr87eRovocO4allS9cM8SSFvmTEQW6tNls0+KkryIyHGsTs1m7manFV9tpkmA//4JmrWCfrc1fnD1oCQvInIcr87eSsvo8JqX9tsxF1LmwflPQER0o8dWH0ryIiK1WLnnMPO3HOT+C7oQG1mlK8bjgVnPQfzpMPB+dwKsB//rQBIR8ROvzt5Gq5gIxg5Jql64fgpkrIXr3oCwyEaPrb7UkhcRqcGK3Yf5dutBxl3QhZiqrfiyYpj7gvPN1j43uRNgPaklLyJSg1dnbyUhJoI7h3SuXrjsLcjeA7d/BiH+3Vb27+hERFywfFcW3207xAMXdiE6okpbuDAbvn0JuoyAbqPcCO+EKMmLiFTx6uxttI6N4PbBNbTiF77qTCl88fONHtfJUJIXEalg2a4sFmw/xAMXdK3eis/Z6ywK0ucm6NDPlfhOlJK8iEgFr8zaSuvYyJpb8fNfBOuBi37V+IGdJCV5ERGvpSmZLNqRyYMXdqFZRGjlwgObYPWHcN590DLJlfhOhpK8iIjXK7O30iaullb87N9ARCyc/7NGj+tUKMmLiACLd2SyJCWLhy7sSlR4lVb8roWw9RsY/hOISXAlvpOlJC8iQc9ayyuzt9I2LpI
"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",
"plt.ylabel(\"#zéros\")\n",
"plt.grid()\n",
"plt.legend()"
2022-03-24 07:10:59 +00:00
]
2022-03-15 07:04:08 +00:00
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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",
2022-04-03 20:18:44 +00:00
"version": "3.10.4"
2022-03-15 07:04:08 +00:00
}
},
"nbformat": 4,
"nbformat_minor": 4
}
