{ "cells": [ { "cell_type": "markdown", "id": "2e90f0ee", "metadata": {}, "source": [ "# Where the smoothness breaks\n", "\n", "An integrator assumes the thing it is integrating is smooth across the step it is taking.\n", "Two things in this toolchain break that assumption on purpose: a **state event**, where a\n", "`when` condition crosses zero and the model changes discontinuously, and a **state-set\n", "switch**, where the set of variables being integrated is swapped out mid-run because the\n", "one in use went singular.\n", "\n", "Both stop the solver, change something, and start it again. This notebook looks at what\n", "that does to a trajectory — and at the numbers that say the engine did it right.\n", "\n", "The sliders need a live kernel, so on this page they are a snapshot. The file itself is\n", "attached here — [events.ipynb](events.ipynb) — but it reads its models from `examples/`, so it\n", "wants a checkout rather than a bare download:\n", "\n", "```bash\n", "git clone https://gitlab.com/jorgeecardona/pymodelica && cd pymodelica\n", "pip install \"modelica[sim]\" matplotlib ipywidgets\n", "jupyter lab docs/notebooks/events.ipynb\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "8e4414fc", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# The models come from `examples/` in the repository rather than from strings in here: a copy\n", "# pasted into a notebook is a copy that goes stale, and the reason this page is executed when the\n", "# site is built is that it then cannot show output from a model the project does not have.\n", "root = next(\n", " (p for p in [Path.cwd(), *Path.cwd().parents] if (p / \"examples\" / \"models\").is_dir()), None\n", ")\n", "if root is None:\n", " message = (\n", " \"This notebook reads its models from `examples/` in the pymodelica repository. \"\n", " \"Run it from a checkout:\\n\"\n", " \" git clone https://gitlab.com/jorgeecardona/pymodelica\\n\"\n", " \" cd pymodelica && jupyter lab docs/notebooks/\"\n", " )\n", " raise RuntimeError(message)\n", "\n", "SURFACE, INK, MUTED, LINE = \"#fcfcfb\", \"#0b0b0b\", \"#52514e\", \"#d4d3cd\"\n", "BLUE, ORANGE = \"#2a78d6\", \"#eb6834\"\n", "\n", "mpl.rcParams.update(\n", " {\n", " \"figure.figsize\": (8.4, 3.4),\n", " \"figure.dpi\": 110,\n", " \"figure.facecolor\": SURFACE,\n", " \"axes.facecolor\": SURFACE,\n", " \"savefig.facecolor\": SURFACE,\n", " \"axes.edgecolor\": LINE,\n", " \"axes.labelcolor\": MUTED,\n", " \"axes.spines.top\": False,\n", " \"axes.spines.right\": False,\n", " \"axes.grid\": True,\n", " \"axes.titlelocation\": \"left\",\n", " \"axes.titlepad\": 10,\n", " \"axes.titlesize\": 11,\n", " \"grid.color\": \"#e8e7e2\",\n", " \"grid.linewidth\": 0.8,\n", " \"lines.linewidth\": 2.0,\n", " \"legend.frameon\": False,\n", " \"text.color\": INK,\n", " \"xtick.color\": MUTED,\n", " \"ytick.color\": MUTED,\n", " \"font.size\": 10,\n", " }\n", ")\n", "%config InlineBackend.figure_formats = [\"svg\"]" ] }, { "cell_type": "markdown", "id": "0a42f852", "metadata": {}, "source": [ "## 1. One event per bounce\n", "\n", "`BouncingBall` is two equations and a `when`: `der(h) = v`, `der(v) = -g`, and\n", "`reinit(v, -e * pre(v))` when the height reaches zero. Every bounce is a state event, and\n", "the run reports them with the Modelica condition that fired." ] }, { "cell_type": "code", "execution_count": null, "id": "a75bb001", "metadata": {}, "outputs": [], "source": [ "from dataclasses import replace\n", "\n", "from modelica.ir import nodes as ir\n", "from modelica.lang import parse\n", "from modelica.sim import causalize, flatten, simulate\n", "\n", "ball_flat = flatten(\n", " [parse((root / \"examples\" / \"models\" / \"BouncingBall.mo\").read_text())], \"BouncingBall\"\n", ")\n", "ball = simulate(causalize(ball_flat), stop=4.0, points=1200)\n", "\n", "print(ball)\n", "for event in ball.events[:4]:\n", " print(f\" {event}\")\n", "print(f\" ... {len(ball.events)} in total\")" ] }, { "cell_type": "code", "execution_count": null, "id": "e6f47f2b", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(ball.time, ball[\"h\"], color=BLUE)\n", "for event in ball.events:\n", " ax.axvline(event.time, color=LINE, linewidth=1, zorder=0)\n", "ax.axhline(0, color=MUTED, linewidth=1)\n", "ax.set_title(f\"Height, and the {len(ball.events)} events that interrupted the integration\")\n", "ax.set_xlabel(\"time [s]\")\n", "ax.set_ylabel(\"h [m]\")\n", "ax.set_ylim(-0.02, 1.05)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b5986b60", "metadata": {}, "source": [ "## 2. What the restart costs\n", "\n", "Asked for no output grid, `simulate` returns the solver's own steps — the honest picture of\n", "where it had to work. Each bounce shows up twice: as a step cut short to land exactly on the\n", "crossing, and as a fresh start on the other side of it." ] }, { "cell_type": "code", "execution_count": null, "id": "0a5c82af", "metadata": {}, "outputs": [], "source": [ "raw = simulate(causalize(ball_flat), stop=4.0)\n", "gap = np.diff(raw.time)\n", "moving = gap > 0\n", "\n", "fig, ax = plt.subplots()\n", "ax.semilogy(raw.time[:-1][moving], gap[moving], color=BLUE, marker=\"o\", markersize=4, linewidth=1.2)\n", "for event in ball.events:\n", " ax.axvline(event.time, color=LINE, linewidth=1, zorder=0)\n", "ax.set_title(\n", " f\"Step size taken, {len(raw.time)} steps — and {int((~moving).sum())} restarts of zero length\"\n", ")\n", "ax.set_xlabel(\"time [s]\")\n", "ax.set_ylabel(\"Δt [s]\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "82585264", "metadata": {}, "source": [ "The steps collapse into each bounce and recover after it. The zero-length entries are the\n", "restarts themselves: the solver stops *at* the crossing and the next segment begins at the\n", "same instant with a different velocity, so the trajectory has two values at one time. That\n", "is not a numerical artifact to be smoothed over — it is what a discontinuity is.\n", "\n", "## 3. The bounces are a geometric series, and the run says so\n", "\n", "With restitution `e`, each bounce keeps a fraction `e` of the impact speed, so apex heights\n", "fall by `e²` and the intervals between bounces by `e`. Both are checkable against what came\n", "out of the integrator." ] }, { "cell_type": "code", "execution_count": null, "id": "f289baeb", "metadata": { "lines_to_next_cell": 1 }, "outputs": [], "source": [ "# Complete flights only: from the start (or a bounce) to the next bounce. The stretch after\n", "# the last event is still in the air when the run ends, so its apex is not one.\n", "from itertools import pairwise\n", "\n", "AIRBORNE = 1e-9 # below this the ball is resting, not flying, and its \"apex\" is arithmetic noise\n", "\n", "apex = []\n", "for start, end in pairwise([0.0, *(event.time for event in ball.events)]):\n", " window = (ball.time >= start) & (ball.time <= end)\n", " if window.any():\n", " apex.append(float(ball[\"h\"][window].max()))\n", "\n", "e = 0.8\n", "apex = np.array([h for h in apex if h > AIRBORNE])\n", "n = np.arange(len(apex))\n", "predicted = apex[0] * e ** (2 * n)\n", "\n", "fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 3.2))\n", "left.semilogy(n, predicted, color=ORANGE, linewidth=1.6, linestyle=(0, (4, 3)))\n", "left.semilogy(n, apex, color=BLUE, marker=\"o\", markersize=6, linewidth=0)\n", "left.annotate(\n", " \"h₀·e²ⁿ\",\n", " xy=(n[-3], predicted[-3]),\n", " xytext=(-38, -6),\n", " textcoords=\"offset points\",\n", " color=ORANGE,\n", " fontsize=9,\n", ")\n", "left.annotate(\n", " \"measured apex\",\n", " xy=(n[4], apex[4]),\n", " xytext=(10, 10),\n", " textcoords=\"offset points\",\n", " color=BLUE,\n", " fontsize=9,\n", ")\n", "left.set_title(\"Apex height per flight\")\n", "left.set_xlabel(\"flight\")\n", "left.set_ylabel(\"apex [m]\")\n", "left.set_xticks(n[::3]) # flights are counted, so the ticks are whole numbers\n", "\n", "intervals = np.diff([event.time for event in ball.events])\n", "ratio = intervals[1:] / intervals[:-1]\n", "right.plot(np.arange(1, len(ratio) + 1), ratio, color=BLUE, marker=\"o\", markersize=6)\n", "right.axhline(e, color=ORANGE, linewidth=1.6, linestyle=(0, (4, 3)))\n", "right.annotate(\n", " f\"e = {e}\", xy=(1, e), xytext=(0, 10), textcoords=\"offset points\", color=ORANGE, fontsize=9\n", ")\n", "right.set_title(\"Ratio of successive flight times\")\n", "right.set_xlabel(\"bounce\")\n", "right.set_ylabel(\"tₙ₊₁ / tₙ\")\n", "right.set_ylim(0.74, 0.86)\n", "right.set_xticks(np.arange(1, len(ratio) + 1, 3))\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "grid = float(np.diff(ball.time).max())\n", "print(f\"mean interval ratio {ratio.mean():.9f} against e = {e}\")\n", "print(f\"worst apex error {np.abs(apex - predicted).max():.2e} m\")\n", "print(\n", " f\"grid can miss a peak by {9.81 / 2 * (grid / 2) ** 2:.2e} m (samples {1e3 * grid:.2f} ms apart)\"\n", ")" ] }, { "cell_type": "markdown", "id": "ad3174de", "metadata": {}, "source": [ "The two numbers are measuring different things, and the difference is the point.\n", "\n", "The interval ratio is the **event locator** being graded: those times came from the solver\n", "stopping at the crossings, not from anything analytic. Nine digits of agreement means each\n", "zero-crossing was found to solver tolerance, not merely bracketed.\n", "\n", "The apex error is not the solver at all — it is the output grid. An apex read off sampled\n", "points is the highest *sample*, and the true peak sits up to half a sample away; for a\n", "parabola that is `g/2·(Δt/2)²`, which is the third number, and it matches. Asking for output\n", "at the apexes themselves would remove it, since requested times are filled from the\n", "integrator's own interpolant rather than replacing its steps.\n", "\n", "## 4. Restitution on a slider\n", "\n", "Dropping `e` makes the bounces die out sooner; raising it toward 1 pushes the run into the\n", "Zeno regime, where infinitely many events pile up before a finite time. The engine gives up\n", "after a hundred consecutive zero-length segments and says so rather than hanging." ] }, { "cell_type": "code", "execution_count": null, "id": "931c71fd", "metadata": {}, "outputs": [], "source": [ "def bounce(e: float, h0: float, stop: float) -> None:\n", " \"\"\"Simulate the ball with that restitution and drop height, and plot it.\"\"\"\n", " values = {\"e\": e}\n", " starts = {\"h\": h0}\n", " tuned = replace(\n", " ball_flat,\n", " variables=tuple(\n", " replace(\n", " v,\n", " binding=ir.RealLiteral(values[v.name]) if v.name in values else v.binding,\n", " start=ir.RealLiteral(starts[v.name]) if v.name in starts else v.start,\n", " )\n", " for v in ball_flat.variables\n", " ),\n", " )\n", " result = simulate(causalize(tuned), stop=stop, points=1200)\n", "\n", " fig, ax = plt.subplots()\n", " ax.plot(result.time, result[\"h\"], color=BLUE)\n", " for event in result.events:\n", " ax.axvline(event.time, color=LINE, linewidth=1, zorder=0)\n", " ax.axhline(0, color=MUTED, linewidth=1)\n", " ax.set_title(\n", " f\"e = {e:.2f}, dropped from {h0:.2f} m — {len(result.events)} bounces in {stop:.0f} s\"\n", " )\n", " ax.set_xlabel(\"time [s]\")\n", " ax.set_ylabel(\"h [m]\")\n", " ax.set_ylim(-0.02, h0 * 1.06)\n", " plt.show()\n", "\n", "\n", "bounce(e=0.8, h0=1.0, stop=4.0)" ] }, { "cell_type": "code", "execution_count": null, "id": "c9f943cf", "metadata": {}, "outputs": [], "source": [ "import ipywidgets as widgets\n", "\n", "_ = widgets.interact(\n", " bounce,\n", " e=widgets.FloatSlider(\n", " value=0.8, min=0.1, max=0.95, step=0.05, description=\"e\", continuous_update=False\n", " ),\n", " h0=widgets.FloatSlider(\n", " value=1.0, min=0.2, max=3.0, step=0.1, description=\"h start\", continuous_update=False\n", " ),\n", " stop=widgets.FloatSlider(\n", " value=4.0, min=1.0, max=12.0, step=1.0, description=\"stop\", continuous_update=False\n", " ),\n", ")" ] }, { "cell_type": "markdown", "id": "cc5adc50", "metadata": {}, "source": [ "## 5. The other kind of interruption\n", "\n", "The pendulum has no `when` in it and still gets interrupted. `x² + y² = L²` makes it an\n", "index-3 DAE, and index reduction leaves two admissible state sets — `(y, vy)` and\n", "`(x, vx)` — neither of which works all the way around the swing. The engine switches\n", "between them while the run is going, and each switch is an event with the set it moved to\n", "as its cause." ] }, { "cell_type": "code", "execution_count": null, "id": "28b6e316", "metadata": {}, "outputs": [], "source": [ "pendulum_flat = flatten(\n", " [parse((root / \"examples\" / \"models\" / \"Pendulum.mo\").read_text())], \"Pendulum\"\n", ")\n", "pendulum = causalize(pendulum_flat)\n", "swing = simulate(pendulum, stop=4.0, points=1600)\n", "\n", "print(f\"state sets: {pendulum.candidates}\")\n", "for event in swing.events:\n", " print(f\" {event}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "9121b0e8", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "ax.plot(swing.time, swing[\"x\"], color=BLUE)\n", "ax.plot(swing.time, swing[\"y\"], color=ORANGE)\n", "ax.annotate(\n", " \"x\",\n", " xy=(swing.time[-1], swing[\"x\"][-1]),\n", " xytext=(6, 0),\n", " textcoords=\"offset points\",\n", " color=BLUE,\n", " fontsize=10,\n", " va=\"center\",\n", ")\n", "ax.annotate(\n", " \"y\",\n", " xy=(swing.time[-1], swing[\"y\"][-1]),\n", " xytext=(6, 0),\n", " textcoords=\"offset points\",\n", " color=ORANGE,\n", " fontsize=10,\n", " va=\"center\",\n", ")\n", "for event in swing.events:\n", " ax.axvline(event.time, color=LINE, linewidth=1, zorder=0)\n", "ax.set_title(f\"The swing, with the {len(swing.events)} state-set switches marked\")\n", "ax.set_xlabel(\"time [s]\")\n", "ax.set_ylabel(\"position [m]\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3162cb15", "metadata": {}, "source": [ "## 6. Why it has to switch\n", "\n", "Each state set causalizes to its own blocks, and each carries a **pivot** — the expression\n", "the solve divides by, kept symbolically so it can be evaluated along the trajectory. For\n", "`(y, vy)` it is `4x²`; for `(x, vx)` it is `4y²`. Where a pivot passes through zero its set\n", "is unusable, and structural analysis cannot see that coming: it is a statement about\n", "numbers, not about incidence.\n", "\n", "Evaluating both along the run puts the whole argument on one chart. They cross where\n", "`|x| = |y|`, which is 45 degrees, which is the answer the literature gives." ] }, { "cell_type": "code", "execution_count": null, "id": "d5cea705", "metadata": {}, "outputs": [], "source": [ "from modelica.ir.evaluate import evaluate\n", "from modelica.lang import unparse_expr\n", "\n", "pivots = []\n", "for mode in pendulum.modes:\n", " print(f\"{mode.states} pivot = {unparse_expr(mode.pivot)}\")\n", " pivots.append(\n", " np.array([abs(float(evaluate(mode.pivot, swing.at(float(t))))) for t in swing.time])\n", " )\n", "\n", "fig, ax = plt.subplots()\n", "for series, mode, colour in zip(pivots, pendulum.modes, (BLUE, ORANGE), strict=True):\n", " ax.semilogy(\n", " swing.time, np.maximum(series, 1e-12), color=colour, label=f\"({', '.join(mode.states)})\"\n", " )\n", "for event in swing.events:\n", " ax.axvline(event.time, color=LINE, linewidth=1, zorder=0)\n", "ax.set_title(\"How well conditioned each state set is, along the trajectory\")\n", "ax.set_xlabel(\"time [s]\")\n", "ax.set_ylabel(\"|pivot|\")\n", "ax.legend(loc=\"lower left\", ncols=2, labelcolor=MUTED)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "60902599", "metadata": {}, "source": [ "The switches land near the crossings, which is the engine choosing the better-conditioned\n", "set — not exactly on them, because a set is only abandoned once it actually stops working,\n", "and \"stops working\" is something the linear algebra reports rather than something a\n", "threshold predicts.\n", "\n", "## 7. The point of all of it\n", "\n", "A constraint that has been differentiated twice is no longer in the equations being\n", "integrated. Nothing in the right-hand side says the rod has a fixed length any more — so if\n", "index reduction and state selection were subtly wrong, the mass would drift off the circle\n", "while satisfying every equation the solver was given.\n", "\n", "It does not:" ] }, { "cell_type": "code", "execution_count": null, "id": "def1bc1e", "metadata": { "lines_to_next_cell": 1 }, "outputs": [], "source": [ "EPS = np.finfo(float).eps\n", "\n", "drift = np.abs(np.hypot(swing[\"x\"], swing[\"y\"]) - 1.0)\n", "worst = np.maximum.accumulate(drift) / EPS # worst so far, in machine epsilons\n", "\n", "fig, ax = plt.subplots()\n", "ax.step(swing.time, worst, where=\"post\", color=BLUE)\n", "ax.axhline(1.0, color=ORANGE, linewidth=1.6, linestyle=(0, (4, 3)))\n", "ax.annotate(\n", " \"one machine epsilon\",\n", " xy=(0.15, 1.0),\n", " xytext=(0, 8),\n", " textcoords=\"offset points\",\n", " color=ORANGE,\n", " fontsize=9,\n", ")\n", "for event in swing.events:\n", " ax.axvline(event.time, color=LINE, linewidth=1, zorder=0)\n", "ax.set_title(\"Worst drift off the rod so far, over the whole run\")\n", "ax.set_ylim(0, 1.6)\n", "ax.set_xlabel(\"time [s]\")\n", "ax.set_ylabel(\"| √(x²+y²) − L | [ε]\")\n", "plt.show()\n", "\n", "print(f\"worst drift over the run: {drift.max():.3e} m = {drift.max() / EPS:.2f} machine epsilons\")" ] }, { "cell_type": "markdown", "id": "1ff39ca0", "metadata": {}, "source": [ "## 8. Amplitude on a slider\n", "\n", "The number of switches is not a property of the model, it is a property of the run. A small\n", "swing stays where `(y, vy)` is fine and never switches at all; a large one crosses 45\n", "degrees four times per period." ] }, { "cell_type": "code", "execution_count": null, "id": "edcc5ba0", "metadata": {}, "outputs": [], "source": [ "def swing_from(degrees: float, stop: float) -> None:\n", " \"\"\"Release the pendulum from that angle off vertical and count the switches.\"\"\"\n", " theta = np.radians(degrees)\n", " starts = {\"x\": float(np.sin(theta)), \"y\": float(-np.cos(theta))}\n", " tuned = replace(\n", " pendulum_flat,\n", " variables=tuple(\n", " replace(v, start=ir.RealLiteral(starts[v.name]) if v.name in starts else v.start)\n", " for v in pendulum_flat.variables\n", " ),\n", " )\n", " result = simulate(causalize(tuned), stop=stop, points=1200)\n", "\n", " fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 3.2))\n", " left.plot(result.time, result[\"x\"], color=BLUE)\n", " left.plot(result.time, result[\"y\"], color=ORANGE)\n", " for event in result.events:\n", " left.axvline(event.time, color=LINE, linewidth=1, zorder=0)\n", " left.set_title(f\"{degrees:.0f}° off vertical — {len(result.events)} switches\")\n", " left.set_xlabel(\"time [s]\")\n", " left.set_ylabel(\"position [m]\")\n", " left.set_ylim(-1.15, 1.15)\n", "\n", " right.plot(result[\"x\"], result[\"y\"], color=BLUE, linewidth=1.4)\n", " right.plot([0.0], [0.0], marker=\"o\", markersize=8, color=MUTED)\n", " right.set_title(\"Path of the mass\")\n", " right.set_xlabel(\"x [m]\")\n", " right.set_ylabel(\"y [m]\")\n", " right.set_xlim(-1.15, 1.15)\n", " right.set_ylim(-1.15, 1.15)\n", " right.set_aspect(\"equal\")\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "swing_from(degrees=90.0, stop=4.0)" ] }, { "cell_type": "code", "execution_count": null, "id": "aafb9e8a", "metadata": {}, "outputs": [], "source": [ "_ = widgets.interact(\n", " swing_from,\n", " degrees=widgets.FloatSlider(\n", " value=90.0, min=5.0, max=175.0, step=5.0, description=\"angle °\", continuous_update=False\n", " ),\n", " stop=widgets.FloatSlider(\n", " value=4.0, min=1.0, max=10.0, step=1.0, description=\"stop\", continuous_update=False\n", " ),\n", ")" ] }, { "cell_type": "markdown", "id": "1eacca92", "metadata": {}, "source": [ "Both interruptions come back to the same design point: the discrete half of a simulation is\n", "*finite*, so it can be enumerated rather than sampled. Which is\n", "[the next notebook](../verifying/) — the same switching logic, handed to a model checker\n", "instead of an integrator." ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "main_language": "python", "notebook_metadata_filter": "-all" } }, "nbformat": 4, "nbformat_minor": 5 }