30:T54e,# Phase Portraits 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...31:T572,{ 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 first = true; const r = new Renderer(w, h).textSize(14); while (first || !equals(thetaDoubleDot, 0) || !equals(thetaDot, 0)) { first = false; t += deltaT; thetaDoubleDot = accel(theta, thetaDot); thetaDot += thetaDoubleDot * deltaT; theta += thetaDot * deltaT; 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); yield Promises.delay(30, r.render()); } }32: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.33: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).34:T417,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). * 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).35:T43f,{ const change = (arr, i, w = 3) => { if (i <= 0 || i >= arr.length - 1) { return 0; } return (arr[i] - arr[i - 1] + (arr[i + 1] - arr[i])) / 2; }; const [w, h, border] = [500, 400, 20]; const nCropped = 1; const smoothed = noise === 0 ? fishSample : regressionGenerator(fishSample) .map(([_, smPop], i) => [fishSample[i][0], Math.max(0, smPop)]) .slice(nCropped, fishSample.length - nCropped); const pop = smoothed.map((f) => f[1]); const maxPop = Math.max(...pop); const deltaPop = pop .map((p, i) => [change(pop, i), p]) .slice(1, pop.length - 1); const xSc = d3 .scaleLinear() .domain(d3.extent(deltaPop.map((dp) => dp[0]))) .range([border, w - border]); const ySc = d3 .scaleLinear() .domain([0, maxPop]) .range([border, h - border]); return new Renderer(w, h) .setScales(xSc, ySc) .textSize(8) .fill() .flipY() .axisBottom(xSc, ySc) .axisLeft(xSc, ySc, [0, maxPop], 0) .stroke("firebrick") .bezierLine(deltaPop, 0.1) .render(); }36:T52a,{ // Lotka–Volterra preditor-prey equations: // const alpha = 1.1; // const beta = 0.4; // const delta = 0.4; // const gamma = 0.1; const preyFn = (prey, pred) => alpha * prey - beta * prey * pred; const predFn = (prey, pred) => gamma * prey * pred - 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 = rungeKutter4(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); } return r.render(); }37:T4a8,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 / 1) * 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)] ]); // r.textSize(6) // .fill("black") // .text(Math.round(reflect), pxScale(col), pyScale(row)); } } return r.pop(); }38:[] 2e:{"title":"Phase Portraits","cells":[{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"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"},{"id":17,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":297,"value":"viewof l = Inputs.range([75, 275], {\n label: \"Pendulum length\",\n step: 5,\n value: 200\n})\n\n// viewof l = slider({\n// title: \"Pendulum length\",\n// min: 75,\n// max: 275,\n// step: 5,\n// value: 200,\n// submit: true\n// })","pinned":false,"hidden":false,"mode":"ojs"},{"id":362,"value":"viewof mu = Inputs.range([0, 2], {\n label: \"Air resistance\",\n step: 0.02,\n value: 0.1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":876,"value":"viewof restart = Inputs.button(\"Restart\", {})","pinned":false,"hidden":false,"mode":"ojs"},{"id":74,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":693,"value":"{\n const r = pendulumPhaseSpace(900, 320);\n addVectorField(r);\n return r.render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"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"},{"id":696,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"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 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"},{"id":2270,"value":"Inputs.bind(\n Inputs.range([75, 275], { label: \"Pendulum length\", step: 5 }),\n viewof l\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":553,"value":"viewof thetaDot0 = Inputs.range([-4, 4], {\n label: \"Initial velocity\",\n step: 0.01,\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":530,"value":"viewof theta0 = Inputs.range([-180, 180], {\n label: \"Initial position\",\n step: 1,\n value: 60\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":985,"value":"Inputs.bind(\n Inputs.range([0, 2], { label: \"Air resistance\", step: 0.02 }),\n viewof mu\n)","pinned":false,"hidden":false,"mode":"ojs"},{"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 return r.render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"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"},{"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"},{"id":726,"value":"{\n const rand = (min = 0, max = 1) => Math.random() * (max - min) + min;\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 return r.render();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2280,"value":"viewof showPendulum = Inputs.toggle({ label: \"Show pendulum\", value: false })","pinned":false,"hidden":false,"mode":"ojs"},{"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-Kutter methods](https://en.wikipedia.org/wiki/Runge–Kutta_methods)).","pinned":false,"hidden":false,"mode":"md"},{"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"},{"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"},{"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\n const ySc = d3\n .scaleLinear()\n .domain([0, maxPop])\n .range([vPad, h - vPad]);\n\n const pop = fish.map((f) => f[1]);\n\n return 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}","pinned":false,"hidden":false,"mode":"ojs"},{"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"},{"id":1866,"value":"{\n const change = (arr, i, w = 3) => {\n if (i <= 0 || i >= arr.length - 1) {\n return 0;\n }\n return (arr[i] - arr[i - 1] + (arr[i + 1] - arr[i])) / 2;\n };\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 return 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}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2320,"value":"$34","pinned":false,"hidden":false,"mode":"md"},{"id":2343,"value":"### NeedABiggerBoat Fisheries\n\nUnfortunately, we tend not to have the resources of BigBuckFishingCo 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"},{"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\n const ySc = d3\n .scaleLinear()\n .domain([0, maxPop])\n .range([vPad, h - vPad]);\n\n const pop = fishSample.map((f) => f[1]);\n\n const smoothed = regressionGenerator(fishSample).map(([_, smPop], i) => [\n fishSample[i][0],\n Math.max(0, smPop)\n ]);\n\n return 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}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2358,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":1904,"value":"viewof noise = Inputs.range([0, 1000], {\n label: \"noise\",\n step: 1,\n value: 70\n})","pinned":false,"hidden":false,"mode":"ojs"},{"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"},{"id":2002,"value":"viewof bandwidth = Inputs.range([0, 0.1], {\n label: \"smoothing\",\n step: 0.001,\n value: 0.045\n})","pinned":false,"hidden":false,"mode":"ojs"},{"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"},{"id":8,"value":"---\n\n### Appendices","pinned":false,"hidden":false,"mode":"md"},{"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"},{"id":856,"value":"#### Constants","pinned":false,"hidden":false,"mode":"md"},{"id":372,"value":"pi = Math.PI","pinned":false,"hidden":false,"mode":"ojs"},{"id":1291,"value":"halfPi = Math.PI / 2","pinned":false,"hidden":false,"mode":"ojs"},{"id":446,"value":"tau = 2 * Math.PI","pinned":false,"hidden":false,"mode":"ojs"},{"id":356,"value":"g = 9.81","pinned":false,"hidden":false,"mode":"ojs"},{"id":1786,"value":"#### Helper functions","pinned":false,"hidden":false,"mode":"md"},{"id":161,"value":"function radians(deg) {\n return (deg * Math.PI) / 180;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1404,"value":"function degrees(rad) {\n return (180 * rad) / Math.PI;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1781,"value":"function rand(mn = 0, mx = 1) {\n return Math.random() * (mx - mn) + mn;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1790,"value":"function gauss() {\n let [u, v] = [0, 0];\n while (u === 0) {\n u = Math.random();\n }\n while (v === 0) {\n v = Math.random();\n }\n return Math.sqrt(-2.0 * Math.log(u)) * Math.cos(2.0 * Math.PI * v);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"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":"ojs"},{"id":863,"value":"#### Pendulum modelling","pinned":false,"hidden":false,"mode":"md"},{"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":"ojs"},{"id":1314,"value":"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 = -tau - pi; theta <= tau + pi; theta += gridRes) {\n grid[row][col++] = changeVec(theta, thetaDot);\n }\n row++;\n }\n return grid;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1352,"value":"pxScale = d3\n .scaleLinear()\n .domain([0, pGrid[0].length - 1])\n .range([-tau - pi, tau + pi])","pinned":false,"hidden":false,"mode":"ojs"},{"id":1356,"value":"pyScale = d3\n .scaleLinear()\n .domain([0, pGrid.length - 1])\n .range([-5, 5])","pinned":false,"hidden":false,"mode":"ojs"},{"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":"ojs"},{"id":1799,"value":"#### Fish Modelling","pinned":false,"hidden":false,"mode":"md"},{"id":1989,"value":"regress = require(\"d3-regression@1\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":2019,"value":"regressionGenerator = regress\n .regressionLoess()\n .x((d) => d[0])\n .y((d) => d[1])\n .bandwidth(bandwidth)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1773,"value":"fish = {\n // const [alpha, beta, delta, gamma] = [1, 0.25, 0.4, 0.05]; // Uncomment for fixed parameters\n const fishFn = (fish, pred) => alpha * fish - beta * fish * pred;\n const predFn = (fish, pred) => gamma * fish * pred - delta * pred;\n const [fishPop, predPop] = [10, 10];\n const [t, dt] = [50, 1 / 12];\n const pops = rungeKutter4(\n fishPop,\n predPop,\n fishFn,\n predFn,\n dt,\n t / dt\n ).map(([fish, pred]) => Math.round(fish * 100));\n\n const times = [];\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":"ojs"},{"id":2353,"value":"fishSample = fish\n .map(([date, pop]) => [date, Math.max(0, Math.round(pop + noise * gauss()))])\n .filter((datum, i) => i % 5 === 0)","pinned":false,"hidden":false,"mode":"ojs"},{"id":2329,"value":"_Lotka-Volterra parameters:_","pinned":false,"hidden":false,"mode":"md"},{"id":1669,"value":"viewof alpha = Inputs.range([0, 8], {\n label: \"alpha\",\n step: 0.01,\n value: 1.45\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1675,"value":"viewof beta = Inputs.range([0.01, 1], {\n label: \"beta\",\n step: 0.01,\n value: 0.25\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1681,"value":"viewof delta = Inputs.range([0.01, 3], {\n label: \"delta\",\n step: 0.01,\n value: 0.75\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1686,"value":"viewof gamma = Inputs.range([0, 2], {\n label: \"gamma\",\n step: 0.01,\n value: 0.15\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1186,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":2288,"value":"viewof showControl = Inputs.toggle({ label: \"Reveal\", value: false })","pinned":false,"hidden":false,"mode":"ojs"},{"id":1747,"value":"#### ODE Solver","pinned":false,"hidden":false,"mode":"md"},{"id":1180,"value":"function rungeKutter4(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 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":"ojs"},{"id":1325,"value":"#### Plotting functions","pinned":false,"hidden":false,"mode":"md"},{"id":638,"value":"function pendulumPhaseSpace(w, h) {\n const xSc = d3\n .scaleLinear()\n .domain([-tau - pi, tau + pi])\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 [-tau - pi, 0],\n [tau + pi, 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":"ojs"},{"id":2123,"value":"function addVectorField(r, isColoured = false) {\n const vSc = 0.15;\n\n r.push().stroke(\"#aac\");\n\n const gridRes = 0.4;\n\n for (let thetaDot = -5; thetaDot < 5 + gridRes; thetaDot += gridRes) {\n for (let theta = -tau - pi; theta <= tau + pi; 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":"ojs"},{"id":1230,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":668,"value":"function addTrajectory(r, theta0, thetaDot0, tLen = 0.3, colour = null) {\n const traj = psTrajectory(radians(theta0), thetaDot0, tLen);\n if (colour) {\n return r.stroke(colour).line(traj);\n }\n return r.stroke(dotProd(traj) < 0 ? \"#99c\" : \"#c99\").line(traj);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":14,"value":"import { Renderer } from \"@jwolondon/renderer\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W38","stdlib":"1","titling":"auto"}