31:Tf60,slope_field = { let f_parsed, f; d3.select(parse_warning).html(''); try { try { f_parsed = math.compile(f1); f = (y, t) => f_parsed.evaluate({ y: y, t: t }); } catch (e1) { d3.select(parse_warning).html( 'Parse error' ); f = () => 0; } f(1, 2); } catch (e2) { f_parsed = math.compile('4*y*(1-y)'); f = (y, t) => f_parsed.evaluate({ y: y, t: t }); let message; if (e2.toString().indexOf('not a function') != -1) { message = 'Parse error, consider using explicit multiplication - for example, write y*(t+1), instead of y(t+1).'; } else { message = e2; } d3.select(parse_warning).html(`${message}`); f = () => 0; } let aspect = 0.618; let xmin = -R; let xmax = R; let ymin = aspect * xmin; let ymax = aspect * xmax; let w = 0.9 * width; let h = aspect * w; let pad = 35; let xScale = d3 .scaleLinear() .domain([xmin, xmax]) .range([pad, w - pad]); let yScale = d3 .scaleLinear() .domain([ymin, ymax]) .range([h - pad, pad]); let pts_to_path = d3 .line() .x(function(d) { return xScale(d[0]); }) .y(function(d) { return yScale(d[1]); }); let dx = R / 10; let dy = dx; let grid = d3 .range(ymin + dy, ymax, dy) .map(y => d3.range(xmin + dx, xmax, dx).map(x => [x, y])); grid = grid.reduce(function(accumulated, currentValue) { return accumulated.concat(currentValue); }, []); let svg = d3 .create('svg') .attr('width', w) .attr('height', h) .style('background-color', 'white'); let slope_field = svg .append('g') .selectAll('path') .data(grid) .join('path') .attr('d', d => pts_to_path(segment(d[0], d[1], R / 30, f))) .style('stroke', 'black') .style('stroke-width', 0.8) .style('fill', 'none'); let tracking = true; svg .on('mousemove', function(evt) { if (tracking) { let [i, j] = d3.pointer(evt); inits.t0 = xScale.invert(i); inits.y0 = yScale.invert(j); draw_solution(i, j); } }) .on('mouseleave', function() { if (tracking) { svg.selectAll('.solution').remove(); inits.t0 = null; inits.y0 = null; } }) .on('click', function(evt) { if (tracking) { tracking = false; } else { tracking = true; let [i, j] = d3.pointer(evt); inits.t0 = xScale.invert(i); inits.y0 = yScale.invert(j); draw_solution(i, j); } }); function draw_solution(i, j) { let pts = rk4_1d(f, inits.y0, inits.t0, -R, R, 500); svg.selectAll('.solution').remove(); let solution_group = svg.append('g').attr('class', 'solution'); solution_group .append('path') .attr('class', 'solution') .attr('d', pts_to_path(pts)) .style('stroke', '#1f77b4') .style('stroke-width', 4) .style('fill', 'none'); solution_group .append('circle') .attr('cx', i) .attr('cy', j) .attr('r', 5) .attr('fill', '#2ca02c') .attr('stroke', 'black'); } // Define a segemnt of specified width 2r centered at (t,y) // and with slope f(t,y) to create the slope field. function segment(t, y, r, f) { let m = f(y, t); let b = r / Math.sqrt(m * m + 1); return [[t - b, y - m * b], [t + b, y + m * b]]; } svg .append('g') .attr('transform', `translate(0, ${h - pad})`) .call(d3.axisBottom(xScale)); svg .append('g') .attr('transform', `translate(${pad})`) .call(d3.axisLeft(yScale)); if (inits.t0) { // Not sure why this guy bugs out sometimes, but // the try/catch seems to make it bearable. try { draw_solution(xScale.invert(inits.t0), yScale.invert(inits.y0)); } catch (e) { ("pass"); } } return svg.node(); }32:T5bb,md`## Geometry The general, first order, ordinary differential equation (or ODE) can be written ${tex.block`y' = f(y,t).`} In this notation, we assume that ${tex`y`} is a function of ${tex`t`} and a solution of the ODE is a function ${tex`y(t)`} that makes the equation true. For example, if ${tex`f(y,t)=t\,y`}, then our differential equation is ${tex.block`y' = t\,y`} and a solution is ${tex`y(t)=e^{t^2/2}`} since, for then, ${tex.block`y'(t) = t\,e^{t^2/2} = t\,y(t).`} Thinking about differential equations from this symbolic perspective is certainly useful but limited since there a lots of important differential equations with no useful symbolic solution. In this notebook, we're going to take a look quick look at geometric and numeric approaches for approaching first order ODEs. The simplest tool explore a first order ODE is the *slope field*. Suppose we're curious about how the graph of a solution ${tex`y(t)`} looks as it passes through a rectangle ${tex.block`a\leq t \leq b, \: \: c \leq y \leq d`} in the ${tex`ty`}-plane. We simply draw a grid of points in that rectangle and, for each point, we draw a small line segment through that point with slope ${tex`f(y,t)`}. If ${tex`y(t)`} is a solution, then ${tex`y'(t)=f(y,t)`}, so that ${tex`f(y_0,t_0)`} tells us the slope of the graph of the graph of ${tex`y`}, as it passes through the point ${tex`(t_0,y_0)`}. Thus, to draw a graph of a soultion, all we have to do is follow the slope field!`33:T915,md`## Numerics Numerical techniques are generally used to solve *initial value problems* (or IVPs). Thus we are given a first order ODE together with an initial condition of the form ${tex`y(t_0) = y_0`}. For example, we might have the IVP ${tex.block` \begin{aligned} y' &= y, \\ y(0) &= 2. \end{aligned}`} Now, any function of the form ${tex`y(t)=c\,e^t`} solves the differential equation; this is called the *general solution*. Since the IVP also requires that ${tex`y(0)=2`}, we must also have ${tex.block` y(0) = c\,e^0 = c = 2.`} Thus, the solution to this specific IVP is ${tex`y(t)=2\,e^t`}. I typically speak of *a* solution to an ODE (because there can be many) but *the* solution to an IVP (because solutions are unique under reasonable assumptions). When you hover over the interactive tool above, your mouse essentially specifies an initial condition and the tool draws the solution to the corresponding IVP obtained via a numerical method. ### Euler's method Numerical techniques to solve IVPs don't find exact symbolic solutions; rather, they find approximate solutions represented as a list of decimal approximations to points on the graph of the solution. To be specific, suppose that we wish to estimate the solution of ${tex`y'=f(y,t)`} over the interval ${tex`a\leq t \leq b`} where ${tex`f(a)=y_0`}. We break the interval up into ${tex`n`} pieces: ${tex.block`a=t_0 < t_1 < t_2 < \cdots < t_n=b.`} Let's suppose that each piece has the same length, say ${tex`\Delta t`}. Thus, ${tex`t_{i+1}-t_i = \Delta t`} for each ${tex`i=0,1,\ldots,n-1`}. Now, we already know one point on the graph: ${tex`(t_0,y_0)`}. Next, to estimate the value of ${tex`y`} at ${tex`t_1`}, we can use the definition of the derivative: ${tex.block`y'(t_0) = f(y_0,t_0) \approx \frac{y(t_0+\Delta t) - y(t_0)}{h}`} We can solve this for ${tex`y(t_0+\Delta t)`} and use the fact that this is just another name for ${tex`y(t_1)`} to get ${tex.block`y(t_1) = y(t+\Delta t) \approx y(t_0) + f(y_1,t_1)\Delta t.`} Of course, once we have our estimate of ${tex`y(t_1)`}, we can use the same process to estimate ${tex`y(t_2)`} and, more generally, ${tex`y(t_{i+1})`} in terms of ${tex`y(t_i)`}. This provides an iterative way to find the list of points on our approximation. It's easy to implement this idea in Javascript: `34:T51d,md`Ultimately, though, Euler's method is a ${tex`1^{\text{st}}`} order method, i.e. the error in decreases linearly with the step size. Thus, if 10 steps gets us a single digit of accuracy, then 100 steps gives us two digits, and 1000 steps gives us three. Getting close to 16 digits (or machine precision) is not really obtainable so we've got to fundamentally improve the order of convergence. Here's another way to think of Euler's method that hints at a better way to improve it. Using the fundamental theorem of calculus, we can write: ${tex.block` \begin{aligned} \int_{t_i}^{t_{i+1}} f(y(t),t)\,dt &= \int_{t_i}^{t_{i+1}} y'(t)\,dt \\ &= y(t_{i+1}) - y(t_i). \end{aligned}`} Thus, ${tex.block`y_{i+1} = y_i + \int_{t_i}^{t_{i+1}} f(y(t),t)\,dt`} so approximating ${tex`y_{i+1}`} boils down to approximating that integral. The simplest approximation is ${tex.block`\int_{t_i}^{t_{i+1}} f(y(t),t)\,dt \approx y(t_i)\Delta t,`} which yields Euler's method. If we use a higher order approximation to the integral, we get a higher order IVP solver. If we use Simpson's rule, we obtain a fourth order technique called the *Runge-Kutta method* (or RK4). This is the technique used by the interactive tool above and is widely considered to be the first option to try. Here's an implementation of RK4: `35:[] 2f:{"title":"The geometry and numerics of first order ODEs","cells":[{"id":0,"value":"md`# The geometry and numerics of first order ODEs`","pinned":false,"hidden":false,"mode":"ojs"},{"id":604,"value":"intro_text = md`Write the general, first order, ODE in the form\n${tex.block`y' = f(y,t).`}\nWe can then draw the slope field for the ODE together with a solution to the ODE through that slope field. This interactive tool illustrates that process.\n\nEnter your expression involving \\`y\\` and \\`t\\` in the box below; the slope field should appear in the plot. Several examples are given in the drop down menu to illustrate. You can then hover over the plot to draw the solution through the point determined by your mouse. Click on the image to lock your solution and click again to unlock. You can also adjust the size of the viewing rectangle with the slider.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":376,"value":"viewof f1 = text({\n title: html`Right side of ODE: ${tex`\\large y'=\\,?`}`,\n value: f, //\"4y*(1-y)\",\n description: html`Enter an expression of ${tex`\\,y`} and ${tex`\\,t`}`,\n submit: 'Plot'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":999,"value":"viewof f = select({\n title: \"Examples\",\n options: [\n { label: \"Exponential\", value: 'y' },\n { label: \"Logistic\", value: '4y*(1-y)' },\n { label: \"Skew line limit\", value: 't-y' },\n { label: \"Sine\", value: 'sin(y)' },\n { label: \"Simple non-autonomous\", value: 't*y' },\n { label: \"Non-autonomous\", value: 'sin(t*y)' }\n ],\n value: '4y*(1-y)'\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":424,"value":"viewof R = slider({\n min: 1,\n max: 10,\n value: 2.5,\n title: \"Range of window\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":448,"value":"ODE_TeXed = tex.block`y' = ${math.parse(f1).toTex()}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":514,"value":"parse_warning = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":294,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1054,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":633,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":712,"value":"// Perform Euler's method with n steps over the interval [a,b]\n// for y' = f(y,t) with y(t0) = y0.\nfunction euler(f, y0, a, b, n) {\n let dt = (b - a) / n;\n let ti = a;\n let yi = y0;\n let pts = [[ti, yi]];\n for (let i = 0; i < n; i++) {\n yi = yi + f(yi, ti) * dt;\n ti = ti + dt;\n pts.push([ti, yi]);\n }\n return pts;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":718,"value":"md`Note that the code expresses ${tex`\\Delta t`} as \\`dt\\` and ${tex`y(t_i)`} as \\`yi\\`. After initializing ${tex`t_i`} and ${tex`y_i`} the code redefines them and pushes the new point onto the list of points. It finally returns that list.\n\nHere's how to use it to solve ${tex`y' = y`}:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":672,"value":"pts = euler((y, t) => y, 1, 0, 1, 10)","pinned":true,"hidden":false,"mode":"ojs"},{"id":824,"value":"md`The first argument, \\`(y, t) => y\\` defines ${tex`f(y,t) = y`} using Javascript's arrow notation for functions. The remaining arguments are simple numbers whose purpose can be read from the defintion of the \\`euler\\` function. Since the equation is ${tex`y' = y`} and the initial condition is ${tex`y(0) = 1`}, the exact solution is ${tex`y(t) = e^t`}. Here's a plot of that solution, together with the approximation obtained by from \\`euler\\`:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":738,"value":"{\n let plot = plotter({\n xDomain: [-0.1, 1.1],\n yDomain: [2.9, 0.9],\n width: 500,\n height: 300,\n grid: false\n });\n plot.func(x => math.exp(x), { stroke: '#1f77b4', width: 4 });\n for (let i = 0; i < pts.length - 1; i++) {\n plot.arrow(...pts[i], ...pts[i + 1]);\n plot.point(...pts[i], { fill: 'yellow' });\n }\n plot.point(...pts[pts.length - 1], { fill: 'yellow' });\n\n return plot.node;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":838,"value":"md`As we see, we get an underestimate, which makes perfect sense considering the “follow the tangent line” approach together with the concavity of the graph.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":848,"value":"md`### Improving Euler's method\n\nThe logical easy way to improve Euler's method is to simply use more steps for a smaller ${tex`\\Delta t`}. The next little interactive illustrates this:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":773,"value":"viewof n = slider({\n min: 1,\n max: 40,\n step: 1,\n value: 2,\n title: \"Number of steps\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":778,"value":"{\n let plot = plotter({\n xDomain: [-0.1, 1.1],\n yDomain: [2.9, 0.9],\n width: 500,\n height: 300,\n grid: false\n });\n let pts = euler(y => y, 1, 0, 1, n);\n plot.func(x => math.exp(x), { stroke: '#1f77b4', width: 4 });\n for (let i = 0; i < pts.length - 1; i++) {\n plot.arrow(...pts[i], ...pts[i + 1]);\n plot.point(...pts[i], { fill: 'yellow' });\n }\n plot.point(...pts[pts.length - 1], { fill: 'yellow' });\n\n return plot.node;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":855,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":103,"value":"// Apply the Runge-Kutta method with n steps over an interval\n// [a,b] containing t0.\nfunction rk4_1d(f, y0, t0, a, b, n) {\n let dt = (b - a) / n;\n if (t0 == a) { // Just go forward\n return rk4_multi_step_1d(f, y0, a, dt, n);\n } else if (b == t0) { // Just go backward\n let result = rk4_multi_step_1d(f, y0, b, -dt, n);\n result.reverse();\n return result;\n } else if (a < t0 && t0 < b) { // Go forward and backward and concatenate\n let i = Math.floor((t0 - a) / dt);\n let backward = rk4_multi_step_1d(f, y0, t0, -dt, i);\n backward.reverse();\n let forward = rk4_multi_step_1d(f, y0, t0, dt, n - i).slice(1);\n return backward.concat(forward);\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":216,"value":"function rk4_multi_step_1d(f, y0, t0, dt, n) {\n let result = [[t0, y0]];\n for (let i = 0; i < n; i++) {\n y0 = rk4_single_step_1d(f, y0, t0, dt);\n t0 = t0 + dt;\n if (Math.abs(y0) < 1000) {\n result.push([t0, y0]);\n } else {\n break;\n }\n }\n return result;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":957,"value":"function rk4_single_step_1d(f, y0, t0, dt) {\n let dt2 = dt / 2;\n let k1 = f(y0, t0);\n let k2 = f(y0 + k1 * dt2, t0 + dt2);\n let k3 = f(y0 + k2 * dt2, t0 + dt2);\n let k4 = f(y0 + k3 * dt, t0 + dt);\n\n return y0 + ((k1 + 2 * k2 + 2 * k3 + k4) * dt) / 6;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":947,"value":"md`### Comparison\nLet's compare \\`euler\\` and \\`rk4_1d\\`. To do so, we'll attack the simple IVP\n${tex.block`\n\\begin{aligned}\n y' &= y, \\\\\n y(0) &= 1.\n\\end{aligned}`}\nThe solution, of course, is ${tex`y(t)=e^t`} so that ${tex`y(1)=e`}. Thus, we can solve the IVP numerically and check the value at ${tex`t=1`} to get a sense of the quality of the technique. Note that, to an accuracy of 10 spots past the decmial, we have\n${tex.block`e \\approx 2.7182818285.`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":962,"value":"md`Euler's method with 10 steps:`","pinned":true,"hidden":false,"mode":"ojs"},{"id":939,"value":"euler_approx10 = euler(y => y, 1, 0, 1, 10).slice(-1)[0][1]","pinned":true,"hidden":false,"mode":"ojs"},{"id":967,"value":"md`Kinda, sorta in the ballpark, in that we've got 2(point)-something kinda close to seven. Here's Runge-Kutta with 10 steps:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":970,"value":"rk4_approx = rk4_1d(y => y, 1, 0, 0, 1, 10).slice(-1)[0][1]","pinned":true,"hidden":false,"mode":"ojs"},{"id":973,"value":"md`That looks like it's got lots of digits! Let's try Euler with 100000 steps:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":952,"value":"euler_approx100000 = euler(y => y, 1, 0, 1, 100000).slice(-1)[0][1]","pinned":true,"hidden":false,"mode":"ojs"},{"id":978,"value":"md`That looks like similar quality to the Runge-Kutta method with only 10 steps. Let's check to see which one is actually closer to ${tex`e`}:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":960,"value":"[Math.E - euler_approx100000, Math.E - rk4_approx]","pinned":true,"hidden":false,"mode":"ojs"},{"id":981,"value":"md`In fact, it looks like Runge-Kutta with only 10 steps is better than Euler than 100000 steps! That's the power of thinking about your numerical techniques!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1010,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":572,"value":"// d3 = require('d3@6')\nd3 = require('d3-selection@2', 'd3-scale@3', 'd3-shape@2', 'd3-array@2', 'd3-axis@2')","pinned":true,"hidden":false,"mode":"ojs"},{"id":575,"value":"math = require('mathjs@7.6.0')","pinned":true,"hidden":false,"mode":"ojs"},{"id":377,"value":"import { slider, text, select } from \"@jashkenas/inputs\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":735,"value":"import { plotter } from '@kjerandp/plotter'","pinned":true,"hidden":false,"mode":"ojs"},{"id":1078,"value":"md`### One global object`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1075,"value":"// Store the values of t0 and y0 in this external object so\n// they don't get rewritten when we resize R.\ninits = ({\n t0: null,\n y0: null\n})","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}