30:Tb39,plot = { const n = 1000, root3 = Math.sqrt(3); const x = d3.range(n+1).map(i => xmin + (i / n) * (xmax - xmin)); const fx = x.map(f); const bin = x.map((x,i) => i < n ? Math.floor((x - xmin) * intervals / (xmax - xmin)) : intervals-1); const odd = fx.map((y,i) => bin[i] % 2 == 1 ? y : NaN); const even = fx.map((y,i) => bin[i] % 2 == 0 ? y : NaN); const approx = x.map((x,i) => { const k = bin[i]; const xlo = xmin + k / intervals * (xmax - xmin); const xhi = xmin + (k + 1) / intervals * (xmax - xmin); if(rule == "midpoint") { const xmid = 0.5 * (xlo + xhi); return f(xmid); } if(rule == "trapezoid") { const s = (x - xlo) / (xhi - xlo); return f(xlo) * (1 - s) + f(xhi) * s; } if(rule == "Simpson") { const xmid = 0.5 * (xlo + xhi); return ( f(xlo) * (x - xhi) * (x - xmid) / (xlo - xhi) / (xlo - xmid) + f(xhi) * (x - xlo) * (x - xmid) / (xhi - xlo) / (xhi - xmid) + f(xmid) * (x - xlo) * (x - xhi) / (xmid - xlo) / (xmid - xhi)); } if(rule == "Gauss") { // 2-point Gauss Legendre on [xmin,xmax] const xmid = 0.5 * (xlo + xhi); const h = 0.5 * (xhi - xlo); const xm = xmid - h / root3, xp = xmid + h / root3; const fm = f(xm), fp = f(xp); return x < xmid ? fm : fp; } }); let points = [ ]; const xp = d3.range(intervals + 1).map(i => xmin + i * (xmax - xmin) / intervals); if(rule == "midpoint") { points = d3.range(intervals).map(i => ({ xpt:(xp[i]+xp[i+1])/2, ypt:f((xp[i]+xp[i+1])/2) })); } else if(rule == "trapezoid") { points = xp.map(x => ({ xpt:x, ypt: f(x) })); } else if(rule == "Simpson") { points = d3.range(intervals).map(i => ({ xpt:(xp[i]+xp[i+1])/2, ypt:f((xp[i]+xp[i+1])/2) })).concat( xp.map(x => ({ xpt:x, ypt: f(x) }))); } else if(rule == "Gauss") { const h = (xmax - xmin) / intervals; for(let k=0; k < intervals; ++k) { const xmid = 0.5 * (xp[k] + xp[k+1]); const x1 = xmid - 0.5 * h / root3, x2 = xmid + 0.5 * h / root3; points.push({ xpt: x1, ypt: f(x1) }); points.push({ xpt: x2, ypt: f(x2) }); } } const showExact = options.includes("exact"), showApprox = options.includes("approx"), showPoints = options.includes("points"); return Plot.plot({ width: width, marks: [ showExact ? Plot.areaY(odd, { x: (_,i) => x[i], y: y=>y, fill: "blue", fillOpacity: 0.15 }) : null, showExact ? Plot.areaY(even, { x: (_,i) => x[i], y: y=>y, fill: "blue", fillOpacity: 0.30 }) : null, //Plot.line(fx, { x: (_,i) => x[i], y: y=>y, stroke: "blue", strokeWidth: 2 }), showApprox ? Plot.line(approx, { x: (_,i) => x[i], y: y=>y, stroke: "red", strokeWidth: 4, strokeOpacity: 0.5 }) : null, showPoints ? Plot.dot(points, { x: "xpt", y:"ypt", fill: "red", r: 5 }) : null, ] }); }31:T559,results = { // Calculate brute force using trapezoid rule on a fine grid. const n = 1000, delta = (xmax - xmin) / n; const x = d3.range(n+1).map(i => xmin + i * delta); const exact = 0.5 * delta * (f(xmin) + f(xmax) + 2 * d3.sum(d3.range(1, n), i => f(x[i]))); const h = (xmax - xmin) / intervals; const xp = d3.range(intervals + 1).map(i => xmin + i * h); let approx = 0; for(let i=0; i <= intervals; ++i) { if(rule == "midpoint") { if(i < intervals) { const xmid = 0.5 * (xp[i] + xp[i+1]); approx += (xp[i+1] - xp[i]) * f(xmid); } } if(rule == "trapezoid") { const coef = (i % intervals) == 0 ? 1 : 2; approx += h / 2 * coef * f(xp[i]); } if(rule == "Simpson") { if(i < intervals) { const a = xp[i], b = xp[i+1]; approx += (b - a) / 6 * (f(a) + 4 * f((a + b)/2) + f(b)); } } if(rule == "Gauss") { if(i < intervals) { const a = xp[i], b = xp[i+1], x0 = 1 / Math.sqrt(3); approx += (b - a) / 2 * ( f((b + a + (b - a) * x0) / 2) + f((b + a - (b - a) * x0) / 2)); } } } const round = d3.format(".5f"); return html` The exact result is ${round(exact)}.
The approximate result with the ${rule} method is ${round(approx)}.`; }32:[] 2e:{"title":"Numerical Integration","cells":[{"id":0,"value":"intro = md`# Numerical Integration\n\nWe wish to evaluate the definite 1D integral\n${tex.block`\nI = \\int_a^b f(x) dx`}\nFor concreteness, we take ${tex`a = -1, b = +1`} (since we can always redefine ${tex`x`} to make this true). Next, pick one of the functions below, a quadrature rule, and the number of intervals to apply the rule on:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":109,"value":"viewof f = Inputs.radio(functions, {label: \"Function f(x)\", value: functions.get(\"cubic\")})","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"viewof rule = Inputs.radio(\n [\"midpoint\",\"trapezoid\",\"Simpson\",\"Gauss\"],\n {label: \"Quadrature rule\", value: \"midpoint\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"viewof intervals = Inputs.range([1, 9], {label: \"Number of intervals\", step: 1, value: 3})","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":70,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":133,"value":"viewof options = {\n const options = [\"exact\", \"approx\", \"points\"];\n return Inputs.checkbox(options, {label: \"Display options\", value: options});\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":93,"value":"For a closer look at the Gaussian quadrature rule, see [this notebook](https://observablehq.com/@dkirkby/gauss-quadrature).\n\n---","pinned":false,"hidden":false,"mode":"md"},{"id":112,"value":"functions = new Map([\n [\"linear\", x => x + 0.5],\n [\"quadratic\", x => (x + 0.4) * (x - 0.6)],\n [\"cubic\", x => 7 * x**3 - 8 * x**2 - 3 * x + 3],\n [\"sinc\", x => x != 0 ? Math.sin(20 * x) / (20 * x) : 1],\n [\"Runge\", x => 1 / (1 + 25 * x**2)],\n])","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"xmin = -1","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"xmax = +1","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W32","stdlib":"1","titling":"auto"}