2e:null 2f:["$","title",null,{"children":"Phase Portraits | Jo Wood | Observable"}] 34:T5af,# Phase Portraits _Originally written in January 2022, this was updated in September 2026 to use Observable 2.0._ Phase portraits are a special case of [connected scatterplots](http://steveharoz.com/research/connected_scatterplot/) that are useful for visualizing dynamic systems. Particularly so when it is easier to say how a system changes from one moment to the next than it is to predict the state of the system at an arbitrary point in time. They provide insight not simply into the state of a system, but the underlying processes that determine that state. An assumption made in what follows is that we have an interest in understanding systems where change over time is in part a property of the state of the system. For example we might be interested in the temperature of something, where the cooling rate is dependent on its absolute temperature. In other words, systems with some kind of feedback (_'dynamical systems'_). While common in physics, I think they are underused in data visualization ([Danny Dorling's use of phase portraits](https://www.dannydorling.org/books/SLOWDOWN/Illustrations.html#_) being a notable exception). They should be particularly useful in data visualization because they are more than descriptors of some state, but rather a characterisation of _process_, and as such closer to a representation of knowledge rather than data. To understand their benefit, let's turn first to a simple example from physics...35:T646,{ restart; const [w, h] = [600, 500]; const rShort = 50; const p0 = [w / 2, 100]; const p0a = [w / 2, 100 + rShort]; const deltaT = 0.05; const equals = (a, b) => Math.abs(a - b) <= 0.005; const accel = (theta, thetaDot) => -mu * thetaDot - (g / (0.02 * l)) * Math.sin(theta); let t = 0; let theta = radians(theta0); let thetaDot = thetaDot0; let thetaDoubleDot = accel(theta, thetaDot); let stopped = false; invalidation.then(() => { stopped = true; }); const r = new Renderer(w, h).textSize(14); const draw = () => { const p1 = [p0[0] + l * Math.sin(theta), p0[1] + l * Math.cos(theta)]; const p2 = [ p0[0] + rShort * Math.sin(theta / 2.5), p0[1] + rShort * Math.cos(theta / 2.5) ]; r.clear() .stroke("black") .fill("rgb(220,150,120)") .line([p0, p1]) .stroke("rgba(0,0,0,0.5)") .line([p0, p0a]) .ellipse(...p1, 15) .stroke() .fill("rgba(0,0,0,0.5)") .textAlign("left") .text(`time t=${Math.round(t)}`, 2, 20) .textAlign("middle") .fill( `rgba(0,0,0,${Math.abs(theta) < 0.5 ? Math.abs(Math.sin(theta)) : 0.5})` ) .text("θ", ...p2); }; draw(); display(r.render()); while ( restart > 0 && !stopped && (!equals(thetaDoubleDot, 0) || !equals(thetaDot, 0)) ) { await new Promise((resolve) => setTimeout(resolve, 30)); if (stopped) { break; } t += deltaT; thetaDoubleDot = accel(theta, thetaDot); thetaDot += thetaDoubleDot * deltaT; theta += thetaDot * deltaT; draw(); } }36:T824,Its position, velocity and acceleration all vary over time. Position is determined by the angle ${tex `\theta`} with the vertical and its length ${tex `l`} as well as the acceleration due to gravity. But given the length of any given pendulum is fixed, as is gravitational acceleration, what we really need to know is the angle ${tex `\theta`} at time ${tex `t`}, or more compactly ${tex `\theta(t)`}. Assuming there is some air resistance that is proportional to the velocity of the pendulum, for any given angle ${tex `\theta`} we can determine the acceleration of the pendulum with some simple geometry: ${tex `a = \dfrac{d^2 \theta}{dt^2} = -\mu \dfrac{d \theta}{dt} - \dfrac{g}{l} sin(\theta)`} where ${tex `g`} is the acceleration due to gravity (9.81 metres per second per second on Earth) and ${tex `\mu`} is a coefficient of air resistance (0 = no resistance). For compactness, we can represent that second derivative of ${tex `\theta`} with respect to time (acceleration) as ${tex `\ddot{\theta}(t)`}, and the first derivative with respect to time (velocity) as ${tex `\dot{\theta}(t)`}: ${tex `\ddot{\theta}(t) = - \mu \dot{\theta}(t) - \dfrac{g}{l} sin(\theta(t))`} But what we have above is a description of how its _acceleration_ changes with time, not what we ultimately need to know: how _position_ changes with time. To find its position we could solve that [ordinary differential equation](https://en.wikipedia.org/wiki/Ordinary_differential_equation) and but even for an apparently simple case like this, it is surprisingly complex to solve algebraically. Instead, consider a visualization of the space of all possible states of the pendulum. From the equation above, we know its state can be described fully by its position ${tex `\theta`} and velocity ${tex `\dot{\theta}`}, from which its acceleration may be determined. So let's see that space where position is plotted on the horizontal axis and velocity on the vertical axis. At each point we represent the pendulum's velocity and acceleration as a vector and visualize the resulting vector field.37:T422,This map of 'phase space' shows what the pendulum will do next given its position and velocity. Following the vector at any given point will generally update the position and velocity to new values. Note the static positions at (-360,0), (-180,0), (0,0), (180,0) and (360,0) where both acceleration and velocity is 0. These correspond to the pendulum resting vertically either at the bottom or top of its swing. Those at (-360,0), (0,0) and (360,0), where the pendulum is resting at the bottom of its swing, are [attractors](https://en.wikipedia.org/wiki/Attractor) in that a small nudge in either velocity or position will soon return back to the same resting point. In contrast, those at (-180,0) and (180,0) are 'repulsers' at the top of their swing. Just a tiny nudge will move it away as it falls from its balance point. We can see the difference between attractors and repulsers because the former is surrounded by vectors that spiral towards them and the latter by vectors that point away (try turning up the air resistance to make this more obvious).38:T496,What does this phase portrait tell us? And does it tell us anything more than the plot of stocks over time? * The flat top to the shape tells us there are periods when stocks are high, but whatever is governing rate of change is itself changing rapidly ("acceleration"). We don't see the same effect when stocks are low. This could be discerned from the shape curve around the minima and maxima of the stocks vs time chart, but is not so apparent. * The spiral tells us the amplitude of oscillation is decreasing over time, which is arguably more intuitive in the stocks vs time chart. In this model, that attenuation is introduced by progressively damping the fish-population update in a modified Runge–Kutta step. * The left-right asymmetry tells us that whatever is governing the decline in stocks is more influential than whatever is governing their recovery. And that as oscillation decreases, that asymmetry is more pronounced. None of this is particularly obvious from the stocks vs time chart. * There is an 'attractor' at around (0,750) suggesting a possible stable population. This is also obvious from the stocks vs time chart (and arguably more intuitive).39:T549,{ // Lotka–Volterra predator-prey equations: // const alpha = 1.1; // const beta = 0.4; // const delta = 0.4; // const gamma = 0.1; const preyFn = (prey, pred) => selLV.alpha * prey - selLV.beta * prey * pred; const predFn = (prey, pred) => selLV.gamma * prey * pred - selLV.delta * pred; const [preyPop, predPop] = [10, 10]; // Solve differential equation over 100t in timesteps of 0.2 const [t, dt] = [50, 1 / 12]; const pops = dampedRungeKutta4(preyPop, predPop, preyFn, predFn, dt, t / dt); const prey = pops.map((pop, i) => [i * dt, pop[0]]); const pred = pops.map((pop, i) => [i * dt, pop[1]]); const maxPop = Math.max(...pops.flat()); const [w, h] = [600, 300]; const xSc = d3.scaleLinear().domain([0, t]).range([0, w]); const ySc = d3.scaleLinear().domain([0, maxPop]).range([0, h]); const r = new Renderer(w, h) .setScales(xSc, ySc) .flipY() .stroke("firebrick") .line(prey); if (showControl) { r.stroke("steelblue") .line(pred) .setScales( d3.scaleLinear().domain([0, 1]).range([0, w]), d3.scaleLinear().domain([0, 1]).range([0, h]) ) .stroke() .textAlign("right") .fill("steelblue") .text("Predator population", 0.9, 0.95) .fill("firebrick") .text("Fish population", 0.9, 0.88); } display(r.render()); }3a:T430,function addRelief(r) { const [nRows, nCols] = [pGrid.length, pGrid[0].length]; r.push().stroke(); const aSc = d3.scaleLinear().domain([0, 8]).range([0, 1]); const sinAlt = Math.sin(radians(35)); const cosAlt = Math.cos(radians(35)); const azim = radians(315); for (let row = 0; row < nRows; row++) { for (let col = 0; col < nCols; col++) { const [vel, accel] = pGrid[row][col]; const slope = halfPi - Math.atan(Math.sqrt(vel * vel + accel * accel)); const aspect = Math.atan2(accel, vel) + halfPi; const reflect = (255 * (Math.max( 0, sinAlt * Math.sin(slope) + cosAlt * Math.cos(slope) * Math.cos(azim - aspect) ) + 1)) / 2; r.fill(`rgba(${reflect},${reflect},${0.2 + reflect},0.15)`).polygon([ [pxScale(col - 0.5), pyScale(row - 0.5)], [pxScale(col - 0.5), pyScale(row + 0.5)], [pxScale(col + 0.5), pyScale(row + 0.5)], [pxScale(col + 0.5), pyScale(row - 0.5)] ]); } } return r.pop(); }3b:[] 32:{"id":"5db0e20023c5a0f1","slug":"phaseportrait","title":"Phase Portraits","thumbnail":"06e53bdde821345f45e7ae3d695236d16d7f55d56a916fac85c1d7b783273dd4","access_level":"public","create_time":"$D2021-12-30T17:35:33.300Z","update_time":"$D2026-09-13T12:00:13.201Z","view_time":null,"first_public_time":"$D2022-01-12T14:01:57.859Z","has_unpublished_edits":null,"fork":false,"stars":26,"starred":false,"owner":{"type":"team","id":"9461f58b209074d5","login":"jwolondon","name":"Jo Wood","avatar_url":"https://avatars-next.observableusercontent.com/27f72cecd1ac20853a533421fab228ce5f611cbc0986c1e420f7d3b1eb49ef50","tier":"starter_2024"},"version":2643,"publish_version":2643,"role":null,"authored":false,"authors":[{"type":"individual","id":"e2432cc13dc447d3","login":"jwolondon","name":"Jo Wood","avatar_url":"https://avatars-next.observableusercontent.com/27f72cecd1ac20853a533421fab228ce5f611cbc0986c1e420f7d3b1eb49ef50","tier":null}],"members":[{"type":"individual","id":"e2432cc13dc447d3","login":"jwolondon","name":"Jo Wood","avatar_url":"https://avatars-next.observableusercontent.com/27f72cecd1ac20853a533421fab228ce5f611cbc0986c1e420f7d3b1eb49ef50","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":6,"collection":{"id":"1986e99b3893465a","slug":"excursions","title":"Excursions","owner":{"id":"9461f58b209074d5","login":"jwolondon","name":"Jo Wood"}},"type":"notebook","body":{"title":"Phase Portraits","cells":[{"id":0,"value":"$34","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":309,"value":"## 1. Known process, unknown state\n\nConsider a swinging pendulum where we wish to predict where it will be at time ${tex `t`}.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":2428,"value":"const restart = view(\n Inputs.button(\"Swing the pendulum\", {\n value: 0,\n reduce: (count) => count + 1\n })\n);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":17,"value":"$35","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":297,"value":"const lInput = Inputs.range([75, 275], {\n label: \"Pendulum length\",\n step: 5,\n value: 200\n});\n\nconst l = view(lInput);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":362,"value":"const muInput = Inputs.range([0, 2], {\n label: \"Air resistance\",\n step: 0.02,\n value: 0.1\n});\n\nconst mu = view(muInput);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":74,"value":"$36","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":693,"value":"{\n const r = pendulumPhaseSpace(900, 320);\n addVectorField(r);\n display(r.render());\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2294,"value":"_An aside: Should we plot position on the x-axis and velocity on the y or the other way around? It is common to place independent variables on the x-axis and dependent ones on the y. However, by definition, we have feedback between the two variables (position determines velocity and velocity determines position), so there appears to be no natural choice here. In this case, there is more interesting variation in position so we will devote more space to this in a landscape orientation._","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":696,"value":"$37","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":718,"value":"Unlike solving the differential equation, we can quite easily compute the evolving position and velocity of the pendulum by updating both in response to the (acceleration, velocity) value of the vector in phase space. If we do so in small increments we can visualise the trajectory of the moving pendulum in phase space (where its position in 'real' space is simply determined from the x-coordinate in phase space):","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":2270,"value":"Inputs.bind(\n Inputs.range([75, 275], {label: \"Pendulum length\", step: 5}),\n lInput\n)","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":553,"value":"const thetaDot0 = view(\n Inputs.range([-4, 4], {\n label: \"Initial velocity\",\n step: 0.01,\n value: 0\n })\n);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":530,"value":"const theta0 = view(\n Inputs.range([-180, 180], {\n label: \"Initial position\",\n step: 1,\n value: 60\n })\n);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":985,"value":"Inputs.bind(\n Inputs.range([0, 2], {label: \"Air resistance\", step: 0.02}),\n muInput\n)","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":643,"value":"{\n const r = pendulumPhaseSpace(900, 320).strokeWidth(0.7);\n // addRelief(r);\n addVectorField(r, true);\n r.strokeWidth(2);\n addTrajectory(r, theta0, thetaDot0, 150, \"black\");\n display(r.render());\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2313,"value":"To help visualize the phase space we can colour vectors according to whether they are curving left or right in the direction of travel. This allows us to see 'cyclonic' and 'anticyclonic' rotations more easily, as well as provide a greater visual distinction between point attractors and repulsers.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1013,"value":"We can further enhance the depiction of phase space by creating a 'phase portrait' of many different starting positions and leaving a trace of the pendulum over time in phase space for each of them. I have also added a subtle shading to the surface based on the magnitude of the local vector at regular points on the surface. This creates a meaningful visual summary of the dynamics of the system.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":726,"value":"{\n const random = d3.randomLcg(\"appleSeed\");\n const rand = (min = 0, max = 1) => random() * (max - min) + min;\n\n const r = pendulumPhaseSpace(900, 320).strokeWidth(0.3);\n addRelief(r);\n for (let i = 0; i < 2000; i++) {\n addTrajectory(r, rand(-560, 560), rand(-6.5, 6.5), 0.4);\n }\n if (showPendulum) {\n r.strokeWidth(1.5);\n addTrajectory(r, theta0, thetaDot0, 150, \"black\");\n }\n display(r.render());\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2280,"value":"const showPendulum = view(\n Inputs.toggle({label: \"Show pendulum\", value: false})\n);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1018,"value":"In summary:\n\n* For dynamical systems it can often be easier to describe the how a system changes than it is to describe the state of a system.\n* If we know how a system changes over time, we can visualise its phase space by plotting its state (e.g. angular position of a pendulum) vs rate of change of state (e.g. velocity of the pendulum).\n* The phase portrait shows many possible states of the system and trajectories showing how states change over time.\n* We can model the state at time ${tex `t`} numerically by incrementing ${tex `t`} in small steps and iteratively updating its state and derivatives (commonly using [Runge-Kutta methods](https://en.wikipedia.org/wiki/Runge–Kutta_methods)).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":315,"value":"## 2. Known state, unknown processes\n\nThe example of the pendulum described a system where we could model its dynamics precisely (we had a differential equation describing the acceleration of the pendulum with respect to time), but we did not know the state of the system until we modelled it numerically.\n\nIn contrast, most data visualizations are descriptions of some state (albeit imperfectly through the lens of the data we have available), but we don't know the underlying dynamics. Can phase plots offer us any insight into these dynamics?","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1033,"value":"Consider a simple example of a population of a particular species of fish in an ocean. \n\n### BigBucksFishingCo\n\nThe BigBucksFishingCo has generously given you temporary access to its long-term oceanographic survey where it has been sampling the population every month for 50 years. Luckily for us this gives us plenty of reliable data on fish stocks:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1805,"value":"{\n const maxPop = Math.max(...fish.map((f) => f[1]));\n const [w, h, hPad, vPad] = [600, 300, 30, 20];\n\n const xSc = d3\n .scaleTime()\n .domain([new Date(1972, 0, 1), new Date(2021, 11, 31)])\n .range([hPad, w]);\n const ySc = d3\n .scaleLinear()\n .domain([0, maxPop])\n .range([vPad, h - vPad]);\n\n display(\n new Renderer(w, h)\n .setScales(xSc, ySc)\n .flipY()\n .textSize(8)\n .axisBottom(xSc, ySc)\n .axisLeft(xSc, ySc)\n .strokeWidth(2)\n .stroke(\"firebrick\")\n .line(fish)\n .render()\n );\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1860,"value":"It looks like we have a regular cyclical pattern in our fish population. But why might this be? Let's try plotting it as a phase portrait where population is plotted on the y-axis and change in population on the x-axis:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1866,"value":"{\n const change = (values, i) => (values[i + 1] - values[i - 1]) / 2;\n\n const [w, h, border] = [500, 400, 20];\n\n const pop = fish.map((f) => f[1]);\n const maxPop = Math.max(...pop);\n\n const deltaPop = pop\n .map((p, i) => [change(pop, i), p])\n .slice(1, pop.length - 1);\n\n const xSc = d3\n .scaleLinear()\n .domain(d3.extent(deltaPop.map((dp) => dp[0])))\n .range([border, w - border]);\n\n const ySc = d3\n .scaleLinear()\n .domain([0, maxPop])\n .range([border, h - border]);\n\n display(\n new Renderer(w, h)\n .fill()\n .setScales(xSc, ySc)\n .textSize(8)\n .flipY()\n .axisBottom(xSc, ySc)\n .axisLeft(xSc, ySc, [0, maxPop], 0)\n .stroke(\"firebrick\")\n .bezierLine(deltaPop)\n .render()\n );\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2320,"value":"$38","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":2343,"value":"### NeedABiggerBoat Fisheries\n\nUnfortunately, we tend not to have the resources of BigBucksFishingCo when gathering data. We are much more likely to be working with a more poorly sampled picture of reality, subject to measurement errors, biases and gaps in data. Can we still gain any benefit from attempting to create phase portraits?\n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":2367,"value":"{\n const maxPop = Math.max(...fish.map((f) => f[1]));\n const [w, h, hPad, vPad] = [600, 300, 30, 20];\n\n const xSc = d3\n .scaleTime()\n .domain([new Date(1972, 0, 1), new Date(2021, 11, 31)])\n .range([hPad, w]);\n const ySc = d3\n .scaleLinear()\n .domain([0, maxPop])\n .range([vPad, h - vPad]);\n\n const smoothed = regressionGenerator(fishSample).map(([_, smPop], i) => [\n fishSample[i][0],\n Math.max(0, smPop)\n ]);\n\n display(\n new Renderer(w, h)\n .setScales(xSc, ySc)\n .flipY()\n .textSize(8)\n .axisBottom(xSc, ySc)\n .axisLeft(xSc, ySc)\n .stroke(\"firebrick\")\n .strokeWidth(2)\n .line(smoothed)\n .strokeWidth(0.5)\n .stroke(\"steelblue\")\n .line(fishSample)\n .stroke()\n .textAlign(\"right\")\n .fill(\"steelblue\")\n .text(\"Sampled data\", new Date(2020, 0, 1), 1300)\n .fill(\"firebrick\")\n .text(\"Smoothed\", new Date(2020, 0, 1), 1210)\n .render()\n );\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2358,"value":"{\n const change = (values, i) => (values[i + 1] - values[i - 1]) / 2;\n\n const [w, h, border] = [500, 400, 20];\n\n const nCropped = 1;\n\n const smoothed =\n noise === 0\n ? fishSample\n : regressionGenerator(fishSample)\n .map(([_, smPop], i) => [fishSample[i][0], Math.max(0, smPop)])\n .slice(nCropped, fishSample.length - nCropped);\n\n const pop = smoothed.map((f) => f[1]);\n const maxPop = Math.max(...pop);\n\n const deltaPop = pop\n .map((p, i) => [change(pop, i), p])\n .slice(1, pop.length - 1);\n\n const xSc = d3\n .scaleLinear()\n .domain(d3.extent(deltaPop.map((dp) => dp[0])))\n .range([border, w - border]);\n\n const ySc = d3\n .scaleLinear()\n .domain([0, maxPop])\n .range([border, h - border]);\n\n display(\n new Renderer(w, h)\n .setScales(xSc, ySc)\n .textSize(8)\n .fill()\n .flipY()\n .axisBottom(xSc, ySc)\n .axisLeft(xSc, ySc, [0, maxPop], 0)\n .stroke(\"firebrick\")\n .bezierLine(deltaPop, 0.1)\n .render()\n );\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1904,"value":"const noise = view(\n Inputs.range([0, 1000], {\n label: \"noise\",\n step: 1,\n value: 70\n })\n);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2400,"value":"This smaller sample (120 rather than 600 data points) with some measurement errors (which can be controlled with the 'noise' slider) is still sufficient to show the asymmetries in phase space as well as the attenuation spiral towards an apparent attractor at around 750.\n\nIf we have more noisy data, the phase portrait soon becomes hard to read. To combat this we can smooth the data with some [local regression](https://en.wikipedia.org/wiki/Local_regression) (try moving the slider below). Oversmoothing while maintaining the cyclicity of the data, tends to remove the more subtle, but important properties of the shape and amplitude of the curves. The becomes one of selecting sufficient smoothing to see reflect structure but not so much we alter the signal in the noise.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":2002,"value":"const bandwidth = view(\n Inputs.range([0.01, 0.1], {\n label: \"Smoothing\",\n step: 0.001,\n value: 0.045\n })\n);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2426,"value":"Ultimately we want to understand what is controlling the dynamics of our fishing data. Those asymmetries give us some clue. Something that has the greatest influence on change when fish stocks are high and that has a greater influence over the rate of decline than rate of growth.\n\nIn this simple example, that cause is a population of [predators](https://en.wikipedia.org/wiki/Lotka–Volterra_equations) that lag behind the fish population. The fish population must be high for the predator population to be able to grow, and therefore be in a position to reduce the fish stocks rapidly. Likewise when the fish population is low, this results in a lagged drop in predator population allowing the fish to grow again. The spiralling of the phase trajectory suggests this might tend towards some dynamic equilibrium over time.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":8,"value":"---\n\n### Appendices","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":630,"value":"#### Acknowledgements\n\nThe pendulum example is adapted largely from Grant Sanderson's [tourist's guide to differential equations](https://youtu.be/p_di4Zn4wz4). The use of phase portraits in information visualization was inspired by Danny Dorling's book [Slowdown](https://www.dannydorling.org/books/SLOWDOWN/index.html). I created this page in response to Andy Kirk's [twitter comment](https://twitter.com/visualisingdata/status/1476518279187415045) on Colin Angus's [Covid phase portrait](https://twitter.com/VictimOfMaths/status/1476258498484871168).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":856,"value":"#### Constants","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":372,"value":"const π = Math.PI;\nconst halfPi = π / 2;\nconst 𝜏 = 2 * π;\nconst g = 9.81;","pinned":true,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1786,"value":"#### Helper functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":161,"value":"const radians = (deg) => (deg * Math.PI) / 180;\nconst degrees = (rad) => (180 * rad) / Math.PI;","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2188,"value":"function dotProd([[ax, ay], [bx, by], [cx, cy]]) {\n const [ux, uy] = [bx - ax, by - ay];\n const [vx, vy] = [cx - bx, cy - by];\n return ux * -vy + uy * vx;\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2429,"value":"function normalSamples(length, seed = 0.42) {\n const normal = d3.randomNormal.source(d3.randomLcg(seed))(0, 1);\n return d3.range(length).map(() => normal());\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":863,"value":"#### Pendulum modelling","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":342,"value":"function changeVec(theta, thetaDot) {\n const vel = thetaDot;\n const accel = -mu * thetaDot - (g / (0.02 * l)) * Math.sin(theta);\n return [vel, accel];\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1314,"value":"const pGrid = (() => {\n const gridRes = 0.4;\n const grid = [];\n let row = 0;\n\n for (let thetaDot = -5; thetaDot < 5 + gridRes; thetaDot += gridRes) {\n grid[row] = [];\n let col = 0;\n for (let theta = -𝜏 - π; theta <= 𝜏 + π; theta += gridRes) {\n grid[row][col++] = changeVec(theta, thetaDot);\n }\n row++;\n }\n return grid;\n})();","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1352,"value":"const pxScale = d3\n .scaleLinear()\n .domain([0, pGrid[0].length - 1])\n .range([-𝜏 - π, 𝜏 + π]);\n\nconst pyScale = d3\n .scaleLinear()\n .domain([0, pGrid.length - 1])\n .range([-5, 5]);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":477,"value":"function psTrajectory(pos0, vel0, t) {\n const accel = (theta, thetaDot) =>\n -mu * thetaDot - (g / (0.02 * l)) * Math.sin(theta);\n let theta = pos0;\n let thetaDot = vel0;\n const traj = [[pos0, vel0]];\n const deltaT = 0.1;\n for (let tInc = 0; tInc <= t; tInc += deltaT) {\n theta += thetaDot * deltaT;\n thetaDot += accel(theta, thetaDot) * deltaT;\n traj.push([theta, thetaDot]);\n }\n return traj;\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1799,"value":"#### Fish Modelling","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1989,"value":"import {regressionLoess} from \"npm:d3-regression@2\";","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2019,"value":"const regressionGenerator = regressionLoess()\n .x((d) => d[0])\n .y((d) => d[1])\n .bandwidth(bandwidth);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2430,"value":"const fishNoise = normalSamples(fish.length);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1773,"value":"const fish = (() => {\n // const [alpha, beta, delta, gamma] = [1, 0.25, 0.4, 0.05]; // Uncomment for fixed parameters\n const fishFn = (fish, pred) => selLV.alpha * fish - selLV.beta * fish * pred;\n const predFn = (fish, pred) => selLV.gamma * fish * pred - selLV.delta * pred;\n const [fishPop, predPop] = [10, 10];\n const [t, dt] = [50, 1 / 12];\n const pops = dampedRungeKutta4(fishPop, predPop, fishFn, predFn, dt, t / dt).map(\n ([fish, pred]) => Math.round(fish * 100)\n );\n\n const tParser = d3.utcParse(\"%Y-%m-%d\");\n return pops.map((pop, i) => [\n tParser(\n 1972 +\n Math.floor(i / 12) +\n \"-\" +\n `${(i % 12) + 1}`.padStart(2, \"0\") +\n \"-01\"\n ),\n Math.max(0, pop)\n ]);\n})();","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2353,"value":"const fishSample = fish\n .map(([date, pop], i) => [\n date,\n Math.max(0, Math.round(pop + noise * fishNoise[i]))\n ])\n .filter((_, i) => i % 5 === 0);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2329,"value":"_Lotka-Volterra parameters:_","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":2427,"value":"const selLV = view(\n Inputs.form({\n alpha: Inputs.range([0, 8], {\n label: \"alpha\",\n step: 0.01,\n value: 1.45\n }),\n beta: Inputs.range([0.01, 1], {\n label: \"beta\",\n step: 0.01,\n value: 0.25\n }),\n delta: Inputs.range([0.01, 3], {\n label: \"delta\",\n step: 0.01,\n value: 0.75\n }),\n gamma: Inputs.range([0, 2], {\n label: \"gamma\",\n step: 0.01,\n value: 0.15\n })\n })\n);","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1186,"value":"$39","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2288,"value":"const showControl = view(Inputs.toggle({label: \"Reveal\", value: false}));","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1747,"value":"#### ODE Solver","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1180,"value":"function dampedRungeKutta4(x, y, xFn, yFn, dt = 0.2, n = 500) {\n const rk4 = (x, y, dt, decl) => {\n const kx0 = xFn(x, y);\n const ky0 = yFn(x, y);\n const kx1 = xFn(x + 0.5 * dt * kx0, y + 0.5 * dt * ky0);\n const ky1 = yFn(x + 0.5 * dt * kx0, y + 0.5 * dt * ky0);\n const kx2 = xFn(x + 0.5 * dt * kx1, y + 0.5 * dt * ky1);\n const ky2 = yFn(x + 0.5 * dt * kx1, y + 0.5 * dt * ky1);\n const kx3 = xFn(x + dt * kx2, y + dt * ky2);\n const ky3 = yFn(x + dt * kx2, y + dt * ky2);\n return [\n x + (1 - decl * decl) * (dt / 6) * (kx0 + 2 * (kx1 + kx2) + kx3),\n y + (dt / 6) * (ky0 + 2 * (ky1 + ky2) + ky3)\n ];\n\n };\n\n const state = [[x, y]];\n let prev = state[0];\n\n for (let i = 1; i < n; i++) {\n const next = rk4(...prev, dt, i / n);\n state.push(next);\n prev = next;\n }\n return state;\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1325,"value":"#### Plotting functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":638,"value":"function pendulumPhaseSpace(w, h) {\n const xSc = d3\n .scaleLinear()\n .domain([-𝜏 - π, 𝜏 + π])\n .range([0, w]);\n const ySc = d3.scaleLinear().domain([-5, 5]).range([0, h]);\n const r = new Renderer(w, h)\n .setScales(xSc, ySc)\n .flipY()\n .line([\n [-𝜏 - π, 0],\n [𝜏 + π, 0]\n ])\n .line([\n [0, -5],\n [0, 4.6]\n ])\n .stroke()\n .textBaseline(\"top\")\n .textAlign(\"middle\")\n .text(\"velocity\", 0, 5)\n .textAlign(\"right\")\n\n .text(\"position\", radians(540), -0.08)\n .textAlign(\"middle\");\n\n return [-360, -180, -90, -45, 45, 90, 180, 360].reduce(\n (r, pos) => r.text(pos, radians(pos), -0.08),\n r\n );\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2123,"value":"function addVectorField(r, isColoured = false) {\n const vSc = 0.15;\n r.push().stroke(\"#aac\");\n const gridRes = 0.4;\n\n for (let thetaDot = -5; thetaDot < 5 + gridRes; thetaDot += gridRes) {\n for (let theta = -𝜏 - π; theta <= 𝜏 + π; theta += gridRes) {\n const [vel, accel] = changeVec(theta, thetaDot);\n const mag = Math.hypot(vel, accel);\n const [dx, dy] = [(vSc * vel) / mag, (vSc * accel) / mag];\n\n const traj = psTrajectory(theta, thetaDot, 0.3);\n if (isColoured) {\n r.stroke(dotProd(traj) < 0 ? \"steelblue\" : \"firebrick\");\n }\n r.line([\n [theta - dx, thetaDot - dy],\n [theta + dx, thetaDot + dy]\n ]).ellipse(theta + dx, thetaDot + dy, 0.01);\n }\n }\n return r.pop();\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1230,"value":"$3a","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":668,"value":"function addTrajectory(r, theta0, thetaDot0, tLen = 0.3, colour = null) {\n const traj = psTrajectory(radians(theta0), thetaDot0, tLen);\n return colour\n ? r.stroke(colour).line(traj)\n : r.stroke(dotProd(traj) < 0 ? \"#99c\" : \"#c99\").line(traj);\n}","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":14,"value":"import {Renderer} from \"observable:@jwolondon/renderer\";","pinned":false,"hidden":false,"mode":"ts","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W3b","stdlib":"2","titling":"auto"}} 33:"$32:body"