50 / 60 : add exponential weight

50_ode/60_Duffing_Oscillator.ipynb

@@ -470,6 +470,7 @@
     "    t_end_sec:float=120.0\n",
     "    x0_0:np.ndarray=np.array((0.0, 0.0, 0.0, 0.0))\n",
     "    x0_1:np.ndarray=np.array((0.01, 0.01, 0.01, 0.01))\n",
+    "    lam:float = 0.1 # exponential weight\n",
     "\n",
     "    def __post_init__(self):\n",
     "        self.t_eval = np.arange(\n",
@@ -496,6 +497,10 @@
     "        )\n",
     "        t0, x0 = sol0.t, sol0.y\n",
     "\n",
+    "        assert len(self.t_eval) == sol0.y.shape[1], (\n",
+    "            f'len(self.t_eval) = {len(self.t_eval)} sol0.y.shape = {sol0.y.shape}'\n",
+    "        )\n",
+    "\n",
     "        sol1 = si.solve_ivp(\n",
     "            fun=Duffing_oscillator,\n",
     "            t_span=(self.t_start_sec, self.t_end_sec),\n",
@@ -513,9 +518,18 @@
     "\n",
     "    def l_exp(self, t0, x0, t1, x1):\n",
     "        assert x0.shape == x1.shape\n",
+    "        assert len(t0) == x0.shape[1], f'length different len(t0) = {len(t0)}, x0.shape = {x0.shape}'\n",
+    "        assert np.allclose(t0, t1), 'time different'\n",
+    "\n",
     "        n = len(x0)\n",
     "        dx = x0 - x1\n",
-    "        return (dx ** 2).sum(axis=1).mean()**0.5\n",
+    "        wdx = np.exp(self.lam*t0) * dx # exponential weight\n",
+    "\n",
+    "        result = (wdx ** 2).sum(axis=0).mean()**0.5\n",
+    "\n",
+    "        assert result > 0.0, f'result = {result}'\n",
+    "\n",
+    "        return result\n",
     "\n"
    ]
   },
@@ -558,9 +572,10 @@
     "    \n",
     "result = so.differential_evolution(\n",
     "    pcost.cost, bounds, popsize=100,\n",
-    "    workers=-1, updating='deferred',\n",
+    "    workers=-1,\n",
+    "    updating='deferred',\n",
     "    init=init,\n",
-    "    maxiter=max_iter\n",
+    "    maxiter=max_iter,\n",
     ")\n",
     "result\n",
     "\n"