2b:null 2c:["$","title",null,{"children":"Comparing Numerical Methods, part 1. | Ron Smith | Observable"}] 31:T917,function compareTheMethods(fun,params,analyticSoln) { const xMin = 0, xMax = tMax, yMin = -10, yMax = 10; const margin = {top: 20, bottom: 20, right: 20, left: 20}; const height = 400; const width = 400; const w = width - margin.right - margin.left; const h = height - margin.top - margin.bottom; const x = d3.scaleLinear().domain([xMin,xMax]).range([0,w]); const y = d3.scaleLinear().domain([yMin,yMax]).range([h,0]); let svg = d3.select(DOM.svg(width, height)) .append("g").attr("transform", "translate(" + margin.left + "," + margin.top + ")"); svg.append("g") .attr("transform", "translate(0," + h + ")") .call(d3.axisBottom(x)); svg.append("g") .call(d3.axisLeft(y)); let path_d, t, temp; // If no analytic solution was provided, skip this block. if (analyticSoln) { // Plot actual solution let dt_an = 0.01; // Overrides the time step for the numerical solution. Aliasing is a neat topic and a sort of approximation failure in itself, but not for this example. let pathA = svg.append("path") .attr("fill", "none") .attr("stroke", "blue") .attr("stroke-width", 3.7); let TA=[],XA=[]; t=0;temp=[0,x0]; while (t<=(tMax+dt_an)) { temp = analyticSoln(temp,params); TA.push(temp[0]); XA.push(temp[1]); t+=dt_an; temp[0]=t; } path_d = make_svg_path_from_data(TA,XA,x,y); pathA.attr("d", path_d); } // Approximate with Euler's method let pathE = svg.append("path") .attr("fill", "none") .attr("stroke", "red") .attr("stroke-width", 1.3); let TE=[],XE=[]; t=0; temp=[0,x0]; TE.push(0); XE.push(x0); while (t<=tMax) { temp = Euler(fun, temp, params, dt); t=temp[0]; TE.push(t); XE.push(temp[1]); } path_d = make_svg_path_from_data(TE,XE,x,y); pathE.attr("d", path_d); // Approximate with Heun's method let pathH = svg.append("path") .attr("fill", "none") .attr("stroke", "orange") .attr("stroke-width", 1.7); let TH=[],XH=[]; t=0; temp=[0,x0]; TH.push(0); XH.push(x0); while (t<=tMax) { temp = Heun(fun, temp, params, dt); t=temp[0]; TH.push(t); XH.push(temp[1]); } path_d = make_svg_path_from_data(TH,XH,x,y); pathH.attr("d", path_d); //return [TA,XA]; return svg.node().parentNode; }32:T91d,function compareTheMethods_ZoomZoom(fun,params,analyticSoln) { const xMin = 4, xMax = 8, yMin = -4, yMax = 4; const margin = {top: 20, bottom: 20, right: 20, left: 20}; const height = 400; const width = 400; const w = width - margin.right - margin.left; const h = height - margin.top - margin.bottom; const x = d3.scaleLinear().domain([xMin,xMax]).range([0,w]); const y = d3.scaleLinear().domain([yMin,yMax]).range([h,0]); let svg = d3.select(DOM.svg(width, height)) .append("g").attr("transform", "translate(" + margin.left + "," + margin.top + ")"); svg.append("g") .attr("transform", "translate(0," + h + ")") .call(d3.axisBottom(x)); svg.append("g") .call(d3.axisLeft(y)); let path_d, t, temp; // If no analytic solution was provided, skip this block. if (analyticSoln) { // Plot actual solution let dt_an = 0.01; // Overrides the time step for the numerical solution. Aliasing is a neat topic and a sort of approximation failure in itself, but not for this example. let pathA = svg.append("path") .attr("fill", "none") .attr("stroke", "blue") .attr("stroke-width", 3.7); let TA=[],XA=[]; t=0;temp=[0,x0]; while (t<=(tMax+dt_an)) { temp = analyticSoln(temp,params); TA.push(temp[0]); XA.push(temp[1]); t+=dt_an; temp[0]=t; } path_d = make_svg_path_from_data(TA,XA,x,y); pathA.attr("d", path_d); } // Approximate with Euler's method let pathE = svg.append("path") .attr("fill", "none") .attr("stroke", "red") .attr("stroke-width", 1.3); let TE=[],XE=[]; t=0; temp=[0,0.5]; TE.push(0); XE.push(x0); while (t<=tMax) { temp = Euler(fun, temp, params, ds); t=temp[0]; TE.push(t); XE.push(temp[1]); } path_d = make_svg_path_from_data(TE,XE,x,y); pathE.attr("d", path_d); // Approximate with Heun's method let pathH = svg.append("path") .attr("fill", "none") .attr("stroke", "orange") .attr("stroke-width", 1.7); let TH=[],XH=[]; t=0; temp=[0,0.5]; TH.push(0); XH.push(x0); while (t<=tMax) { temp = Heun(fun, temp, params, ds); t=temp[0]; TH.push(t); XH.push(temp[1]); } path_d = make_svg_path_from_data(TH,XH,x,y); pathH.attr("d", path_d); //return [TA,XA]; return svg.node().parentNode; }33:[] 2f:{"id":"0e84dcb35e29e8fe","slug":"comparing-numerical-methods-part-1","title":"Comparing Numerical Methods, part 1.","thumbnail":"439bc7961a2e438612aa77426127b916c1c2959e45c97cef1814698877d7b1f2","access_level":"public","create_time":"$D2019-01-27T15:49:02.241Z","update_time":"$D2019-03-08T15:31:24.212Z","view_time":null,"first_public_time":"$D2019-01-27T22:00:12.968Z","fork":false,"stars":1,"starred":false,"owner":{"type":"team","id":"8dab63897b6adf13","login":"ronzor","name":"Ron Smith","avatar_url":"https://avatars-next.observableusercontent.com/72c20176307cf53407d67a22afcc4d7d49bbdd419afc85de1cf01fe8aa2e320f","tier":"starter_2024"},"version":912,"publish_version":912,"role":null,"authored":false,"authors":[{"type":"individual","id":"a5eeba816274b983","login":"ronzor","name":"Ron Smith","avatar_url":"https://avatars-next.observableusercontent.com/72c20176307cf53407d67a22afcc4d7d49bbdd419afc85de1cf01fe8aa2e320f","tier":null}],"members":[{"type":"individual","id":"a5eeba816274b983","login":"ronzor","name":"Ron Smith","avatar_url":"https://avatars-next.observableusercontent.com/72c20176307cf53407d67a22afcc4d7d49bbdd419afc85de1cf01fe8aa2e320f","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":2,"collection":null,"type":"notebook","body":{"title":"Comparing Numerical Methods, part 1.","cells":[{"id":0,"value":"md`# Comparing Numerical Methods, part 1.\n\nThis is part of an ongoing exercise where I practice JavaScript and try to leave something useful behind in the process. The goal is: You can completely ignore the code and still possibly learn something, or you can read the code and learn that part too (albeit from a newbie).\n\nWe'll compare ${html`Euler's method`} and ${html`Heun's method`}. The goal is to develop a visual intuition for how they might deviate from a true solution.\n\n\n\nWe'll use the differential equation\n${tex.block `\\frac{dx}{dt} = a_1\\cos(w_1t),`}\nthe analytic solution of which is\n${tex.block `x(t) = x_0 + \\frac{a_1}{w_1}\\sin(w_1t).`}\n\nAs a reminder, you can click on the vertical dots to the left of each cell to see the code for that cell. You can also experiment and edit the code.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":295,"value":"md`--- ---\nFirst, let's define the differential equation,`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":304,"value":"function dxdt(x,params) {\n // Pass all inputs and outputs as arrays, even if just one entry.\n // x = [t,x1,x2,...,xN];\n // t = x[0] is the zeroth state, and dt/dt=1,\n // so we're returning [dt/dt=1, dx/dt=a*cos(w*t)],\n return [1, params[0]*Math.cos(params[1]*x[0])];\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":310,"value":"md`and, in this case, our analytic solution,`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":124,"value":"function xActual(x, params) {\n // All inputs and outputs as arrays\n // Even though this example is a function of one variable, we're still returning an array where [x0=t,x1=x]. \n // It'll help keep the code looking more consistent later and in future examples with systems of equations.\n return [x[0], params[2] + (params[0]*Math.sin(params[1]*x[0]))/params[1]];\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":312,"value":"md`--- ---\nNext, define Euler's method. The function will evaluate one step; to get the entire solution we'll use a loop later.\n\nEuler's method and, more specifically, what the next function should be doing, is computing\n${tex.block`\\mathbf{x}_{t+\\Delta t} = \\mathbf{x}_t + f(\\mathbf{x}_t)\\Delta t`}\nwhere ${tex`\\mathbf{x}_t = \\begin{bmatrix} t \\\\ x \\end{bmatrix}`} is a vector representing the state of the system at time ${tex`t`}, and ${tex.block`f(\\mathbf{x}_t) = \\dot{\\mathbf{x}}_t = \\begin{bmatrix} 1 \\\\ \\dot{x}(t)\\end{bmatrix}`}\n`\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":25,"value":"function Euler(fun,state,params,dt) {\n // All inputs and output as ararys\n // state = [t,x1,x2,...,xN];\n // Evaluate Euler's method at the current state, with step size dt, and return the next state.\n let deriv = fun(state,params);\n let y = [];\n for (let i=0; iblue`}, Euler's method is ${html`red`}, and Heun's method is ${html`orange`}. The analytic solution is ${html` ~~evaluated at the same time steps as the numerical approximations and each are connected by straight lines~~`} evaluated with constant step size ${tex`dt=0.01`} to avoid ${html`aliasing`}.\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":111,"value":"viewof x0 = slider({\n min: -2, \n max: 2, \n step: 0.1, \n value: 0, \n title: tex `x_0`\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":85,"value":"viewof a1 = slider({\n min: -10, \n max: 10, \n step: 0.1, \n value: 1, \n title: tex `a_1`\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":101,"value":"viewof w1 = slider({\n min: -10, \n max: 10, \n step: 0.1, \n value: 1, \n title: tex `w_1`\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":108,"value":"viewof dt = slider({\n min: 0.01, \n max: 1, \n step: 0.01, \n value: 0.05, \n title: tex `dt`\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":118,"value":"viewof tMax = slider({\n min: 1, \n max: 50, \n step: 1, \n value: 25, \n title: tex `t_{max}`\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":330,"value":"compareTheMethods(dxdt,[a1,w1,x0],xActual)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":388,"value":"md`Have some fun with that one. You might already know that some numerical methods are better than others, and that Euler's is one of the most basic and least accurate... The point is to visually explore the interesting ways in which these methods deviate from the true answer. We often use numerical systems to study dynamical systems, but a numerical method is a dynamical system in and of itself! \n\n--- ---\n# Exercise\nTry the following parameters: ${tex`(x_0,a_1,w_1,dt,t_{max}) = (0,10,10,0.33,50)`}. \n\nCheck what happens when you \"wiggle\" them (especially ${tex`dt`}) a little from there. Although ${tex`dt=0.33`} is an impractically large time step in most applications, and words like unstable might come to mind... did you expect to see things like that? \n\nAside from knowing that our time step is too large, would you be able to diagnose this problem if you didn't have the analytic solution as a reference?\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":582,"value":"md` Enjoy,\n\n_ronzor_ ~ 2019 Jan 27`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":323,"value":"md`--- ---\nBelow are some other stuffs used to make things like plots and slider bars happen:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":39,"value":"d3 = require(\"d3@5\");","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":42,"value":"function make_svg_path_from_data(xt, yt, X, Y) {\n //Inputs:\n // xt,yt are arrays of the same size containing data for the horizontal and vertical axes respectively.\n // X,Y are d3.linearScale objects that map our coordinates to the target SVG's coordinates.\n \n // Function returns a string containing an SVG path.\n var path_d = 'M' + X(xt[0]) + ',' + Y(yt[0]);\n for (var i = 1; i < xt.length; i++) {\n path_d = path_d + ' L' + X(xt[i]) + ',' + Y(yt[i]);\n }\n return (path_d)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":87,"value":"import {slider} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":826,"value":"md`--- ---\nundocumented bonus materials \n---\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":842,"value":"ds = {\n let t = 0;\n let dt = 0.0005;\n while (true) {\n t = 0.6 + (t + dt)%0.6;\n yield t.toFixed(4);\n }\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":836,"value":"compareTheMethods_ZoomZoom(dxdt,[-4,2,0.5],xActual)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":831,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W33","stdlib":"1","titling":"auto"}} 30:"$2f:body"