30:T488,## Problem Statement We wish to find a numerical solution to the first-order ordinary differential equation: ${tex.block`\frac{dy}{dt} = f(t, y)`} that passes through the initial condition ${tex`(t_0, y_0)`}. The examples below use ${tex.block` f(t, y) = y^2 - t `} which is a nonlinear, non-autonomous [Riccati equation](https://en.wikipedia.org/wiki/Riccati_equation) that has no elementary closed-form solutions. ## Problem Visualization The function ${tex`f(t,y)`} contains all of the information about the possible solutions passing through any point ${tex`(t_0, y_0)`}. Although it returns a scalar value at each ${tex`(t,y)`}, it is helpful to visualize this value as a vector. How does this work? Draw a small right-angle triangle with sides ${tex`\Delta t`} and ${tex`\Delta y`} that satisfy ${tex.block` \frac{\Delta y}{\Delta t} = f(t, y)`} then the hypotenuse at ${tex`(t,y)`} is our vector, and is tangent to the solution passing through this point. The plot below shows these triangles: we can visualize the problem of solving the ODE as simply connecting together the small tangent lines passing through the initial point ${tex`(t_0,y_0)`}.31:T474,## Numerical Methods for a Single Time Step The Euler method uses: ${tex.block`y_{i+1} = y_i + k_1`} The RK2 method uses: ${tex.block`y_{i+1} = y_i + k_2`} The RK4 method uses the weighted average: ${tex.block`y_{i+1} = y_i + \frac{1}{6}\left( k_1 + 2 k_2 + 2 k_3 + k_4\right)`} with ${tex.block` \begin{aligned} k_1 &= h\, f\left(t_i, y_i\right) \\ k_2 &= h\, f\left(t_i + \frac{h}{2}, y_i + \frac{k_1}{2}\right) \\ k_3 &= h\, f\left(t_i + \frac{h}{2}, y_i + \frac{k_2}{2}\right) \\ k_4 &= h\, f\left(t_i + h, y_i + k_3\right) \end{aligned} \; . `} This table summarizes how the 3 different methods build their solutions from the 4 ${tex`k_i`} estimates: Method | ${tex`k_1`} | ${tex`k_2`} | ${tex`k_3`} | ${tex`k_4`} -------|----|----|----|---- |${tex`t_i`} |${tex`t_i+h/2`} |${tex`t_i+h/2`} |${tex`t_i+h`} Euler | 1 | 0 | 0 | 0 RK2 | 0 | 1 | 0 | 0 RK4 |1/6 |2/6 |2/6 |1/6 Use the controls below to visualize how each ${tex`k_i`} is calculated from the tangent lines. Since ${tex`k_{i+1}`} is based on ${tex`k_i`}, it is helpful to select adjacent pairs to see how this works (e.g. ${tex`k_1`} and ${tex`k_2`}).32:T1210,function odePlot(params) { const {t0,y0,dt,viz} = params; const {k1,k2,k3,k4} = numericSolutions; const w = width; const h = Math.round(w / 1.618); const arrowLength = 0.8 * (tgrid[ngrid-1] - tgrid[0]) / ngrid; const arrows = ytgrid.map(({t,y,yp}) => { const s = Math.sqrt(1 + yp**2); const dx = arrowLength / s; const dy = dx * yp; const x1 = t - 0.5*dx, y1 = y - 0.5*dy, x2 = t + 0.5*dx, y2 = y + 0.5*dy; return { x1, x2, y1, y2 }; }); const marks = [ Plot.link(arrows, {x1:"x1", y1:"y1", x2:"x2", y2:"y2", stroke:"white"}), ]; if(viz.includes("triangles")) marks.push( Plot.link(arrows, {x1:"x1", y1:"y1", x2:"x2", y2:"y1", stroke:"white", strokeOpacity:0.5 }), Plot.link(arrows, {x1:"x2", y1:"y1", x2:"x2", y2:"y2", stroke:"white", strokeOpacity:0.5 }), ); if(viz.includes("step")) marks.push( Plot.ruleX([params], { x:d=>d.t0, stroke:"gray", strokeWidth:4, strokeOpacity:0.25 }), Plot.ruleX([params], { x:d=>d.t0+d.dt, stroke:"gray", strokeWidth:4, strokeOpacity:0.25 }), Plot.dot([params], { x:d=>d.t0, y:d=>d.y0, fill:"black", r:3 }), ); if(viz.includes("exact")) marks.push( Plot.line(getSmoothSolution(t0,y0), { x:"t", y:"y", stroke:"black", clip:"frame" }), ); if(viz.includes("solutions")) marks.push( solutionsGrid.map(soln => Plot.line(soln, { x:"t", y:"y", stroke:"white", strokeWidth:3, strokeOpacity:0.5, clip:"frame" })) ); let c = "red", sw = 2; if(viz.includes("k1")) marks.push( Plot.link([{x1:t0,x2:t0+dt,y1:y0,y2:y0+k1}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:sw}), Plot.link([{x1:t0,x2:t0+dt,y1:y0,y2:y0+k1}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:12, strokeOpacity:0, title:d=>"k1"}), Plot.dot([{x:t0,y:y0}], {x:"x", y:"y", stroke:c, r:6, strokeWidth:sw, fill:c, fillOpacity:0, title:d=>"f(t0,y0)"}), Plot.dot([{x:t0+dt,y:y0+k1}], {x:"x", y:"y", fill:c, r:3}), ); c = "green"; if(viz.includes("k2")) marks.push( //Plot.areaY([{x:t0, y1:y0+k1/2-k2/2, y2:y0}, {x:t0+dt, y1:y0+k1/2+k2/2, y2:y0+k2}], // {x:"x", y1:"y1", y2:"y2", fill:c, opacity:0.2}), Plot.link([{x1:t0, y1:y0+k1/2-k2/2, x2:t0+dt, y2:y0+k1/2+k2/2}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:2, strokeDasharray:[2,4]}), Plot.link([{x1:t0, y1:y0, x2:t0+dt, y2:y0+k2}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:2}), Plot.dot([{x:t0+dt/2,y:y0+k1/2}], {x:"x", y:"y", stroke:c, r:6, strokeWidth:2, fill:c, fillOpacity:0, title:d=>"f(t0+dt/2,y0+k1/2)"}), Plot.dot([{x:t0+dt,y:y0+k2}], {x:"x", y:"y", fill:c, r:3}), ); c = "magenta"; if(viz.includes("k3")) marks.push( //Plot.areaY([{x:t0, y1:y0+k2/2-k3/2, y2:y0}, {x:t0+dt, y1:y0+k2/2+k3/2, y2:y0+k3}], // {x:"x", y1:"y1", y2:"y2", fill:c, opacity:0.2}), Plot.link([{x1:t0, x2:t0+dt, y1:y0+k2/2-k3/2, y2:y0+k2/2+k3/2}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:2, strokeDasharray:[2,4]}), Plot.link([{x1:t0, x2:t0+dt, y1:y0, y2:y0+k3}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:2}), Plot.dot([{x:t0+dt/2,y:y0+k2/2}], {x:"x", y:"y", stroke:c, r:6, strokeWidth:2}), Plot.dot([{x:t0+dt,y:y0+k3}], {x:"x", y:"y", fill:c, r:3}), ); c = "orange", sw = 2; if(viz.includes("k4")) marks.push( //Plot.areaY([{x:t0, y1:y0+k3-k4, y2:y0}, {x:t0+dt, y1:y0+k3, y2:y0+k4}], // {x:"x", y1:"y1", y2:"y2", fill:c, opacity:0.2}), Plot.link([{x1:t0, x2:t0+dt, y1:y0+k3-k4, y2:y0+k3}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:sw, strokeDasharray:[2,4]}), Plot.link([{x1:t0, x2:t0+dt, y1:y0, y2:y0+k4}], {x1:"x1", x2:"x2", y1:"y1", y2:"y2", stroke:c, strokeWidth:sw}), Plot.dot([{x:t0+dt,y:y0+k3}], {x:"x", y:"y", stroke:c, r:6, strokeWidth:sw}), Plot.dot([{x:t0+dt,y:y0+k4}], {x:"x", y:"y", fill:c, r:3}), ); if(["k1","k2","k3","k4"].every(val => viz.includes(val))) { const yRK4 = y0 + (k1 + 2*k2 + 2*k3 + k4) / 6; marks.push( Plot.dot([{x:t0+dt,y:yRK4}], {x:"x", y:"y", stroke:"black", r:8, symbol:"square"}), Plot.text(["RK4"], { x:t0+dt, dx:24, y:yRK4, text:d=>d, fontSize:14 }), ); } return Plot.plot({ x: { domain:[tgrid[0], tgrid[tgrid.length-1]], label:"t" }, y: { domain:[ygrid[0], ygrid[ygrid.length-1]], label:"y" }, width: w, height: h, style: { backgroundColor: "lightblue" }, marks }); }33:[] 2e:{"title":"Numerical Solution of Ordinary Differential Equations","cells":[{"id":0,"value":"# Numerical Solution of Ordinary Differential Equations","pinned":false,"hidden":false,"mode":"md"},{"id":292,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"id":385,"value":"viewof viz1 = Inputs.checkbox([\"triangles\", \"solutions\"], { value:[\"triangles\"] })","pinned":false,"hidden":false,"mode":"ojs"},{"id":376,"value":"trianglesPlot = odePlot(Object.assign({viz:viz1}, defaults))","pinned":false,"hidden":false,"mode":"ojs"},{"id":437,"value":"## Divide and Conquer\n\nWe break the solution into small time steps, indexed by ${tex`i`}, of size ${tex`h`}, so\n${tex.block`t_i = t_0 + i h`}\n\nUse the controls below to visualize a single time step. Remember that the \"cost\" of our algorithm is the number of times we evaluate ${tex`f(t,y)`} per step times the number of steps. Our accuracy for each step will generally improve when the exact solution is \"smother\" over the step (the ideal would be a straight line). Using a larger number of smaller steps generally increases the accuracy but also the cost.","pinned":false,"hidden":false,"mode":"md"},{"id":442,"value":"viewof stepParams = Inputs.form({\n t0: Inputs.range([-1.5,1], {label:\"Initial Time\", value:-0.6, step: 0.01}),\n y0: Inputs.range([-1,1], { label:\"Initial Position\", value:0.3, step: 0.01}),\n dt: Inputs.range([0.1,2], {label: \"Step Size (s)\", value:2.0, step: 0.01}),\n viz: Inputs.input([\"step\",\"exact\"]),\n},\n{\n template: (I) => html`\n

