2c:null 2d:["$","title",null,{"children":"Estimate the Brachistochrone using Monte Carlo Simulation | Chuncheng | Observable"}] 32:T553,# Estimate the Brachistochrone using Monte Carlo Simulation The problem is finding a path that a ball falling with the shortest time costing. ## The ground truth Mathematically, the solution fits the parameter equation: ${tex.block`\begin{aligned} x &= R \cdot (\theta - sin(\theta)) \\ y &= R \cdot (1 - cos(\theta)) \end{aligned}`} where the ${tex`\theta \in [0, 2\pi)`} and ${tex`R \in \mathcal{R}^+`} ## Monte Carlo simulation There is a question, > Can we find the path without mathematical equations? So I make the x, y grids system, and fill it with gravity field. The graph is generated by the rules that the cells connect with their _8_ immediate neighbors. The weights of the edges are the time costing of passing it. As a result, the shortest path finding problem is converted into dijkstra path searching problem, which is fast and simple. ${tex.block`argmin_{p \in \mathcal{P}} T_{p}`} where ${tex`\mathcal{P}`} is the collection of the point inside the found path. **And I perform Monte Carlo simulation by randomizing the time costing between the grid cells.** It turns out that the simulation finds several near solutions. ## Conclusion It describes two things. The first is that the rect grids can not simulate the system of mathematical solution. The second is the Monte Carlo simulation can be used to find a comparable solution.33:T74c,canvas = { const context = DOM.context2d(width, height); const opacityPostfix = (k = 2) => ("0" + parseInt((255 / repeat) * k).toString(16)).slice(-2); var pathPoints; // Background { // context.fillStyle = "hsl(216deg 100% 13%)"; // context.fillRect(0, 0, width, height); context.canvas.style.background = "hsl(216deg 50% 13%)"; } // Draw points { function draw(pts, yOffset = 0, radius = 5, k = 2) { const maxV = d3.max(pts, (e) => e.v); pts.map((pts) => { const { x, y, v } = pts; (context.fillStyle = d3.color(colorMap(v / maxV)).hex() + opacityPostfix(k)), context.beginPath(), context.arc( constants.xScale(x), constants.yScale(y + yOffset), radius, 0, Math.PI * 2 ), context.fill(); }); } const radius = width / constants.numX / 3; if (showGridToggle) { draw(data.points, 0, radius / 2, 3); } pathList.map((path) => { pathPoints = path.map((str) => { return data.map[str]; }); draw(pathPoints, 0, radius); }); } // Draw ground truth { context.strokeStyle = "#ff0000" + "80"; context.lineWidth = 3; context.beginPath(); context.moveTo(constants.xScale(0), constants.yScale(0)); // context.lineTo(constants.xScale(0.2), constants.yScale(0.3)); const r = 1, n = 100; const thetaScale = d3 .scaleLinear() .domain([0, n * 2]) .range([0, Math.PI * 2]); var x, y, theta; for (let i = 0; i < n + 1; ++i) { theta = thetaScale(i); x = (r * (theta - Math.sin(theta))) / Math.PI; y = (r * (1 - Math.cos(theta))) / 2; context.lineTo(constants.xScale(x), constants.yScale(y)); } context.stroke(); } return context.canvas; }34:[] 30:{"id":"2b415b1a01bd9398","slug":"estimate-the-brachistochrone-using-monte-carlo-simulatio","title":"Estimate the Brachistochrone using Monte Carlo Simulation","thumbnail":"6e266720f14f7f026614638d4213ba00c326e9278674125a95eeca252df50866","access_level":"public","create_time":"$D2022-11-10T10:47:46.627Z","update_time":"$D2022-11-11T09:24:37.678Z","view_time":null,"first_public_time":"$D2022-11-10T13:41:25.335Z","fork":{"id":"7a0980291ec2d753","slug":"interactive-html-canvas","title":"Interactive HTML Canvas","owner":{"id":"f35c755083683fe5","login":"observablehq","name":"Observable"}},"stars":1,"starred":false,"owner":{"type":"team","id":"bcd65cef1cb4d2ed","login":"listenzcc","name":"Chuncheng","avatar_url":"https://avatars-next.observableusercontent.com/bb3c68e77a708a587ffc2d663bb83421f20d67991dbe6c1623c6e5d2fb73d15e","tier":"starter_2024"},"version":471,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"16328058352682a8","login":"listenzcc","name":"Chuncheng","avatar_url":"https://avatars-next.observableusercontent.com/bb3c68e77a708a587ffc2d663bb83421f20d67991dbe6c1623c6e5d2fb73d15e","tier":null}],"members":[{"type":"individual","id":"16328058352682a8","login":"listenzcc","name":"Chuncheng","avatar_url":"https://avatars-next.observableusercontent.com/bb3c68e77a708a587ffc2d663bb83421f20d67991dbe6c1623c6e5d2fb73d15e","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":3,"collection":{"id":"bb5b16d3c3bab332","slug":"monte-carlo-simulation","title":"Monte Carlo Simulation","owner":{"id":"bcd65cef1cb4d2ed","login":"listenzcc","name":"Chuncheng"}},"type":"notebook","body":{"title":"Estimate the Brachistochrone using Monte Carlo Simulation","cells":[{"id":0,"value":"$32","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":261,"value":"viewof repeat = Inputs.select([10, 15, 20, 50], { label: \"Repeat\", value: 10 })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":388,"value":"viewof paraRange = Inputs.range([1, 3], {\n label: \"Simulation range\",\n step: 0.1,\n value: 2\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":107,"value":"viewof numX = Inputs.range([50, 200], { label: \"numX\", value: 100, step: 10 })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":401,"value":"viewof showGridToggle = Inputs.toggle({ label: \"Grid toggle\", value: false })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":125,"value":"viewof palette = Inputs.select(palettes, {\n label: \"Color\",\n value: \"Viridis\"\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":292,"value":"{\n repeat;\n numX;\n paraRange;\n return html`\n

