421 lines
61 KiB
Plaintext
421 lines
61 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 1,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"import mpl_toolkits.mplot3d\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import scipy.sparse as scps\n",
|
||
|
|
"import scipy.sparse.linalg as ssl\n",
|
||
|
|
"import math"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"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",
|
||
|
|
" n (int): nombre de points par cote du care => Npts points de discretisation au total\n",
|
||
|
|
"\n",
|
||
|
|
" Returns:\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",
|
||
|
|
" 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",
|
||
|
|
" \"\"\"\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",
|
||
|
|
"execution_count": 3,
|
||
|
|
"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",
|
||
|
|
" elements3 : matrice a trois colonnes contenant les elements triangles de la discretisation, identifies par les indices de leurs trois sommets.\n",
|
||
|
|
" coordinates : matrice a deux colonnes contenant les coordonnes 2D des points de la discretisation.\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",
|
||
|
|
"\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": [
|
||
|
|
"**Partie I : maillage triangulaire et conditions de Dirichlet**\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 4,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"def f(x, y):\n",
|
||
|
|
" # return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\n",
|
||
|
|
" return np.ones(x.shape[0])\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])\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def g(x):\n",
|
||
|
|
" return np.cos(x)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 5,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAADyCAYAAACLfbNuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACEsElEQVR4nO29d3hkd3X//7rTm0Yzo96llbb36oZtSgBj3Kg2JAZTY0I3JIGEJ3EK5QGSfPPFlN8XnJhQjLEN2GCzAQw2GBt7d70rrXZXq97LjKb3en9/jO7dGWmatJJ21zvv5/EDq7lz752Zez7nfM55n/cRRFGkjDLKePlDcaFvoIwyylgflI29jDIuE5SNvYwyLhOUjb2MMi4TlI29jDIuE5SNvYwyLhOoirxersuVUcbaQ1iPi5Q9exllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigb+wWAKIrEYjESiQRlKe8y1gvFxCvKWGWkUilisRiRSET+m1KpRK1Wo1KpUCqVCMK6aBmUcZlBKOJZym5nlSCKIolEgkQigSAIxONx+e+iKJJKpWQjj0ajVFRUoNFoysZ/eWBdfuCyZ18HSGF7pkFLEAQBQRBQKBTysUNDQ7S3t2MwGICy5y9jdVA29jVGIpFgcnKSZDJJU1MTgiDI3jyX0UrGr1QqUSqVstcPh8Py8SqVSv6vbPxllIqysa8RMsP2VColh+/LRS7Pn0wmSSQS8jEqlUr2/AqFomz8ZeRE2djXAKlUing8LoftkjcvFYWOl84nYbHxC4KQ5fnLxl+GhLKxryIkw5OSb5I3zme8+UL55SCX8ScSCfkepKhApVKh0WjKxn8Zo2zsqwRRFInH4ySTySUGuNjYi3n75UYCi9+72PgnJycBaGhoKHv+yxhlY18FSLVzyVPnyrhfKPJM5v1ICb94PJ7l+aU9v1KpLBv/yxhlYz8PLK6dS2H7Yqzmnv18IWX6JeQyfinZp1Kpci5eZVyaKBv7CrG4dl7IIBYb7+zsLP39/ajVaiwWC1arlcrKyiVGuJrId3+5jD8WixGNRoF03kGtVsuev2z8ly7Kxr4CSEm4fGH7YkjGnkwm6evrIxaLsX//fkRRxOv1Mj8/z9DQEEqlEqvVKmfyLwQKGX9msi8z7C/j0kCZLrsMSGF7T08PXV1d6HS6kt7ndDqZnp4mEAjQ1NRES0tL1mIhIRqN4vF4GBkZQRRFDAYDVqsVq9WKyWRasUednJxEoVDQ2Ni4ovdLkJ6VVCrF6dOn2bJlixzul43/vFCmy15MyKydS8m4UuF0OnE4HBw4cACz2Zz3OK1WS11dHR6Ph/r6ejQaDR6Ph4mJCfx+f5bxGwyGdQ+nMxN9kUhEjgBisRixWAyg7PkvYpSNvQgW186lULYUY08kEpw5c4ZIJEJtbW1BQ8+EFPbr9Xr0ej0NDQ2Iokg4HMbtdjMyMkIwGMRkMsnGr9frz+tzrgSZ1F4419Sz2Pgzef1l479wKBt7AeSrnQuCUHRP7ff7OXnyJK2trRiNRqanp8/rXgRBwGAwYDAYaGpqQhRFgsEgbrebgYEBIpEIFRUVsvFrtdrzul6p97T434tr/KIoEo1GlyT8ysa//igbex4Uqp0XKo2JosjU1BTj4+Ps3LmTiooKvF7vqpfeBEHAZDJhMploaWkhlUrh9/txu92cPn2aeDxOZWUlVquVZDJ5QYyqmPGLopgV8kulvjLWBmVjX4TMsD1f7VyhUOT07IlEglOnTqFQKDh06BAqVfrrXQ9SjUKhoLKyksrKStrb20mlUni9XtxuN3Nzc3IkYLVasVgs8r2tJ3IZfyqVkoU8pqenaW5uRqPRlDv61gBlY89AqbXzXMbr8/no7e2lvb19SdY71/FrRZeVoFAo5JBeo9EAoNfrcbvdjI6OIghC3hr/emHxdzw7O0tjY2NZxWeNUDb2BRSjvGYiM0EniiITExNMTU2xe/dujEbjkuMvJF1WgkKhoKqqiqqqKgDi8Tgej2dJjd9qtWI2my/YXnpxjX9xL3/Z+FeOy97YS6W8ZkJK0MXjcXp7e9FqtRw6dCivd7yY6LIS1Go1NTU11NTUAOnymdvtltl9Wq12VWr854Ncvfxl4185LmtjT6VS2O12RFHEarWW/KAIgoDf76evr4/Ozk7q6+uLHn+hPXuxz6bRaKirq6Ourg6ASCSC2+3OW+O/ECjF+MsqPvlxWRp7ZhLO7/cjiiI2m63k9/p8PiKRCPv37y/pwb/Qnn0l59LpdDQ0NOSt8UejUaanp7Fareh0uovG82cmVyVoNBq0Wu1l39F32Rn74rBdqVTKhJliiMVinDx5kkQiQVdXV8ke7mLw7OeDXDX+F154gWQyecFq/Pnuc7Hxj4+Po1arqa2tzWrnvRx7+S8rY88lF5WvjLYYUv1648aNBAKBZT0kF9qzrzakRbKlpUWu8QcCAVwu15Iav8VikasBF+I+pVq+SqW67CW8LgtjL1Q7L0Z9FUWR4eFh5ufn2bdvH3q9nmAwuKyutJebsS+GQqHAbDZjNpuX1PglZV2pzLfeNf5UKpUlD1ZMwuvlbPwve2MvJBcFhamv0WiUkydPUlFRwcGDB4tqyuXDpWa854vMGj+kW4I9Hs8FqfFnGvti5DL+l7OKz8va2EupnecL451OJ319fWzatEkuTxV7Tz4sNvZgMEhPTw+CIGC1WrHZbFRUVKx4MbnYoVQqS6rxS7Lbq1njX875lqPicyka/8vS2JdTO19suKIoMjg4iMfjYf/+/Tl71s/Hs8/OzjI8PMyWLVtQq9V4PB6mp6fx+/3odDqsViuxWOyCJbnWA7lq/C6Xi0QiwdGjR9FoNHJkUFFRcV4GdT6LRyEhj+npaerq6jAYDJeMhNfLztiXIxcF2cYeiUTo6enBZrNx4MCBgnTZZDJZ8j1JWwWp3VXaEiQSCbm2nVnecrlc2O12XC6X7Plfzsav0WioqalhamqK/fv3yzX+yclJAoGAvAharVaMRuOyDGo1I4VM43e5XNTV1WWp+Eie/2Lt5X9ZGbsUcklffCkPhWTsDoeD/v5+tmzZIoebhd6zHM8eiUQIBoM0NjayZcsWBEHImugC2eWtZDKJWq3GZDLJGe5EIiFnuK1W6wVpZFlLZBplvhr/6OgowWAQo9GY1cdf6Hde7W2BhGQyKY/ognNcBqmX/yMf+Qif+9zn2LJly6pfe6V4WTwxUtje19eHxWKhtrZ2We/3eDxEo1EOHDhQkgddThgvLSJarZb29vaS3iM9vBUVFVRUVNDW1kYymZQz3GNjY/J+X0pyFXugL+bwEvIPzCjUxz84OEgkEpFFPHJFQGtl7KlUKivEz6TvQtrzXwhBkUK45I09s3a+3MRZOBzm1KlTCILA/v37SzaIUq4j7f29Xi8HDx7k6NGjJd9XLiiVSmw2m8z0i8fjuN1u7HY7AwMDaDQabDbbunDZ1yJ5KIpiyX0Ji/v4A4EAbrebM2fOEIvFsmr8a2Xs0r3kg6QkdDHhkjX2XKOWlEplycYuGUlnZyfT09OrSpKJRqP09PRgtVoLLiKFvFkxg5JYYVIUE4lEcLlcjI+PEwgEMBqNsvGvtnGuxtiqxcg1zroUZNb429raSKVS+Hw+ec/v9/sZGRmhqqpqVWv8xb7TUChUNvbVQL7aeSkeN5VKcfbsWcLhMAcPHswaj1QqCtXmJabd5s2bqa6uXtZ5zwc6nY7GxkYaGxvlUNflctHf34/f78doNMr1b7VavW73VSpWawFRKBRYLBYsFgsdHR0cOXKEqqoqvF4vY2NjAOtS44/H4xeMOZgPl5yxF6qdK5XKglnyUChET08P9fX1cqJsJRrtuRJ0oigyOjqK3W7PW7IrFedbZ88MdVtbWxkbGyORSBAIBJiYmEAURSwWCzabbUUP/Fp59rUKt6urq+WFN5FI4PF4cDqd59XHf7HnQHLhkjH2UmrnCoUib1PL7OwsQ0ND7Ni
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 432x288 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"n = 10\n",
|
||
|
|
"coords, elems3, dirichlet, neumann = maillage_carre(n)\n",
|
||
|
|
"show(coords, f(coords[:, 0], coords[:, 1]))\n",
|
||
|
|
"# show(coords, u_ex(coords[:, 0], coords[:, 1]))\n",
|
||
|
|
"# print(maillage_carre(3))"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 6,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"def raideur(triangle):\n",
|
||
|
|
" M = np.zeros((3, 3))\n",
|
||
|
|
" x = triangle[:, 0]\n",
|
||
|
|
" y = triangle[:, 1]\n",
|
||
|
|
"\n",
|
||
|
|
" # calcul de alpha\n",
|
||
|
|
" 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",
|
||
|
|
" alpha = np.linalg.det(mat_alpha)\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
" 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",
|
||
|
|
" return M / alpha / 2"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 7,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"ename": "IndexError",
|
||
|
|
"evalue": "tuple index out of range",
|
||
|
|
"output_type": "error",
|
||
|
|
"traceback": [
|
||
|
|
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
|
||
|
|
"\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)",
|
||
|
|
"\u001b[1;32m/home/laurent/Documents/Cours/ENSEEIHT/S8 - Équations aux Dérivées Partielles/TP-EDP.ipynb Cell 8'\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=59'>60</a>\u001b[0m coords, elems3, dirichlet, neumann \u001b[39m=\u001b[39m maillage_carre(n)\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=61'>62</a>\u001b[0m A \u001b[39m=\u001b[39m assemblage(coords, elems3)\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=62'>63</a>\u001b[0m b \u001b[39m=\u001b[39m second_membre(coords, elems3, f)\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=63'>64</a>\u001b[0m U_d \u001b[39m=\u001b[39m calcul_Ud(coords, dirichlet)\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=64'>65</a>\u001b[0m b \u001b[39m-\u001b[39m\u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mdot(A, U_d)\n",
|
||
|
|
"\u001b[1;32m/home/laurent/Documents/Cours/ENSEEIHT/S8 - Équations aux Dérivées Partielles/TP-EDP.ipynb Cell 8'\u001b[0m in \u001b[0;36msecond_membre\u001b[0;34m(coords, elem3, f)\u001b[0m\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=21'>22</a>\u001b[0m mat_alpha \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39marray([[x[\u001b[39m1\u001b[39m] \u001b[39m-\u001b[39m x[\u001b[39m0\u001b[39m], x[\u001b[39m2\u001b[39m] \u001b[39m-\u001b[39m x[\u001b[39m0\u001b[39m]], [y[\u001b[39m1\u001b[39m] \u001b[39m-\u001b[39m y[\u001b[39m0\u001b[39m], y[\u001b[39m2\u001b[39m] \u001b[39m-\u001b[39m y[\u001b[39m0\u001b[39m]]])\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=22'>23</a>\u001b[0m alpha \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mlinalg\u001b[39m.\u001b[39mdet(mat_alpha)\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=24'>25</a>\u001b[0m b[triangle] \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m alpha \u001b[39m/\u001b[39m \u001b[39m6\u001b[39m \u001b[39m*\u001b[39m f(centre[\u001b[39m0\u001b[39;49m], centre[\u001b[39m1\u001b[39;49m])\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000007?line=26'>27</a>\u001b[0m \u001b[39mreturn\u001b[39;00m b\n",
|
||
|
|
"\u001b[1;32m/home/laurent/Documents/Cours/ENSEEIHT/S8 - Équations aux Dérivées Partielles/TP-EDP.ipynb Cell 5'\u001b[0m in \u001b[0;36mf\u001b[0;34m(x, y)\u001b[0m\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000004?line=0'>1</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mf\u001b[39m(x, y):\n\u001b[1;32m <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000004?line=1'>2</a>\u001b[0m \u001b[39m# return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\u001b[39;00m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/laurent/Documents/Cours/ENSEEIHT/S8%20-%20%C3%89quations%20aux%20D%C3%A9riv%C3%A9es%20Partielles/TP-EDP.ipynb#ch0000004?line=2'>3</a>\u001b[0m \u001b[39mreturn\u001b[39;00m np\u001b[39m.\u001b[39mones(x\u001b[39m.\u001b[39;49mshape[\u001b[39m0\u001b[39;49m])\n",
|
||
|
|
"\u001b[0;31mIndexError\u001b[0m: tuple index out of range"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"def assemblage(coords, elems3):\n",
|
||
|
|
" Ns = len(coords)\n",
|
||
|
|
" A = np.zeros((Ns, Ns))\n",
|
||
|
|
" for triangle in elems3:\n",
|
||
|
|
" M = raideur(coords[triangle])\n",
|
||
|
|
" for i, a in enumerate(triangle):\n",
|
||
|
|
" for j, b in enumerate(triangle):\n",
|
||
|
|
" A[a, b] += M[i, j]\n",
|
||
|
|
" return A\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def second_membre(coords, elem3, f):\n",
|
||
|
|
" Ns = len(coords)\n",
|
||
|
|
" b = np.zeros(Ns)\n",
|
||
|
|
" for triangle in elem3:\n",
|
||
|
|
" coords_triangle = coords[triangle]\n",
|
||
|
|
" centre = np.mean(coords_triangle, 0)\n",
|
||
|
|
"\n",
|
||
|
|
" # calcul de alpha\n",
|
||
|
|
" x = coords_triangle[:, 0]\n",
|
||
|
|
" y = coords_triangle[:, 1]\n",
|
||
|
|
" mat_alpha = np.array([[x[1] - x[0], x[2] - x[0]], [y[1] - y[0], y[2] - y[0]]])\n",
|
||
|
|
" alpha = np.linalg.det(mat_alpha)\n",
|
||
|
|
"\n",
|
||
|
|
" b[triangle] += alpha / 6 * f(centre[0], centre[1])\n",
|
||
|
|
"\n",
|
||
|
|
" return b\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def calcul_Ud(coords, dirichlet):\n",
|
||
|
|
" Ns = len(coords)\n",
|
||
|
|
" U = np.zeros(Ns)\n",
|
||
|
|
" # for d in dirichlet:\n",
|
||
|
|
" # x, y = coords[d].flatten()\n",
|
||
|
|
" # U[d] = u_d(x, y)\n",
|
||
|
|
" U[dirichlet.T] = u_d(coords[dirichlet, 0], coords[dirichlet, 1])\n",
|
||
|
|
"\n",
|
||
|
|
" return U\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def tildage(A, b, coords, dirichlet):\n",
|
||
|
|
" A_tild = np.delete(A, dirichlet, 0)\n",
|
||
|
|
" A_tild = np.delete(A_tild, dirichlet, 1)\n",
|
||
|
|
" b_tild = np.delete(b, dirichlet, 0)\n",
|
||
|
|
" coords_tild = np.delete(coords, dirichlet, 0)\n",
|
||
|
|
"\n",
|
||
|
|
" return A_tild, b_tild, coords_tild\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"def untildage(x, dirichlet, U_d):\n",
|
||
|
|
" x_untild = np.zeros(U_d.shape[0])\n",
|
||
|
|
" not_dirichlet = np.setdiff1d(range(n*n), dirichlet)\n",
|
||
|
|
"\n",
|
||
|
|
" x_untild[dirichlet] = U_d[dirichlet]\n",
|
||
|
|
" x_untild[not_dirichlet] = x\n",
|
||
|
|
"\n",
|
||
|
|
" return x_untild\n",
|
||
|
|
"\n",
|
||
|
|
"n = 50\n",
|
||
|
|
"coords, elems3, dirichlet, neumann = maillage_carre(n)\n",
|
||
|
|
"\n",
|
||
|
|
"A = assemblage(coords, elems3)\n",
|
||
|
|
"b = second_membre(coords, elems3, f)\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",
|
||
|
|
"\n",
|
||
|
|
"# show(coords_tild, x)\n",
|
||
|
|
"# print(coords.shape, x_untild.shape)\n",
|
||
|
|
"# print(coords, x_untild)\n",
|
||
|
|
"show(coords, x_untild)\n",
|
||
|
|
"show(coords, u_ex(coords[:, 0], coords[:, 1]))"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"**Partie II : maillage mixte et ajoût des conditions de Neumann**\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"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",
|
||
|
|
" [0.0, 0.0],\n",
|
||
|
|
" [1 / 3, 0],\n",
|
||
|
|
" [0.53333333333333, 0.0],\n",
|
||
|
|
" [2 / 3, 1 / 3],\n",
|
||
|
|
" [1.0, 0.47],\n",
|
||
|
|
" [1, 2 / 3],\n",
|
||
|
|
" [1.0, 1.0],\n",
|
||
|
|
" [2 / 3, 1.0],\n",
|
||
|
|
" [1 / 3, 1.0],\n",
|
||
|
|
" [0.0, 1.0],\n",
|
||
|
|
" [0.0, 2 / 3],\n",
|
||
|
|
" [0.0, 1 / 3],\n",
|
||
|
|
" [1 / 3, 1 / 3],\n",
|
||
|
|
" [1 / 3, 2 / 3],\n",
|
||
|
|
" [2 / 3, 2 / 3],\n",
|
||
|
|
" [1.0, 0.0],\n",
|
||
|
|
" ]\n",
|
||
|
|
")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": []
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"**Compléments d'analyse du système**\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": []
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": []
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"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",
|
||
|
|
"version": "3.10.2"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 4
|
||
|
|
}
|