${I['t0']}

${I['dt']}

${I['y0']}

`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":439,"value":"stepPlot = odePlot(stepParams)","pinned":false,"hidden":false,"mode":"ojs"},{"id":373,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":17,"value":"viewof params = Inputs.form({\n t0: Inputs.range([-1.5,1], {label:\"Initial Time\", value:-0.6, step: 0.01}),\n y0: Inputs.range([-1,1], { label:\"Initial Position\", value:0.3, step: 0.01}),\n dt: Inputs.range([0.1,2], {label: \"Step Size (s)\", value:2.0, step: 0.01}),\n ksel: Inputs.checkbox([\"k1\", \"k2\", \"k3\", \"k4\"], { value:[\"step\"] }),\n},\n{\n template: (I) => html`\n

${I['ksel']}

${I['t0']}

${I['dt']}

${I['y0']}

`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":394,"value":"methodsPlot = {\n const viz = [\"step\",\"exact\"].concat(params.ksel);\n return odePlot(Object.assign({viz}, params));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":362,"value":"---","pinned":false,"hidden":false,"mode":"md"},{"id":380,"value":"defaults = ({ t0: -0.6, y0: 0.3, dt: 2.0 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":54,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":166,"value":"numericSolutions = {\n const { t0, y0, dt } = params;\n const k1 = dt * dydt(t0, y0);\n const k2 = dt * dydt(t0 + dt/2, y0 + k1/2);\n const k3 = dt * dydt(t0 + dt/2, y0 + k2/2);\n const k4 = dt * dydt(t0 + dt, y0 + k3);\n return { k1, k2, k3, k4 };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":409,"value":"solutionsGrid = {\n const t0 = 1.0;\n const ymin = ygrid[0], ymax = ygrid[ngrid-1];\n const n = 15, dy = (ymax - ymin) / (n + 1);\n return d3.range(n).map(k => {\n const y0 = ymin + (k + 1) * dy;\n return getSmoothSolution(t0, y0);\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":104,"value":"function getSmoothSolution(t0, y0) {\n const tmin = tgrid[0], tmax = tgrid[ngrid-1], deltat = 0.001;\n const ymin = ygrid[0], ymax = ygrid[ngrid-1];\n const np = Math.ceil((tmax - t0) / 0.001);\n const hp = (tmax - t0) / (np - 1);\n let t = t0, y = y0;\n const solnp = [ { t, y } ];\n for(let i=1; i < np; ++i) {\n y = RK4step(dydt, hp, t, y);\n t = t0 + i * hp;\n solnp.push({ t, y });\n if(y < ymin || y > ymax) break;\n }\n const nm = Math.ceil((t0 - tmin) / 0.001);\n const hm = (t0 - tmin) / (nm - 1);\n t = t0, y = y0;\n const solnm = [ { t, y } ];\n for(let i=1; i < nm; ++i) {\n y = RK4step(dydt, -hm, t, y);\n t = t0 - i * hm;\n solnm.push({ t, y });\n if(y < ymin || y > ymax) break;\n }\n return solnm.slice(1).reverse().concat(solnp);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":102,"value":"function RK4step(f, h, ti, yi) {\n const k1 = h * f(ti, yi);\n const k2 = h * f(ti + h/2, yi + k1/2);\n const k3 = h * f(ti + h/2, yi + k2/2);\n const k4 = h * f(ti + h, yi + k3);\n return yi + (k1 + 2 * k2 + 2 * k3 + k4) / 6;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":41,"value":"ytgrid = {\n const {atan,PI} = Math;\n const tstep = 4 / ngrid;\n return d3.cross(tgrid, ygrid).map(([t,y]) => {\n const yp = dydt(t, y);\n const rot = 90 - 180/PI * atan(yp);\n return { t, y, yp, rot };\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"ygrid = {\n const y1 = -2, y2 = 2;\n return d3.range(ngrid).map(k => y1 + k * (y2 - y1) / (ngrid - 1));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"tgrid = {\n const t1 = -1.5, t2 = 2.5;\n return d3.range(ngrid).map(k => t1 + k * (t2 - t1) / (ngrid - 1));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"ngrid = 50","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"function dydt(t, y) {\n return y * y - t;\n}","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}