Simulation for ${repeat} times costs ${constants.simulationTime} ms

\n`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":4,"value":"$33","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":119,"value":"points","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":152,"value":"points = data.points","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":229,"value":"path = mkPath(mkLinks(data.points, data.map))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":244,"value":"pathList = {\n const t = new Date();\n const lst = [...new Array(repeat)];\n const pathLst = lst.map((e) => mkPath(mkLinks(data.points, data.map)));\n constants[\"simulationTime\"] = new Date() - t;\n return pathLst;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":226,"value":"function mkPath(links) {\n const { numX, numY } = constants;\n return dijkstra.find_path(links, \"0-0\", ij2str(numX, numY));\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":219,"value":"function mkLinks(points, map) {\n const links = {};\n var _pnt, s;\n points.map((pnt) => {\n const { i, j } = pnt;\n\n links[ij2str(i, j)] = {};\n\n for (let a = -1; a < 2; ++a) {\n for (let b = -1; b < 2; ++b) {\n if (a === 0 && b === 0) continue;\n\n _pnt = map[ij2str(i + a, j + b)];\n if (!_pnt) continue;\n\n s = a * b === 0 ? 1 : constants.sqrt2;\n\n links[ij2str(i, j)][ij2str(_pnt.i, _pnt.j)] =\n d3.randomUniform(1, paraRange)() * s * _pnt.dt;\n }\n }\n });\n return links;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":104,"value":"data = {\n const { numX, numY } = constants;\n\n const points = [];\n const map = {};\n\n var pnt, x, y, v, dt;\n\n // Build dataset with points and map\n for (let i = 0; i < numX + 1; ++i) {\n for (let j = 0; j < numY + 1; ++j) {\n x = i / numX;\n y = j / numY;\n v = computeV(y);\n dt = v === 0 ? 100 : 1 / v;\n\n pnt = {\n i,\n j,\n x,\n y,\n dt,\n v\n };\n\n map[ij2str(i, j)] = pnt;\n\n points.push(pnt);\n }\n }\n\n return { points, map };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":165,"value":"function str2ij(str) {\n const pair = str.split(\"-\").map((e) => parseInt(e));\n return { i: pair[0], j: pair[1] };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":162,"value":"function ij2str(i, j) {\n return i + '-' + j\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":95,"value":"function computeV(h) {\n const { mass, gravity } = constants;\n const v = Math.sqrt(2 * h * gravity);\n return v;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":99,"value":"constants = {\n return {\n sqrt2: Math.sqrt(2),\n mass: 1,\n gravity: 9.8,\n numX: numX,\n numY: parseInt((numX / width) * height),\n xScale: d3.scaleLinear().domain([-0.1, 1.1]).range([0, width]),\n yScale: d3.scaleLinear().domain([-0.1, 1.1]).range([0, height])\n };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":97,"value":"g = 9.8","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":6,"value":"height = (width * 9) / 16","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":123,"value":"colorMap = d3[`interpolate${palette}`]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":121,"value":"palettes = [\n \"RdBu\",\n \"PuOr\",\n \"Viridis\",\n \"Inferno\",\n \"Turbo\",\n \"Magma\",\n \"Plasma\",\n \"Cividis\",\n \"Warm\",\n \"Cool\",\n \"CubehelixDefault\",\n \"BuPu\",\n \"BuGn\",\n \"GnBu\",\n \"PuBu\",\n \"PuRd\",\n \"RdPu\",\n \"PuBuGn\",\n \"OrRd\",\n \"YlGnBu\",\n \"YlGn\",\n \"YlOrRd\",\n \"YlOrBr\"\n]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":142,"value":"dijkstra = require(\"https://bundle.run/dijkstrajs@1.0.2\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":91,"value":"d3 = require(\"d3\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W34","stdlib":"1","titling":"auto"}} 31:"$30:body"