2c:null 2d:["$","title",null,{"children":"A Prime on Web-Based Simulations | Icaro Fonseca | Observable"}] 32:T40a,md`# A Prime on Web-Based Simulations Support material for the paper: http://www.scs-europe.net/dlib/2019/2019-0023.htm, coauthored with [@hmgaspar](https://observablehq.com/@hmgaspar). We formulate a simple problem of a block vibrating in one direction while attached to a spring and damper, i.e., a damped harmonic oscillator. The block is subject to an external force, as shown in figure below. ![](https://shiplab.github.io/vesseljs/images/mass-spring-damper.jpg) The motion of the block over time can be expressed with the equation: F = mẍ + cẋ + kx This equation can be solved in JavaScript using the [Runge-Kutta Dormand-Prince solver (RKDP)](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method), implemented in the [Numeric.js library](https://github.com/sloisel/numeric). To do so, we need to provide the solver the initial time, final time, initial conditions (position and velocity) and a function containing the equation to be solved. To model the oscillator, we start defining the constant parameters.`33:T535,{ function forcedMotion(t, x) { var F = 10 * Math.sin(t / 5); var xdotdot = (F - c * x[1] - k * x[0]) / m; return [x[1], xdotdot]; }; var initCond = [0, 0]; var sol_forced = numeric.dopri(0, tmax, initCond, forcedMotion); var f_motion_force = numeric.rep([sol_forced.x.length, 2], 0); for (var i = 0; i < sol_forced.x.length; i++) { f_motion_force[i][0] = sol_forced.x[i]; f_motion_force[i][1] = 10 * Math.sin(sol_forced.x[i] / 5); } $.plot("#f_motion_force", [f_motion_force], { xaxis: { tickFormatter: function(val, axis) { return val < axis.max ? val.toFixed(2) : "time (s)"; } }, yaxis: { tickFormatter: function(val, axis) { return val < axis.max ? val.toFixed(2) : "force (N)"; }, } }); var x_motion_force = numeric.rep([sol_forced.x.length, 2], 0); for (var i = 0; i < sol_forced.x.length; i++) { x_motion_force[i][0] = sol_forced.x[i]; x_motion_force[i][1] = sol_forced.y[i][0]; } $.plot("#x_motion_force", [x_motion_force], { xaxis: { tickFormatter: function(val, axis) { return val < axis.max ? val.toFixed(2) : "time (s)"; } }, yaxis: { tickFormatter: function(val, axis) { return val < axis.max ? val.toFixed(2) : "position (m)"; }, } }); }34:[] 30:{"id":"3cb45332585d9e92","slug":"a-prime-on-web-based-simulations","title":"A Prime on Web-Based Simulations","thumbnail":null,"access_level":"public","create_time":"$D2020-01-15T17:52:36.308Z","update_time":"$D2021-08-08T18:25:34.140Z","view_time":null,"first_public_time":"$D2020-01-15T18:37:27.185Z","has_unpublished_edits":null,"fork":false,"stars":1,"starred":false,"owner":{"type":"team","id":"45fe56d5093061eb","login":"icarofonseca","name":"Icaro Fonseca","avatar_url":"https://avatars-next.observableusercontent.com/74086e355d9bde53fd85227b6dd08934b5747470968c89184341654aa47c5050","tier":"starter_2024","bio":""},"version":550,"publish_version":550,"role":null,"authored":false,"authors":[{"type":"individual","id":"a91dfff54ae76357","login":"hmgaspar","name":"Henrique Gaspar","avatar_url":"https://avatars-next.observableusercontent.com/d261fd8e8091b2374febc872ea401770bf95cd78a7f0f0dbf1aaf85e023f116e","tier":null,"bio":""},{"type":"individual","id":"e2c2510ce71e416b","login":"icarofonseca","name":"Icaro Fonseca","avatar_url":"https://avatars-next.observableusercontent.com/74086e355d9bde53fd85227b6dd08934b5747470968c89184341654aa47c5050","tier":null,"bio":""}],"members":[{"type":"individual","id":"e2c2510ce71e416b","login":"icarofonseca","name":"Icaro Fonseca","avatar_url":"https://avatars-next.observableusercontent.com/74086e355d9bde53fd85227b6dd08934b5747470968c89184341654aa47c5050","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":1,"collection":null,"type":"notebook","body":{"titling":"auto","title":"A Prime on Web-Based Simulations","cells":[{"id":0,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":8,"value":"m = 30 // kg","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":15,"value":"k = 20 // N/m","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":18,"value":"c = 2.9 // (N*s)/m","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":205,"value":"md`Now the function representing the differential equation needs to be defined. To do that, the previous equation is rearranged in a system of two equations.\n\nd(x)/dt = ẋ\n\nd(ẋ)/dt = ẍ = (F - cẋ - kx)/m\n\nA function is defined to receive as input an array with the position and velocity, respectively x and ẋ (the arguments of the derivatives on the left side of the equations) and output the evaluated derivatives (as on the right side).`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":24,"value":"function freeMotion(t, x) {\n var F = 0; // N\n var xdotdot = (F - c * x[1] - k * x[0]) / m;\n return [x[1], xdotdot];\n};","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":261,"value":"md`The initial conditions are set to x = 1 m and ẋ = 0 m/s and the maximum simulation time to 100 seconds.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":158,"value":"initX = [1, 0] // x = [initial position in m, initial velocity in m/s]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":171,"value":"tmax = 100 // s","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":273,"value":"md`Finally, the following lines solve the equation. The RKDP method uses adaptive time-step, so no standard time-step size is specified. If necessary, the default tolerance can be tuned.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":36,"value":"sol_free = numeric.dopri(0, tmax, initX, freeMotion);","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":262,"value":"md`This yields an object \\`\\`sol_free\\`\\` with all the positions, velocities and accelerations calculated by the solver.\n\n### Graphical Visualization\n\nBesides inspecting the simulation results in the cell above, it is also possible to create visualizations of relevant simulation variables. The plots presented in this section use [Flot](https://github.com/flot/flot) due to its good compatibility with Numeric.js. Other useful visualization libraries are [D3.js](https://d3js.org) and [plotly.js](https://plot.ly/). The cell below shows a plot of the block position x in relation to the neutral position over time as calculated for the free oscillator.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":297,"value":"html`

x position (free oscillation)

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":305,"value":"{\n const x_motion = numeric.rep([sol_free.x.length, 2], 0);\n \n for (var i = 0; i < sol_free.x.length; i++) {\n x_motion[i][0] = sol_free.x[i];\n x_motion[i][1] = sol_free.y[i][0];\n }\n \n $.plot(\"#x_motion\", [x_motion], {\n xaxis: {\n tickFormatter: function(val, axis) {\n return val < axis.max ? val.toFixed(2) : \"time (s)\";\n }\n },\n yaxis: {\n tickFormatter: function(val, axis) {\n return val < axis.max ? val.toFixed(2) : \"position (m)\";\n },\n }\n });\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":316,"value":"md`It is also possible to simulate and visualize the motion of a forced oscillator. The force is defined as a sinusoidal with the equation:\n\nF = 10 sin(t/5)\n\nBy making both the initial position and velocity equal to zero and solving the problem again with the new force, we obtain the plots for the force and position.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":149,"value":"html`

x position (forced oscillation)

force (forced oscillation)

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":245,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":287,"value":"md`### 3D Visualization\n\nThe WebGL library allows rendering of 2D and 3D graphics and textures in any compatible browser. The [Three.js](https://threejs.org/) library is used to aid the creation and manipulation of 3D objects with WebGL, reducing the effort necessary to build a meaningful visualization.\n\nYou can **restart the simulation** by triggering the clock cell below.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":498,"value":"clock = new THREE.Clock()","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":501,"value":"mutable time = 0","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":46,"value":"{\n const renderer = new THREE.WebGLRenderer({antialias: true});\n invalidation.then(() => renderer.dispose());\n renderer.setSize(width, height);\n renderer.setPixelRatio(devicePixelRatio);\n const controls = new THREE.OrbitControls(camera, renderer.domElement);\n\n var tprev = 0;\n var posX = initX;\n\n while (true) {\n mutable time = clock.getElapsedTime();\n if (time !== tprev) {\n // solve ODE for the time step\n posX = numeric.dopri(tprev, time, posX, freeMotion).at(time);\n cube.position.x = posX[0];\n tprev = time;\n }\n\n spring.position.x = (cube.position.x + base.position.x) / 2;\n spring.scale.y = (cube.position.x - base.position.x) / 2;\n\n damper.position.x = (cube.position.x + base.position.x) / 2;\n damper.scale.y = (cube.position.x - base.position.x) / 2;\n\n renderer.render(scene, camera);\n yield renderer.domElement;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":410,"value":"md`The scene contains all the pieces of the oscillator: a cube attached to a base with a spring and a damper represented with a cylindrical geometry.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":336,"value":"scene = {\n const scene = new THREE.Scene();\n scene.background = new THREE.Color(0xA9CCE3);\n scene.add(cube);\n scene.add(base);\n scene.add(spring);\n scene.add(damper);\n return scene;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":434,"value":"md`The individual pieces are created with geometries and positioned on the scene.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":367,"value":"cube = {\n //Add cube\n const geometryCube = new THREE.BoxGeometry(1, 1, 1);\n const materialCube = new THREE.MeshBasicMaterial({\n color: \"blue\"\n });\n return new THREE.Mesh(geometryCube, materialCube);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":369,"value":"base = {\n //Add base\n const geometryBase = new THREE.BoxGeometry(0.05, 2, 2);\n const materialBase = new THREE.MeshBasicMaterial({\n color: \"green\"\n });\n const base = new THREE.Mesh(geometryBase, materialBase);\n base.position.x = -2;\n return base;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":381,"value":"geometryCyl = new THREE.CylinderGeometry(0.05, 0.05, 2, 32);","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":371,"value":"spring = {\n const materialSpring = new THREE.MeshBasicMaterial({\n color: \"yellow\"\n });\n const spring = new THREE.Mesh(geometryCyl, materialSpring);\n spring.position.x = -1;\n spring.position.z = 0.25;\n spring.rotation.z = Math.PI / 2;\n return spring;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":373,"value":"damper = {\n const materialDamper = new THREE.MeshBasicMaterial({\n color: \"red\"\n });\n const damper = new THREE.Mesh(geometryCyl, materialDamper);\n damper.position.x = -1;\n damper.position.z = -0.25;\n damper.rotation.z = Math.PI / 2;\n return damper;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":431,"value":"md`The camera is adjusted to show the oscillator on the horizontal position.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":339,"value":"camera = {\n const fov = 50;\n const aspect = width / height;\n const near = 1;\n const far = 1000;\n const camera = new THREE.PerspectiveCamera(fov, aspect, near, far);\n camera.up.set(0, 0, 1);\n scene.add(camera);\n camera.position.set(0.5, -5, 0);\n return camera;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":423,"value":"md`The height of the rendered scene is fixed and the width is adjusted automatically.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":341,"value":"height = 600","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":61,"value":"md`### Libraries`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":50,"value":"THREE = {\n const THREE = window.THREE = await require(\"three@0.99.0/build/three.min.js\");\n await require(\"three@0.99.0/examples/js/controls/OrbitControls.js\").catch(() => {});\n await require(\"three@0.99.0/examples/js/loaders/STLLoader.js\").catch(() => {});\n return THREE;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":68,"value":"numeric = require(\"numeric@1.2.6/numeric-1.2.6.js\").catch(() => window[\"numeric\"]);","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":131,"value":"$$ = {\n const $ = window.$ = await require('jquery@1.11.1/dist/jquery.js').catch(() => window[\"$\"]);\n await require(\"flot@4.1.1/dist/es5/jquery.flot.js\").catch(() => window[\"_typeof\"]);\n return $;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W34","stdlib":"1"}} 31:"$30:body"