From eeb83fe2682e9451d363ac70fafc47d3af9d73e8 Mon Sep 17 00:00:00 2001 From: Paul Brinkmeier Date: Wed, 21 Jun 2023 09:56:19 +0200 Subject: [PATCH] Fix broken equation --- experiments/ELM.ipynb | 133 +++++++++++++++++++++++++----------- experiments/ELM_Debug.ipynb | 59 +++++++++------- 2 files changed, 128 insertions(+), 64 deletions(-) diff --git a/experiments/ELM.ipynb b/experiments/ELM.ipynb index 759e710..d36322d 100644 --- a/experiments/ELM.ipynb +++ b/experiments/ELM.ipynb @@ -16,7 +16,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 72, "id": "b39df788-2971-4974-9ced-1518a47181e8", "metadata": {}, "outputs": [], @@ -33,16 +33,32 @@ " # c\n", " bond_bending_constant: float\n", " \n", - " def __init__(self, n_rows: int, n_cols: int):\n", + " def make(n_rows: int, n_cols: int) -> ELM:\n", + " return ELM(\n", + " n_rows,\n", + " n_cols,\n", + " numpy.zeros((n_rows, n_cols, 2)),\n", + " numpy.zeros((n_rows, n_cols, 2)),\n", + " numpy.zeros((n_rows, n_cols, 2))\n", + " )\n", + " \n", + " def __init__(\n", + " self,\n", + " n_rows: int,\n", + " n_cols: int,\n", + " displacement: numpy.ndarray,\n", + " velocity: numpy.ndarray,\n", + " force: numpy.ndarray\n", + " ) -> None:\n", " self.n_rows = n_rows\n", " self.n_cols = n_cols\n", - " self.displacement = numpy.zeros((n_rows, n_cols, 2))\n", - " self.velocity = numpy.zeros((n_rows, n_cols, 2))\n", - " self.force = numpy.zeros((n_rows, n_cols, 2))\n", + " self.displacement = displacement\n", + " self.velocity = velocity\n", + " self.force = force\n", "\n", " self.elastic_spring_constant = 1.0\n", - " self.bond_bending_constant = 1.0\n", - "\n", + " self.bond_bending_constant = 2.0\n", + " \n", " # u_i(t)\n", " def get_displacement(self, row: int, col: int) -> numpy.ndarray:\n", " return self.displacement[row, col, :]\n", @@ -102,39 +118,72 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 77, "id": "b4ed52bd-02aa-411a-9af4-be9bb28227bc", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAk1UlEQVR4nO3dfXCU1d3/8U9CyAMP2ZA47CY1gdQyg4gPPChEnLa35C6oQ0EZW5zYwYeBqkFFWlGq4FDFKG3VYhEqo4BTlOqMoDgVhztYrGNIIIAVsYg1P0iFDVXMLqAEzJ7fHwyrFydogps918b3a2ZnzHWd7J4cR8+Z/e7nu2nGGCMAAAAfSXc9AQAAgJNxQAEAAL7DAQUAAPgOBxQAAOA7HFAAAIDvcEABAAC+wwEFAAD4DgcUAADgOxmuJ3A6YrGY9u7dq969eystLc31dAAAQDsYY3Tw4EEVFRUpPf3r3yNJyQPK3r17VVxc7HoaAADgNDQ2NurMM8/82jEpeUDp3bu3pON/YG5uruPZAACA9ohGoyouLo7v418nJQ8oJ8o6ubm5HFAAAEgx7fl4Bh+SBQAAvsMBBQAA+A4HFAAA4DscUAAAgO9wQAEAAL7DAQUAAPhOhw8ob7zxhsaNG6eioiKlpaVp9erVnvvGGM2ZM0eFhYXKyclReXm5du3a5Rlz4MABVVRUKDc3V3l5ebrxxht16NChb/WHAACArqPDB5TDhw/r/PPP18KFC9u8P3/+fC1YsECLFy9WbW2tevbsqTFjxujIkSPxMRUVFXr33Xe1bt06vfLKK3rjjTc0derU0/8rAABAl5JmjDGn/ctpaVq1apUmTJgg6fi7J0VFRfrVr36lX//615KkSCSiYDCoZcuWadKkSXrvvfc0aNAgbdq0ScOHD5ckrV27Vpdffrn+85//qKio6BtfNxqNKhAIKBKJdHqjttaYUV3DAe0/eER9e2drWL8+qt/9afzni0rzJYkxjGEMYxjDGM8Y16+fyL+jW3pivveuI/t3QjvJNjQ0KBwOq7y8PH4tEAhoxIgRqqmp0aRJk1RTU6O8vLz44USSysvLlZ6ertraWl155ZXW87a0tKilpSX+czQaTeS0T2nt9n2au2aH9kW+fPcnPU2KfeVIl9ejuySp+bNjjGEMYxjDGMb45vUTNaYwkK37xg3S2MGFSqaEHlDC4bAkKRgMeq4Hg8H4vXA4rL59+3onkZGh/Pz8+JiTVVVVae7cuYmc6jdau32fbv7LFp389lLspAtf/ZfIGMYwhjGMYYxfXj9RY8KRI7r5L1u06NqhST2kpESKZ9asWYpEIvFHY2Njp75ea8xo7pod1uEEAIDvmhN74dw1O9R68ommEyX0gBIKhSRJTU1NnutNTU3xe6FQSPv37/fc/+KLL3TgwIH4mJNlZWXFvxgwGV8QWNdwwFPWAQDgu8xI2hc5orqGA0l7zYQeUEpLSxUKhVRdXR2/Fo1GVVtbq7KyMklSWVmZmpubVV9fHx+zfv16xWIxjRgxIpHTOW37D3I4AQDgZMncHzv8GZRDhw7pgw8+iP/c0NCgbdu2KT8/XyUlJZo+fboeeOABDRgwQKWlpZo9e7aKioriSZ+zzz5bY8eO1ZQpU7R48WIdO3ZM06ZN06RJk9qV4EmGvr2zXU8BAADfSeb+2OEDyubNm/U///M/8Z9nzJghSZo8ebKWLVummTNn6vDhw5o6daqam5t1ySWXaO3atcrO/vKPWrFihaZNm6bRo0crPT1dEydO1IIFCxLw5yTGsH59rE82AwDwXZaednx/TJZv1QfFlc7ug1Lz7090zZKNCX9eAABS2XNTRqrsrILT/v2O7N8pkeJJNj6DAgCALZn7IweUNpzRK8v1FAAA8J1k7o8cUNqSckUvAACSIIn7IweUNnx8uOWbBwEA8B2TzP2RA0obiBkDAGBL5v7IAaUNJ2LGAADguGTHjDmgtKF+96f0QAEA4Cti5vj+mCwcUNpAzBgAABsxY8eIGQMAYCNm7BrlHQAAbMSM3SJmDACAjZixY5R4AACwUeJxjRIPAAA2SjxuUeIBAMBGiccxOskCAGCjk6xjdJIFAMCLTrI+QCdZAAC86CTrA3SSBQDARidZx4gZAwBgI2bsGuUdAABsxIzdImYMAICNmLFjxIwBALARM3aMmDEAAF7EjH2AmDEAAF7EjH2AmDEAADZixo4RMwYAwEbM2DXKOwAA2IgZu0XMGAAAGzFjxyjxAABgo8TjGiUeAABslHjcosQDAICNEo9jdJIFAMBGJ1nH6CQLAIAXnWR9gE6yAAB40UnWB+gkCwCAjU6yjhEzBgDARszYNco7AADYiBm7RcwYAAAbMWPHiBkDAGAjZuwYMWMAALyIGfsAMWMAALyIGfsAMWMAAGzEjB0jZgwAgI2YsWuUdwAAsBEzdouYMQAANmLGjlHiAQDARonHNUo8AADYKPG4RYkHAAAbJR7H6CQLAICNTrKO0UkWAAAvOsn6AJ1kAQDwopOsD9BJFgAAG51kHSNmDACAjZixa5R3AACwETN2i5gxAAA2YsaOETMGAMBGzNgxYsYAAHilfMy4tbVVs2fPVmlpqXJycnTWWWfp/vvvlzFfFq6MMZozZ44KCwuVk5Oj8vJy7dq1K9FTOW3EjAEA8Er5mPHDDz+sRYsW6U9/+pPee+89Pfzww5o/f74ef/zx+Jj58+drwYIFWrx4sWpra9WzZ0+NGTNGR474I95LzBgAAFsy98eMRD/hW2+9pfHjx+uKK66QJPXv31/PPfec6urqJB1/9+Sxxx7Tvffeq/Hjx0uSnnnmGQWDQa1evVqTJk1K9JQ6jJgxAAC2lI4ZX3zxxaqurtb7778vSXr77bf15ptv6rLLLpMkNTQ0KBwOq7y8PP47gUBAI0aMUE1NTZvP2dLSomg06nl0Kso7AADYkrg/JvwdlLvvvlvRaFQDBw5Ut27d1Nraqnnz5qmiokKSFA6HJUnBYNDze8FgMH7vZFVVVZo7d26ip3pKxIwBALCldMz4+eef14oVK/Tss89qy5YtWr58uX7/+99r+fLlp/2cs2bNUiQSiT8aGxsTOGMbJR4AAGzJ3B8T/g7KnXfeqbvvvjv+WZJzzz1Xu3fvVlVVlSZPnqxQKCRJampqUmFhYfz3mpqadMEFF7T5nFlZWcrKSuKhgRIPAAC2VO4k+9lnnyk93fu03bp1UywWkySVlpYqFAqpuro6fj8ajaq2tlZlZWWJns5pocQDAIAtmftjwt9BGTdunObNm6eSkhKdc8452rp1qx555BHdcMMNkqS0tDRNnz5dDzzwgAYMGKDS0lLNnj1bRUVFmjBhQqKnc1roJAsAgC2Z+2PCDyiPP/64Zs+erVtuuUX79+9XUVGRfvnLX2rOnDnxMTNnztThw4c1depUNTc365JLLtHatWuVne2Pg8GJTrI0awMA4Lhkd5JNM19t8ZoiotGoAoGAIpGIcnNzE/78Nf/+RNcs2Zjw5wUAIJU9N2Wkys4qOO3f78j+zXfxtIFOsgAA2JK5P3JAaQMxYwAAbCndSbZLSLmiFwAASZDKMeOugJgxAAC2lO4k2xUQMwYAwJbM/ZEDShtOxIwBAMBxyY4Zc0BpQ/3uT+mBAgDAV8TM8f0xWTigtIGYMQAANmLGjhEzBgDARszYNco7AADYiBm7RcwYAAAbMWPHKPEAAGCjxOMaJR4AAGyUeNyixAMAgI0Sj2N0kgUAwEYnWcfoJAsAgBedZH2ATrIAAHjRSdYH6CQLAICNTrKOETMGAMBGzNg1yjsAANiIGbtFzBgAABsxY8eIGQMAYCNm7BgxYwAAvIgZ+wAxYwAAvIgZ+wAxYwAAbMSMHSNmDACAjZixa5R3AACwETN2i5gxAAA2YsaOUeIBAMBGicc1SjwAANgo8bhFiQcAABslHsfoJAsAgI1Oso7RSRYAAC86yfoAnWQBAPCik6wP0EkWAAAbnWQdI2YMAICNmLFrlHcAALARM3aLmDEAADZixo5R4gEAwEaJxzVKPAAA2CjxuEWJBwAAGyUex+gkCwCAjU6yjtFJFgAALzrJ+gCdZAEA8KKTrA/QSRYAABudZB0jZgwAgI2YsWuUdwAAsBEzdouYMQAANmLGjhEzBgDARszYMWLGAAB4ETP2AWLGAAB4ETP2AWLGAADYiBk7RswYAAAbMWPXKO8AAGAjZuwWMWMAAGzEjB2jxAMAgI0Sj2uUeAAAsKV6ieejjz7Stddeq4KCAuXk5Ojcc8/V5s2b4/eNMZozZ44KCwuVk5Oj8vJy7dq1qzOmcloo8QAAYEvpEs+nn36qUaNGqXv37nr11Ve1Y8cO/eEPf1CfPl82d5k/f74WLFigxYsXq7a2Vj179tSYMWN05Ig/4r10kgUAwJbM/TEj0U/48MMPq7i4WEuXLo1fKy0tjf+zMUaPPfaY7r33Xo0fP16S9MwzzygYDGr16tWaNGlSoqfUYSc6ydKsDQCA41K+k+zLL7+s4cOH6+qrr1bfvn01ZMgQLVmyJH6/oaFB4XBY5eXl8WuBQEAjRoxQTU1Nm8/Z0tKiaDTqeXQmOskCAOCV8p1kP/zwQy1atEgDBgzQa6+9pptvvlm33Xabli9fLkkKh8OSpGAw6Pm9YDAYv3eyqqoqBQKB+KO4uDjR0/agkywAALaU7iQbi8U0dOhQPfjggxoyZIimTp2qKVOmaPHixaf9nLNmzVIkEok/GhsbEzhjGzFjAABsKR0zLiws1KBBgzzXzj77bO3Zs0eSFAqFJElNTU2eMU1NTfF7J8vKylJubq7n0ako7wAAYEvlmPGoUaO0c+dOz7X3339f/fr1k3T8A7OhUEjV1dXx+9FoVLW1tSorK0v0dE4LMWMAAGzJ3B8TnuK54447dPHFF+vBBx/Uz372M9XV1enJJ5/Uk08+KUlKS0vT9OnT9cADD2jAgAEqLS3V7NmzVVRUpAkTJiR6OqeFmDEAALaUjhlfeOGFWrVqlWbNmqXf/va3Ki0t1WOPPaaKior4mJkzZ+rw4cOaOnWqmpubdckll2jt2rXKzvbHwYCYMQAAXsmOGacZY1JuG45GowoEAopEIp3yeZSaf3+ia5ZsTPjzAgCQyp6bMlJlZxWc9u93ZP/mu3jaQMwYAABbSseMuwJixgAA2FI6ZtwlpFzRCwCAJEjlmHFXQMwYAABbSn+bcVdAiQcAABslHtco8QAAYKPE4xYlHgAAbJR4HKOTLAAAtmTujxxQ2nCikywAADgu2Z1kOaC0oX73p7S5BwDgK2Lm+P6YLBxQ2vB/O8KupwAAgO/QSdah1pjRqm0fuZ4GAAC+Q8zYobqGAzpw+JjraQAA4D/EjN3hiwIBAGgbMWOH/t/Hh11PAQAAXyJm7EhrzOi5uj2upwEAgO8QM3aoruGAwlG6yAIAcDJixg4RLwYA4NSIGTvQGjN6cSvxYgAAToWYsQN/Wr9Ln35GvBgAgFMiZpxca7fv06P/t8v1NAAA8DVixknUGjOau2aH62kAAOB7lHiSqK7hgPZFaM4GAMA3osSTPHSOBQCgfSjxJFEyu+IBAJDK6CSbRMP69VF62jePa88YAAC6KjrJJln97k8Va0dNrT1jAADoqugkm2R8BgUAgPahk2wSJTMyBQBAKiNmnEyUbgAAaB9ixsmTzMgUAACpjJhxEhEzBgCgfYgZJ9FFpfkqDGSLFDEAAKdGzDjJuqWn6b5xg/goCgAAX4OYsQNjBxfqhlH9XU8DAABfI2bswP8OCrmeAgAAvkbM2IGLSvOV3zPT9TQAAPAvYsbJ1y09TRMuKHI9DQAAfIuYsSOUeQAAODVKPI5Q5gEA4GtQ4nGDMg8AAKdGicchyjwAALSNTrIODevXR+m0lQUAwINOso7V7/5UMdrKAgDgQSdZx5LZJQ8AgFRCJ1mH+HZjAADaRszYIaLGAACcAjFjd4gaAwDQNmLGjhE1BgDARszYMaLGAAB4ETP2AaLGAAB4ETP2AaLGAADYiBk7lswYFQAAqYKYsWuUdwAAsBEzdiuZMSoAAFIFMWPHKPEAAGCjxOMaJR4AAGxdqcTz0EMPKS0tTdOnT49fO3LkiCorK1VQUKBevXpp4sSJampq6uyptBslHgAAbF2mxLNp0yb9+c9/1nnnnee5fscdd2jNmjV64YUXtGHDBu3du1dXXXVVZ06lQ/jCQAAAbF2ik+yhQ4dUUVGhJUuWqE+fLzvPRSIRPfXUU3rkkUd06aWXatiwYVq6dKneeustbdy4sbOm0yF0kgUAwKvLdJKtrKzUFVdcofLycs/1+vp6HTt2zHN94MCBKikpUU1NTZvP1dLSomg06nl0JjrJAgDglexOshmd8aQrV67Uli1btGnTJuteOBxWZmam8vLyPNeDwaDC4XCbz1dVVaW5c+d2xlTbRCdZAABsKd1JtrGxUbfffrtWrFih7OzE1KpmzZqlSCQSfzQ2NibkeU+FmDEAALaUjhnX19dr//79Gjp0qDIyMpSRkaENGzZowYIFysjIUDAY1NGjR9Xc3Oz5vaamJoVCoTafMysrS7m5uZ5Hp6K8AwCALYn7Y8JLPKNHj9Y777zjuXb99ddr4MCBuuuuu1RcXKzu3bururpaEydOlCTt3LlTe/bsUVlZWaKnc1qIGQMAYEvm/pjwA0rv3r01ePBgz7WePXuqoKAgfv3GG2/UjBkzlJ+fr9zcXN16660qKyvTyJEjEz2d00LMGAAAWzL3x075kOw3efTRR5Wenq6JEyeqpaVFY8aM0RNPPOFiKm06ETMmyQMAwHHJjhmnGWNSbhuORqMKBAKKRCKd8nmUmn9/omuW+KMnCwAAfvHclJEqO6vgtH+/I/s338XTBmLGAADYUjpm3BUQMwYAwJbSMeMuIeWKXgAAJEFX+jbjVETMGAAAW5f5NuNURYkHAAAbJR7XKPEAAGCjxOMWJR4AAGyUeByjkywAALZk7o8cUNpwopMsAAA4LtmdZDmgtKF+96e0uQcA4Cti5vj+mCwcUNpAJ1kAAGx0knWMmDEAADZixq5R3gEAwEbM2C1ixgAA2IgZO0aJBwAAGyUe1yjxAABgo8TjFiUeAABslHgco5MsAAA2Osk6RidZAAC86CTrA3SSBQDAi06yPkAnWQAAbHSSdYyYMQAANmLGrlHeAQDARszYLWLGAADYiBk7RswYAAAbMWPHiBkDAOBFzNgHiBkDAOBFzNgHiBkDAGAjZuwYMWMAAGzEjF2jvAMAgI2YsVvEjAEAsBEzdowSDwAANko8rlHiAQDARonHLUo8AADYKPE4RidZAABsdJJ1jE6yAAB40UnWB+gkCwCAF51kfYBOsgAA2Ogk6xgxYwAAbMSMXaO8AwCAjZixW8SMAQCwETN2jJgxAAA2YsaOETMGAMCLmLEPEDMGAMCLmLEPEDMGAMBGzNgxYsYAANiIGbtGeQcAABsxY7eIGQMAYCNm7BglHgAAbJR4XKPEAwCAjRKPW5R4AACwUeJxjE6yAADY6CTrGJ1kAQDwopOsD9BJFgAALzrJ+gCdZAEAsNFJ1jFixgAA2FI6ZlxVVaULL7xQvXv3Vt++fTVhwgTt3LnTM+bIkSOqrKxUQUGBevXqpYkTJ6qpqSnRUzl9lHcAALClcsx4w4YNqqys1MaNG7Vu3TodO3ZMP/nJT3T48OH4mDvuuENr1qzRCy+8oA0bNmjv3r266qqrEj2V00bMGAAAWzL3x4xEP+HatWs9Py9btkx9+/ZVfX29fvjDHyoSieipp57Ss88+q0svvVSStHTpUp199tnauHGjRo4cmegpdRgxYwAAbF0qZhyJRCRJ+fn5kqT6+nodO3ZM5eXl8TEDBw5USUmJampqOns67ULMGAAAry4VM47FYpo+fbpGjRqlwYMHS5LC4bAyMzOVl5fnGRsMBhUOh9t8npaWFkWjUc+jMxEzBgDAq0vFjCsrK7V9+3atXLnyWz1PVVWVAoFA/FFcXJygGbaNmDEAALYuETOeNm2aXnnlFb3++us688wz49dDoZCOHj2q5uZmz/impiaFQqE2n2vWrFmKRCLxR2NjY2dNWxIxYwAA2pLSMWNjjKZNm6ZVq1Zp/fr1Ki0t9dwfNmyYunfvrurq6vi1nTt3as+ePSorK2vzObOyspSbm+t5dCrKOwAA2JK4PyY8xVNZWalnn31WL730knr37h3/XEkgEFBOTo4CgYBuvPFGzZgxQ/n5+crNzdWtt96qsrIyXyR4JGLGAAC0JaVjxosWLZIk/fjHP/ZcX7p0qa677jpJ0qOPPqr09HRNnDhRLS0tGjNmjJ544olET+W0UeIBAMCWzP0x4QcUY775/Z/s7GwtXLhQCxcuTPTLJwYlHgAAbKncSbYroMQDAIAtmfsjB5Q20EkWAABbl+okm4roJAsAgFeX6iSbqugkCwCAV5fqJJuq6CQLAICtS3SSTWXEjAEAsKV0J9kugfIOAAA2YsZuETMGAMBGzNgxYsYAANiIGTtGzBgAAC9ixj5AzBgAAC9ixj5AzBgAABsxY8eIGQMAYCNm7BrlHQAAbMSM3SJmDACAjZixY5R4AACwUeJxjRIPAAA2SjxuUeIBAMBGiccxOskCAGCjk6xjdJIFAMCLTrI+QCdZAAC86CTrA3SSBQDARidZx4gZAwBgI2bsGuUdAABsxIzdImYMAICNmLFjxIwBALARM3aMmDEAAF7EjH2AmDEAAF7EjH2AmDEAADZixo4RMwYAwEbM2DXKOwAA2IgZu0XMGAAAGzFjxyjxAABgo8TjGiUeAABslHjcosQDAICNEo9jdJIFAMBGJ1nH6CQLAIAXnWR9gE6yAAB40UnWB+gkCwCAjU6yjhEzBgDARszYNco7AADYiBm7RcwYAAAbMWPHiBkDAGAjZuwYMWMAALyIGfsAMWMAALyIGfsAMWMAAGzEjB0jZgwAgI2YsWuUdwAAsBEzdouYMQAANmLGjlHiAQDARonHNUo8AADYKPG4RYkHAAAbJR7H6CQLAICNTrKO0UkWAAAvOsn6AJ1kAQDwopOsD9BJFgAAG51kHSNmDACA7TsTM164cKH69++v7OxsjRgxQnV1dS6n8yXKOwAA2L4LMeO//vWvmjFjhu677z5t2bJF559/vsaMGaP9+/e7mlIcMWMAAGzfiZjxI488oilTpuj666/XoEGDtHjxYvXo0UNPP/20qynFUeIBAMDW5Us8R48eVX19vcrLy7+cSHq6ysvLVVNTY41vaWlRNBr1PDoVJR4AAGxdvcTz8ccfq7W1VcFg0HM9GAwqHA5b46uqqhQIBOKP4uLizp0fJR4AACzfiRJPR8yaNUuRSCT+aGxs7NTXo5MsAAC2ZO6PGUl7pa8444wz1K1bNzU1NXmuNzU1KRQKWeOzsrKUlZW8utdFpfkqDGQrHDlCtQcA8J2XJikUyNZFpflJe00n76BkZmZq2LBhqq6ujl+LxWKqrq5WWVmZiyl5dEtP033jBkk6/i8FAIDvqhP74H3jBqlbEr8HxlmJZ8aMGVqyZImWL1+u9957TzfffLMOHz6s66+/3tWUPMYOLtSia4cqFPC+nXXyv5u8Ht2V16M7YxjDGMYwhjG+ev1EjQkFsrXo2qEaO7hQyeSkxCNJP//5z/Xf//5Xc+bMUTgc1gUXXKC1a9daH5x1aezgQv3voJDqGg5o/8Ej6ts7W8P69VH97k/jP594u4sxjGEMYxjDmK+Ocf36ifw7kvnOyQlpxpiU+5hFNBpVIBBQJBJRbm6u6+kAAIB26Mj+nRIpHgAA8N3CAQUAAPgOBxQAAOA7HFAAAIDvcEABAAC+wwEFAAD4DgcUAADgOxxQAACA73BAAQAAvuOs1f23caL5bTQadTwTAADQXif27fY0sU/JA8rBgwclScXFxY5nAgAAOurgwYMKBAJfOyYlv4snFotp79696t27t9LSOv8LjKLRqIqLi9XY2Mh3/3Qi1jk5WOfkYJ2Tg3VOjkStszFGBw8eVFFRkdLTv/5TJin5Dkp6errOPPPMpL9ubm4u/wEkAeucHKxzcrDOycE6J0ci1vmb3jk5gQ/JAgAA3+GAAgAAfIcDSjtkZWXpvvvuU1ZWluupdGmsc3KwzsnBOicH65wcLtY5JT8kCwAAujbeQQEAAL7DAQUAAPgOBxQAAOA7HFAAAIDvcED5BgsXLlT//v2VnZ2tESNGqK6uzvWUUlpVVZUuvPBC9e7dW3379tWECRO0c+dOz5gjR46osrJSBQUF6tWrlyZOnKimpiZHM+4aHnroIaWlpWn69Onxa6xzYnz00Ue69tprVVBQoJycHJ177rnavHlz/L4xRnPmzFFhYaFycnJUXl6uXbt2OZxx6mltbdXs2bNVWlqqnJwcnXXWWbr//vs93+fCOnfcG2+8oXHjxqmoqEhpaWlavXq153571vTAgQOqqKhQbm6u8vLydOONN+rQoUOJmaDBKa1cudJkZmaap59+2rz77rtmypQpJi8vzzQ1NbmeWsoaM2aMWbp0qdm+fbvZtm2bufzyy01JSYk5dOhQfMxNN91kiouLTXV1tdm8ebMZOXKkufjiix3OOrXV1dWZ/v37m/POO8/cfvvt8eus87d34MAB069fP3PdddeZ2tpa8+GHH5rXXnvNfPDBB/ExDz30kAkEAmb16tXm7bffNj/96U9NaWmp+fzzzx3OPLXMmzfPFBQUmFdeecU0NDSYF154wfTq1cv88Y9/jI9hnTvub3/7m7nnnnvMiy++aCSZVatWee63Z03Hjh1rzj//fLNx40bzj3/8w/zgBz8w11xzTULmxwHla1x00UWmsrIy/nNra6spKioyVVVVDmfVtezfv99IMhs2bDDGGNPc3Gy6d+9uXnjhhfiY9957z0gyNTU1rqaZsg4ePGgGDBhg1q1bZ370ox/FDyisc2Lcdddd5pJLLjnl/VgsZkKhkPnd734Xv9bc3GyysrLMc889l4wpdglXXHGFueGGGzzXrrrqKlNRUWGMYZ0T4eQDSnvWdMeOHUaS2bRpU3zMq6++atLS0sxHH330redEiecUjh49qvr6epWXl8evpaenq7y8XDU1NQ5n1rVEIhFJUn5+viSpvr5ex44d86z7wIEDVVJSwrqfhsrKSl1xxRWe9ZRY50R5+eWXNXz4cF199dXq27evhgwZoiVLlsTvNzQ0KBwOe9Y5EAhoxIgRrHMHXHzxxaqurtb7778vSXr77bf15ptv6rLLLpPEOneG9qxpTU2N8vLyNHz48PiY8vJypaenq7a29lvPISW/LDAZPv74Y7W2tioYDHquB4NB/etf/3I0q64lFotp+vTpGjVqlAYPHixJCofDyszMVF5enmdsMBhUOBx2MMvUtXLlSm3ZskWbNm2y7rHOifHhhx9q0aJFmjFjhn7zm99o06ZNuu2225SZmanJkyfH17Kt/4+wzu139913KxqNauDAgerWrZtaW1s1b948VVRUSBLr3Anas6bhcFh9+/b13M/IyFB+fn5C1p0DCpyprKzU9u3b9eabb7qeSpfT2Nio22+/XevWrVN2drbr6XRZsVhMw4cP14MPPihJGjJkiLZv367Fixdr8uTJjmfXdTz//PNasWKFnn32WZ1zzjnatm2bpk+frqKiIta5C6PEcwpnnHGGunXrZqUampqaFAqFHM2q65g2bZpeeeUVvf766zrzzDPj10OhkI4eParm5mbPeNa9Y+rr67V//34NHTpUGRkZysjI0IYNG7RgwQJlZGQoGAyyzglQWFioQYMGea6dffbZ2rNnjyTF15L/j3w7d955p+6++25NmjRJ5557rn7xi1/ojjvuUFVVlSTWuTO0Z01DoZD279/vuf/FF1/owIEDCVl3DiinkJmZqWHDhqm6ujp+LRaLqbq6WmVlZQ5nltqMMZo2bZpWrVql9evXq7S01HN/2LBh6t69u2fdd+7cqT179rDuHTB69Gi988472rZtW/wxfPhwVVRUxP+Zdf72Ro0aZcXk33//ffXr10+SVFpaqlAo5FnnaDSq2tpa1rkDPvvsM6Wne7erbt26KRaLSWKdO0N71rSsrEzNzc2qr6+Pj1m/fr1isZhGjBjx7SfxrT9m24WtXLnSZGVlmWXLlpkdO3aYqVOnmry8PBMOh11PLWXdfPPNJhAImL///e9m37598cdnn30WH3PTTTeZkpISs379erN582ZTVlZmysrKHM66a/hqiscY1jkR6urqTEZGhpk3b57ZtWuXWbFihenRo4f5y1/+Eh/z0EMPmby8PPPSSy+Zf/7zn2b8+PHEXzto8uTJ5nvf+148Zvziiy+aM844w8ycOTM+hnXuuIMHD5qtW7earVu3GknmkUceMVu3bjW7d+82xrRvTceOHWuGDBliamtrzZtvvmkGDBhAzDhZHn/8cVNSUmIyMzPNRRddZDZu3Oh6SilNUpuPpUuXxsd8/vnn5pZbbjF9+vQxPXr0MFdeeaXZt2+fu0l3EScfUFjnxFizZo0ZPHiwycrKMgMHDjRPPvmk534sFjOzZ882wWDQZGVlmdGjR5udO3c6mm1qikaj5vbbbzclJSUmOzvbfP/73zf33HOPaWlpiY9hnTvu9ddfb/P/x5MnTzbGtG9NP/nkE3PNNdeYXr16mdzcXHP99debgwcPJmR+acZ8pRUfAACAD/AZFAAA4DscUAAAgO9wQAEAAL7DAQUAAPgOBxQAAOA7HFAAAIDvcEABAAC+wwEFAAD4DgcUAADgOxxQAACA73BAAQAAvsMBBQAA+M7/B3qz1s2r8+4UAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "def neighbors(row: int, col: int) -> Iterator[Tuple[float, int, int]]:\n", + "def neighbors(elm: ELM, row: int, col: int) -> Iterator[Tuple[float, int, int]]:\n", " for row_offset in [-1, 0, 1]:\n", " for col_offset in [-1, 0, 1]:\n", + " nb_row = row + row_offset\n", + " nb_col = col + col_offset\n", + " \n", + " if nb_row < 0 or nb_row >= elm.n_rows or nb_col < 0 or nb_col >= elm.n_cols:\n", + " continue\n", " if row_offset == 0 and col_offset == 0:\n", " continue\n", - " else:\n", - " k = 1.0 if row_offset == 0 or col_offset == 0 else 2.0\n", - " yield (k, row + row_offset, col + col_offset)\n", + " \n", + " k = 1.0 if row_offset == 0 or col_offset == 0 else 2.0\n", + " yield (k, nb_row, nb_col)\n", "\n", - "def calculate_force(elm: ELM, row: int, col: int) -> numpy.ndarray:\n", + "def calculate_force(elm: ELM, displacement: numpy.ndarray, row: int, col: int) -> numpy.ndarray:\n", " K = elm.elastic_spring_constant\n", + " c = elm.bond_bending_constant\n", "\n", " total = numpy.array([0.0, 0.0])\n", - " for _, nb_row, nb_col in neighbors(row, col):\n", - " u_ij = elm.get_displacement(row, col) - elm.get_displacement(nb_row, nb_col)\n", + " for k, nb_row, nb_col in neighbors(elm, row, col):\n", + " u_ij = displacement[row, col, :] - displacement[nb_row, nb_col, :]\n", " x_ij = elm.get_lattice_position(row, col) - elm.get_lattice_position(nb_row, nb_col)\n", - " r_ij = elm.get_actual_position(row, col) - elm.get_actual_position(nb_row, nb_col)\n", + " r_ij = x_ij + u_ij\n", " \n", - " total -= K * numpy.dot(u_ij, x_ij) * (r_ij / numpy.linalg.norm(r_ij))\n", + " total -= K * k * numpy.dot(u_ij, x_ij) * (r_ij / numpy.linalg.norm(r_ij)) + (c * u_ij) / (numpy.linalg.norm(x_ij) ** 2)\n", "\n", " return total\n", "\n", "def gaussian(mu, sigma_sq, x):\n", " return numpy.exp(-(((x - mu) / numpy.sqrt(sigma_sq)) ** 2) / 2) / numpy.sqrt(2 * numpy.pi * sigma_sq)\n", "\n", - "e = ELM(20, 20)\n", - "for row in range(1, e.n_rows - 1):\n", - " e.set_displacement(row, 1, numpy.array([gaussian((e.n_rows - 1) / 2, 5.0, row) * -5, 0.0]))" + "e = ELM.make(100, 100)\n", + "n_plucked = 10\n", + "base = e.n_rows // 2 - n_plucked // 2\n", + "for row in range(n_plucked):\n", + " e.set_displacement(base + row, 0, numpy.array([\n", + " gaussian(n_plucked // 2, n_plucked / 3, row) * -5,\n", + " 0.0\n", + " ]))\n", + "\n", + "x, y = numpy.meshgrid(\n", + " numpy.linspace(0, e.n_rows - 1, e.n_rows, dtype=int),\n", + " numpy.linspace(0, e.n_cols - 1, e.n_cols, dtype=int)\n", + ")\n", + "\n", + "plt.scatter(\n", + " x + e.displacement[:, :, 0],\n", + " y + e.displacement[:, :, 1]\n", + ")\n", + "plt.show()" ] }, { @@ -158,26 +207,32 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 82, "id": "ff7adc30-191e-4b7c-b146-2a533463d6ea", "metadata": {}, "outputs": [], "source": [ "def velocity_verlet_step(elm: ELM, delta_t: float) -> ELM:\n", - " new_elm = ELM(elm.n_rows, elm.n_cols)\n", + " new_displacement = elm.displacement + elm.velocity * delta_t + 0.5 * elm.force * delta_t * delta_t\n", + " \n", + " new_force = numpy.zeros((elm.n_rows, elm.n_cols, 2))\n", + " for row in range(elm.n_rows):\n", + " for col in range(elm.n_cols):\n", + " new_force[row, col, :] = calculate_force(\n", + " elm,\n", + " new_displacement,\n", + " row,\n", + " col\n", + " )\n", + " \n", + " new_velocity = elm.velocity + 0.5 * (elm.force + new_force) * delta_t\n", "\n", - " new_elm.displacement = elm.displacement + elm.velocity * delta_t + 0.5 * elm.force * delta_t * delta_t\n", - " for row in range(1, elm.n_rows - 1):\n", - " for col in range(1, elm.n_cols - 1):\n", - " new_elm.set_force(row, col, calculate_force(new_elm, row, col))\n", - " new_elm.velocity = elm.velocity + 0.5 * (elm.force + new_elm.force) * delta_t\n", - "\n", - " return new_elm" + " return ELM(elm.n_rows, elm.n_cols, new_displacement, new_velocity, new_force)" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 85, "id": "bea964bf-6e94-4b72-a458-50c2964bbba1", "metadata": {}, "outputs": [], @@ -189,14 +244,14 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 86, "id": "421c0fd2-2dba-4fd5-89ff-0968227b08b9", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "32e0ea7e326e4b54badc0ba9b1e24e85", + "model_id": "0e9725fa01fe4ff6a520690877b4f5eb", "version_major": 2, "version_minor": 0 }, @@ -213,17 +268,12 @@ " None>" ] }, - "execution_count": 6, + "execution_count": 86, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "x, y = numpy.meshgrid(\n", - " numpy.linspace(0, e.n_rows - 1, e.n_rows, dtype=int),\n", - " numpy.linspace(0, e.n_cols - 1, e.n_cols, dtype=int)\n", - ")\n", - "\n", "def do_plot(i: int) -> None:\n", " fig, ax = plt.subplots(figsize=(10, 5), ncols=2)\n", " ax[0].scatter(\n", @@ -235,7 +285,8 @@ " y + es[i].displacement[:, :, 1],\n", " es[i].force[:, :, 0],\n", " es[i].force[:, :, 1],\n", - " scale = 10\n", + " angles = \"xy\",\n", + " scale = 100\n", " )\n", " ax[1].imshow(numpy.linalg.norm(es[i].displacement, axis=2), cmap=\"hot\")\n", " plt.show()\n", @@ -246,7 +297,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d419c328-a466-42fa-a8d9-14807d7100e4", + "id": "38e341e1-de34-41e7-b8e7-1322d27db18b", "metadata": {}, "outputs": [], "source": [] diff --git a/experiments/ELM_Debug.ipynb b/experiments/ELM_Debug.ipynb index b255e5b..5bef416 100644 --- a/experiments/ELM_Debug.ipynb +++ b/experiments/ELM_Debug.ipynb @@ -40,7 +40,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 1, "id": "05d4e22c-3d7d-4ac2-bf22-59b78aa40b84", "metadata": {}, "outputs": [], @@ -64,7 +64,7 @@ }, { "cell_type": "code", - "execution_count": 61, + "execution_count": 18, "id": "82bb4a02-1ecf-4c84-8b9d-b753dee75337", "metadata": {}, "outputs": [ @@ -78,10 +78,10 @@ { "data": { "text/plain": [ - "-0.20000000000000007" + "array([-0.14798801, -0.13453456])" ] }, - "execution_count": 61, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -91,27 +91,32 @@ " displacement: np.ndarray,\n", " i: Tuple[int, int],\n", " j: Tuple[int, int]\n", - ") -> float:\n", + ") -> np.ndarray:\n", " u_ij = displacement[i] - displacement[j]\n", " x_ij = np.array(i, dtype=float) - np.array(j, dtype=float)\n", " K_ij = np.linalg.norm(x_ij) ** 2\n", + " r_ij = x_ij + u_ij\n", "\n", " print(f\"{j=} {K_ij=} {u_ij=} {x_ij=}\")\n", "\n", - " return -K_ij * np.dot(u_ij, x_ij)\n", + " return -K_ij * np.dot(u_ij, x_ij) * r_ij / np.linalg.norm(r_ij)\n", "\n", "def calculate_force_bond(\n", " displacement: np.ndarray,\n", " i: Tuple[int, int],\n", " j: Tuple[int, int]\n", - ") -> float:\n", - " return 0.0\n", + ") -> np.ndarray:\n", + " u_ij = displacement[i] - displacement[j]\n", + " x_ij = np.array(i, dtype=float) - np.array(j, dtype=float)\n", + " c_ij = 1.0\n", + "\n", + " return (-c_ij * u_ij) / (np.linalg.norm(x_ij) ** 2)\n", "\n", "def calculate_force(\n", " displacement: np.ndarray,\n", " i: Tuple[int, int],\n", " j: Tuple[int, int]\n", - ") -> float:\n", + ") -> np.ndarray:\n", " return calculate_force_elastic(displacement, i, j) + calculate_force_bond(displacement, i, j)\n", "\n", "displacement = np.array([\n", @@ -124,19 +129,19 @@ }, { "cell_type": "code", - "execution_count": 71, + "execution_count": 19, "id": "f7a019c2-8d20-4e84-9e0c-ef788781a07a", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "89618a87c5514c6698d62e3344092fdf", + "model_id": "78145432b809422bb171f675b00d05a6", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "interactive(children=(FloatSlider(value=0.0, description='dp_x', max=1.0, min=-1.0, step=0.05), FloatSlider(va…" + "interactive(children=(FloatSlider(value=0.0, description='dp_x', max=1.5, min=-1.5, step=0.05), FloatSlider(va…" ] }, "metadata": {}, @@ -148,7 +153,7 @@ " None>" ] }, - "execution_count": 71, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } @@ -172,16 +177,14 @@ " ys.append(y + displacement[j][1])\n", "\n", " if j != i:\n", - " F_ij = calculate_force_elastic(displacement, i, j)\n", - " F_ij_vec = (\n", - " (i + displacement[i])\n", - " - (j + displacement[j])\n", - " ) * F_ij\n", - " print(f\"{F_ij_vec=}\")\n", - " elastic_xs.append(F_ij_vec[0])\n", - " elastic_ys.append(F_ij_vec[1])\n", + " F_ij = calculate_force(displacement, i, j)\n", + " print(f\"{F_ij=}\")\n", + " elastic_xs.append(F_ij[0])\n", + " elastic_ys.append(F_ij[1])\n", "\n", " pi = i + displacement[i]\n", + " plt.xlim(-1, 3)\n", + " plt.ylim(-1, 3)\n", " plt.scatter(xs, ys)\n", " plt.quiver(\n", " [pi[0]] * len(elastic_xs),\n", @@ -189,6 +192,7 @@ " elastic_xs,\n", " elastic_ys,\n", " scale = 5,\n", + " angles = \"xy\",\n", " pivot = \"tip\",\n", " width = 0.005\n", " )\n", @@ -198,14 +202,15 @@ " sum(elastic_ys),\n", " pivot = \"tip\",\n", " scale = 5,\n", + " angles = \"xy\",\n", " color = \"red\"\n", " )\n", " \n", "\n", "interact(\n", " do_plot,\n", - " dp_x = (-1.0, 1.0, 0.05),\n", - " dp_y = (-1.0, 1.0, 0.05)\n", + " dp_x = (-1.5, 1.5, 0.05),\n", + " dp_y = (-1.5, 1.5, 0.05)\n", ")" ] }, @@ -216,6 +221,14 @@ "metadata": {}, "outputs": [], "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f283012a-6b7c-48c3-8bdc-f45f5b458edc", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": {