31:T48d,The Shanks transformation doesn't magically make divergent series converge, and of course *divergent series don't converge*! This extraneous "correctness" (by the yardstick of more sophisticated methods) happens because assumed *convergence* but didn't actually enforce whether ${tex`A`} was the limit to which it was *converging* or from which it was *diverging*. Analytic continuation provides more rigorous tools with which to explain these series, as explained by Mathologer in his video [Numberphile v. Math: the truth about 1+2+3+...=-1/12](https://www.youtube.com/watch?v=YuIIjLr6vUA). The eponymous Daniel Shanks also explores these specific issues with more precision and rigor in his 1955 paper [*Non-linear Transformations of Divergent and Slowly Convergent Sequences*](https://onlinelibrary.wiley.com/doi/abs/10.1002/sapm19553411). It's a moderately approachable paper for the non-mathematician, though certainly somewhat involved. I'll leave this exploration here for now, though I'd highly recommend [digging up a copy of the paper](https://en.wikipedia.org/wiki/Sci-Hub) and at the very least browsing through some of the more entertaining results!32:T154a,function makePlot( series, { xDomain = [0, 20], yDomain = [1e-9, 10], showIndex = 0, w = width, h = Math.max(500, Math.floor(width * 0.7)), caption } = {} ) { const plot = Plot.plot({ width: w, height: h, grid: true, x: { domain: xDomain, label: "n" }, y: { domain: yDomain, type: "log", label: "|A(n) - π|" } }); const svg = d3.select(plot); const margin = { t: 20, r: 10, b: 20, l: 40 }; const xScale = plot.scale("x"); const yScale = plot.scale("y"); const line = d3 .line() .x(({ n }) => xScale.apply(n)) .y(({ absError }) => yScale.apply(absError)); const gArrows = svg.append("g").attr("class", "arrows"); const gLines = svg .append("g") .attr("class", "lines") .selectAll("path.line") .data(series) .join((enter) => enter .append("path") .attr("class", "line") .attr("d", (c) => line(c)) .attr("fill", "none") .attr("stroke", "#888") .attr("stroke-dasharray", "2, 2") .attr("stroke-width", 1) ); const gPoints = svg .append("g") .attr("class", "points") .selectAll("g.trace") .data(series) .join((enter) => enter.append("g").attr("class", "trace")) .selectAll("g.circle") .data((s) => s) //.filter(({ n }) => n > 0)) .join((enter) => { const g = enter .append("g") .attr("class", "circle") .attr( "transform", ({ n, absError }) => `translate(${xScale.apply(n)},${yScale.apply(absError)})` ); g.append("circle") .attr("r", 3) .style("fill", (d) => d3.schemeSet1[d.seriesIdx]) .style("stroke", "black") .style("stroke-width", 1); g.append("circle") .attr("r", 20) .style("fill", "transparent") .style("cursor", "pointer") .on("click", hover) .on("mouseover", hover); }); function hover(event, { seriesIdx = null, n } = {}) { if (seriesIdx === null) { gArrows.selectAll("*").remove(); return; } const arrows = []; for (let s = 0; s < seriesIdx; s++) { for (let m = n - seriesIdx; m < n + seriesIdx - s * 2 - 1; m++) { arrows.push([series[s][m], series[s + 1][m]]); arrows.push([series[s][m + 1], series[s + 1][m]]); arrows.push([series[s][m + 2], series[s + 1][m]]); } } const hierarchy = []; for (let s = 0; s <= seriesIdx; s++) { const h = []; for (let m = n - seriesIdx; m <= n + seriesIdx - s * 2; m++) { h.push(series[s][m]); } hierarchy.push(h); } figure.value = hierarchy; figure.dispatchEvent(new CustomEvent("input")); gArrows .selectAll("g.arrows") .data(arrows) .join( (enter) => enter.append("g").attr("class", "arrows").call(arrow), (update) => update.call(arrow), (remove) => remove.remove() ); } function makeArrow(el) { function computeArrow([d0, d1]) { const x1 = xScale.apply(d0.n); const x2 = xScale.apply(d1.n); const y1 = yScale.apply(d0.absError); const y2 = yScale.apply(d1.absError); const rx = x2 - x1; const ry = y2 - y1; const l = Math.hypot(rx, ry); const tx = rx / l; const ty = ry / l; const nx = ty; const ny = -tx; const inset = 10; const al = 7 + inset; const aw = 2.5; const hw = 0.75; return ( d3.line()([ [x1 + hw * nx + tx * inset, y1 + hw * ny + ty * inset], [x2 + hw * nx - tx * al, y2 + hw * ny - ty * al], [x2 + aw * nx - tx * al, y2 + aw * ny - ty * al], [x2 - tx * inset, y2 - ty * inset], [x2 - aw * nx - tx * al, y2 - aw * ny - ty * al], [x2 - hw * nx - tx * al, y2 - hw * ny - ty * al], [x1 - hw * nx + tx * inset, y1 - hw * ny + ty * inset] ]) + "Z" ); } return el .attr("d", (d) => computeArrow(d)) .style("fill", ([{ seriesIdx }]) => d3.schemeSet1[seriesIdx + 1]) .style("stroke", "black") .style("stroke-width", 0.1); } function arrow(el) { const data = el.data(); if (!data.length) return; el.selectAll("path") .data(data) .join( (enter) => enter.append("path").call(makeArrow), (update) => update.call(makeArrow), (remove) => remove.remove() ); } const gLabels = svg.append("g").attr("class", "labels"); function label(i) { if (i === 0) return "An"; const script = "⁰¹²³⁴⁵⁶⁷⁸⁹"; return `S${i === 1 ? "" : script[i]}(An)`; } gLabels .selectAll("text") .data(series) .join((enter) => enter.append("text")) .text((d) => label(d[0].seriesIdx)) .attr("font-size", 16) .attr("font-family", "serif") .attr("font-style", "italic") .attr("text-anchor", "middle") .attr("dy", ([{ seriesIdx }]) => (seriesIdx > 3 ? 20 : -10)) .attr("x", (d) => xScale.apply( d[Math.max(0, Math.min(12 - d[0].seriesIdx, d.length - 1))].n ) ) .attr("y", (d) => yScale.apply( d[Math.max(0, Math.min(12 - d[0].seriesIdx, d.length - 1))].absError ) ); const figure = htl.html`

${svg.node()}${ caption ? htl.html`

${caption}

` : "" }`; //hover(null, { seriesIdx: series.length - 1, n: showIndex }); return figure; }33:[] 2f:{"title":"Shanks Transformation","cells":[{"id":0,"value":"# Shanks Transformation","pinned":false,"hidden":false,"mode":"md"},{"id":2434,"value":"*[Read this notebook ➡️](https://rreusser.github.io/notebooks/shanks-transformation/) on my personal website via [Observable Notebook Kit](https://observablehq.com/notebook-kit/kit).*","pinned":false,"hidden":false,"mode":"md"},{"id":855,"value":"I've encountered some lovely mathematical results in my time, but the *[Shanks Transformation](https://en.wikipedia.org/wiki/Shanks_transformation)* for nonlinear series acceleration is truly one of the most delightful.","pinned":false,"hidden":false,"mode":"md"},{"id":859,"value":"Let's start with a very slowly converging series for ${tex`\\pi`} based on the Taylor series of the arctangent function,\n${tex.block`\n\\begin{aligned}\n\\pi &= 4 \\arctan(1) \\\\\n\\pi &= 4 \\left( 1 - \\frac{1}{3} + \\frac{1}{5} - \\frac{1}{7} + \\cdots \\right)\n\\end{aligned}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":1073,"value":"We define partial sums of this series as the function ${tex`A(n),`}\n${tex.block`\nA_n = 4 \\sum_{k = 0}^{n} \\frac{(-1)^k}{2k + 1}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":2249,"value":"I'll define a little helper for ${tex`\\sum_{k=0}^n`} since we'll use it repeatedly.","pinned":false,"hidden":false,"mode":"md"},{"id":1920,"value":"sum = (f) => (n) => {\n let sum = 0;\n for (let k = 0; k <= n; k++) sum += f(k);\n return sum;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1911,"value":"A = sum((k) => (4 * (k % 2 ? -1 : 1)) / (2 * k + 1))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1078,"value":"How slowly does this series converge? Very, very slowly!","pinned":false,"hidden":false,"mode":"md"},{"id":1014,"value":"| ${tex`n`} | ${tex`A_n`} |\n| -------- | ------- |\n| ${tex`0`} | ${tex`${piToTex(A(0))}`} |\n| ${tex`1`} | ${tex`${piToTex(A(1))}`} |\n| ${tex`2`} | ${tex`${piToTex(A(2))}`} |\n| ${tex`{10}`} | ${tex`${piToTex(A(10))}`} |\n| ${tex`{100}`} | ${tex`${piToTex(A(100))}`} |\n| ${tex`{1000}`} | ${tex`${piToTex(A(1000))}`} |\n| ${tex`{10000}`} | ${tex`${piToTex(A(10000))}`} |\n| ${tex`{100000}`} | ${tex`${piToTex(A(100000))}`} |\n| ${tex`{1000000}`} | ${tex`${piToTex(A(1000000))}`} |","pinned":false,"hidden":false,"mode":"md"},{"id":1090,"value":"Achieving error on the order of ${tex`10^{-6}`} requires about ${tex`10^6`} terms. Until I learned about Shanks transformation, I simply accepted that not every series is destined to be useful.","pinned":false,"hidden":false,"mode":"md"},{"id":1397,"value":"## A first attempt","pinned":false,"hidden":false,"mode":"md"},{"id":1151,"value":"Upon closer inspection, this is a silly conclusion. Even cursory visual inspection shows that any two consecutive terms bracket the correct answer.","pinned":false,"hidden":false,"mode":"md"},{"id":1101,"value":"Plot.plot({\n grid: true,\n width: Math.min(640, width),\n height: 200,\n y: { label: \"A(n)\" },\n x: { label: \"n\" },\n marks: [\n Plot.line(tabulate(A, [0, 20]), {\n x: \"n\",\n y: \"An\",\n stroke: \"#888\",\n strokeDasharray: \"1, 3\",\n strokeWidth: 1.25\n }),\n Plot.dot(tabulate(A, [0, 20]), {\n x: \"n\",\n y: \"An\",\n stroke: \"black\",\n strokeWidth: 1.5,\n fill: d3.schemeSet1[0]\n }),\n Plot.ruleY([Math.PI], {}),\n Plot.text([{ x: 19.5, y: Math.PI, text: \"π\" }], {\n x: \"x\",\n y: \"y\",\n text: \"text\",\n dy: -12,\n fontSize: 15\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1192,"value":"What happens if we just average two consecutive terms?","pinned":false,"hidden":false,"mode":"md"},{"id":1158,"value":"| ${tex`n`} | ${tex`A_n`} | ${tex`\\frac{1}{2}(A_n + A_{n+1})`} |\n| -------- | ------- | ------- |\n| ${tex`0`} | ${tex`${piToTex(A(0))}`} | ${tex`${piToTex((A(0) + A(1))/2)}`} |\n| ${tex`1`} | ${tex`${piToTex(A(1))}`} | ${tex`${piToTex((A(1) + A(2))/2)}`} |\n| ${tex`2`} | ${tex`${piToTex(A(2))}`} | ${tex`${piToTex((A(2) + A(3))/2)}`} |\n| ${tex`10`} | ${tex`${piToTex(A(10))}`} | ${tex`${piToTex((A(10) + A(11))/2)}`} |\n| ${tex`100`} | ${tex`${piToTex(A(100))}`} | ${tex`${piToTex((A(100) + A(101))/2)}`} |\n| ${tex`1000`} | ${tex`${piToTex(A(1000))}`} | ${tex`${piToTex((A(1000) + A(1001))/2)}`} |","pinned":false,"hidden":false,"mode":"md"},{"id":1419,"value":"At a thousand terms, we go from three digits of accuracy to six, so it would seem we're onto something! In general we get about twice as many digits out of the average of two consecutive terms, but it still takes a *lot* of terms to get a decent answer.","pinned":false,"hidden":false,"mode":"md"},{"id":1422,"value":"## Shanks Transformation\n","pinned":false,"hidden":false,"mode":"md"},{"id":1413,"value":"The essence of Shanks Transformation is to use consecutive terms to extrapolate as far as we can toward the right answer, only with a better choice than our averaging above. So instead we use three consecutive terms to fit a geometric sequence of three parameters, then solve for its limit.","pinned":false,"hidden":false,"mode":"md"},{"id":1428,"value":"To get the ball rolling, we write our series in the form of a general geometric sequence,\n\n${tex.block`\n\\begin{aligned}\nA_{n-1} &= A + \\alpha \\cdot q^{n-1} \\\\\nA_n &= A + \\alpha \\cdot q^{n} \\\\\nA_{n+1} &= A + \\alpha \\cdot q^{n+1}.\n\\end{aligned}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":1513,"value":"Solving the first two equations for ${tex`q`} yields\n\n${tex.block`\nq = \\frac{A_n - A}{A_{n - 1} - A}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":1516,"value":"Performing the same operation on the second and third equations, we have\n\n${tex.block`\nq = \\frac{A_n - A}{A_{n - 1} - A} = \\frac{A_{n + 1} - A}{A_n - A}.\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":1521,"value":"Cross-multiplying and solving for ${tex`A`} gives us\n\n${tex.block`A = \\frac{A_{n+1}A_{n-1} - A_n^2}{A_{n+1} - 2A_n + A_{n - 1}}.`}","pinned":false,"hidden":false,"mode":"md"},{"id":1532,"value":"We could backtrack and solve for ${tex`\\alpha`} and ${tex`q`}, but if the sequence is convergent then ${tex`|q| < 1,`} so that as ${tex`n \\rarr \\infty`}, our second equation for ${tex`A_n`} gives\n\n${tex.block`\n\\lim_{n \\rarr \\infty} A_n = A.\n`}\n\n${tex`A`} is in fact the limit we're looking for, so we define this as the *Shanks transformation* of series ${tex`A_n`} and denote it\n\n${tex.block`\nS(A_n) = \\frac{A_{n+1}A_{n-1} - A_n^2}{A_{n+1} - 2A_n + A_{n - 1}}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":1559,"value":"We can implement this pretty easily as a higher order JavaScript function (although please note that this is a very, very inefficient implementation for the sake of clarity, with many repeated, unnecessary function evaluations!)","pinned":false,"hidden":false,"mode":"md"},{"id":51,"value":"//S = (A) => (n) =>\n// (A(n + 1) * A(n - 1) - A(n) ** 2) / (A(n + 1) - 2 * A(n) + A(n - 1))\n\n// An exactly equivalent and more numerically stable form:\nS = (A) => (n) =>\n A(n + 1) - (A(n + 1) - A(n)) ** 2 / (A(n + 1) - 2 * A(n) + A(n - 1))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1579,"value":"The plot below shows the error of ${tex`A_n`} and the Shanks transformation ${tex`S(A_n)`}. Hover over various points and observe in the plot and table below how much the transformation accelerates convergence. Suddenly the twentieth term yields accuracy the original series takes nearly ten thousand terms to produce!","pinned":false,"hidden":false,"mode":"md"},{"id":129,"value":"tabulate = (A, [start, end], { expected = Math.PI, seriesIdx } = {}) =>\n range(start, end)\n .map((n) => ({ n, An: A(n) }))\n .map(({ n, An }) => ({\n seriesIdx,\n n,\n An,\n absError: Math.abs(An - expected)\n }))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1201,"value":"viewof firstShanks = makePlot(\n [\n tabulate(Shanks(A, 0), [0, 30], { seriesIdx: 0 }),\n tabulate(Shanks(A, 1), [1, 30], { seriesIdx: 1 })\n ],\n {\n yDomain: [1e-5, 4],\n showIndex: 4,\n w: Math.min(width, 640),\n h: Math.floor(0.7 * Math.min(width, 640)),\n caption:\n \"Hover over a point to observe the values it depends upon. Inspect corresponding values in the table below.\"\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1623,"value":"makeTable(firstShanks)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1437,"value":"## Iterated Shanks Transformation","pinned":false,"hidden":false,"mode":"md"},{"id":1562,"value":"But wait, you say! The Shanks transformation ${tex`S(A_n)`} gives us a *new* series which presumably converges to the same limit. Can we apply the Shanks transformation to the Shanks transformation and achieve an even better approximation?","pinned":false,"hidden":false,"mode":"md"},{"id":1564,"value":"Remarkably, the answer is *yes*! We can simply apply the Shanks transformation over and over, each time yielding a better and better approximation, each successive iteration requiring only two additional terms of ${tex`A_n`}.","pinned":false,"hidden":false,"mode":"md"},{"id":292,"value":"Shanks = (A, n = 1) => [...Array(n).keys()].reduce(S, A)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1728,"value":"Suddenly, we're at a pretty remarkable state of affairs. Consider the fifth iterate, ${tex`S^5(A_n).`} The very first term is already ${tex`${piToTex(Shanks(A, 4)(4))}`}, and really only depends upon the first nine terms of the original series, the most accurate of which was ${tex`${piToTex(A(10))}`}!","pinned":false,"hidden":false,"mode":"md"},{"id":2431,"value":"Let's be clear about what happened here. We took the first eleven terms of the original series. These values only contain one or two digits of accuracy. But they contain a lot of _information_. By plugging them over and over into the Shanks transformation, we extract that information and obtain _seven digits_ of accuracy!","pinned":false,"hidden":false,"mode":"md"},{"id":1815,"value":"Hover over results in the table below and observe the values in the table which follows, keeping in mind that the series itself only enters in the terms of ${tex`A_n`} itself and that the Shanks transformation introduces no additional information. I've plotted up to ${tex`S^5(A_n)`}, though I've truncated the series once they run into machine precision and yield the exactly correct answer.","pinned":false,"hidden":false,"mode":"md"},{"id":930,"value":"viewof iteratedShanks = makePlot(\n [\n tabulate(Shanks(A, 0), [0, 30], { seriesIdx: 0 }),\n tabulate(Shanks(A, 1), [1, 30], { seriesIdx: 1 }),\n tabulate(Shanks(A, 2), [2, 30], { seriesIdx: 2 }),\n tabulate(Shanks(A, 3), [3, 30], { seriesIdx: 3 }),\n tabulate(Shanks(A, 4), [4, 30], { seriesIdx: 4 }),\n tabulate(Shanks(A, 5), [5, 30], { seriesIdx: 5 }),\n tabulate(Shanks(A, 6), [6, 12], { seriesIdx: 6 }),\n tabulate(Shanks(A, 7), [7, 11], { seriesIdx: 7 })\n ],\n {\n h: Math.max(300, Math.floor(width * 0.8)),\n showIndex: 5,\n yDomain: [1e-16, 4],\n //xDomain: [0, 13],\n caption:\n \"Hover over a point to observe the values it depends upon. Inspect corresponding values in the table below.\"\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1655,"value":"makeTable(iteratedShanks)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1853,"value":"Since no new information enters the picture via the Shanks transformation, the stunning aspect to me is just *how much information* is contained even in the original, slowly converging, inaccurate series ${tex`A_n`}.","pinned":false,"hidden":false,"mode":"md"},{"id":1937,"value":"## Taking the transformation too far","pinned":false,"hidden":false,"mode":"md"},{"id":1939,"value":"Now we just get reckless and greedy. How about the [comically divergent series](https://en.wikipedia.org/wiki/1_%E2%88%92_2_%2B_3_%E2%88%92_4_%2B_%E2%8B%AF)\n\n${tex.block`\n1 - 2 + 3 - 4 - \\cdots = \\frac{1}{4}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":2209,"value":"B = sum((k) => k * Math.pow(-1, k - 1))","pinned":true,"hidden":false,"mode":"ojs"},{"id":2306,"value":"The series diverges, as is should.","pinned":false,"hidden":false,"mode":"md"},{"id":2267,"value":"B(10000)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1950,"value":"However, the first term of the fourth iterate ${tex`S^4(A_n)`} shows that the Shanks-transformed series converges just fine to ${tex`1/4`}.","pinned":false,"hidden":false,"mode":"md"},{"id":1860,"value":"Shanks(B, 4)(4)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2005,"value":"Or perhaps the similarly divergent series\n${tex.block`1 - 2 + 4 - 8 + \\cdots = \\frac{1}{3}`}","pinned":false,"hidden":false,"mode":"md"},{"id":2217,"value":"C = sum((k) => Math.pow(-2, k))","pinned":true,"hidden":false,"mode":"ojs"},{"id":2311,"value":"The divergent series obviously diverges.","pinned":false,"hidden":false,"mode":"md"},{"id":2270,"value":"C(100)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2009,"value":"Again though, a single iteration of the Shanks transformation ${tex`S(A_n)`} nails it.","pinned":false,"hidden":false,"mode":"md"},{"id":1966,"value":"Shanks(C)(1)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2200,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":1957,"value":"## Helpers","pinned":false,"hidden":false,"mode":"md"},{"id":427,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1595,"value":"function makeTable(series) {\n return html`\n \n \n ${range(0, series.length - 1).map(\n (i) =>\n htl.html``\n )}\n \n ${series[0].map(\n ({ n }, i) => htl.html`\n \n ${series.map(\n (_, j) => htl.html``\n )}\n `\n )}\n \n

${tex`n`}	${\n i === 0 ? tex`A_n` : tex`S^{${i === 1 ? \"\" : i}}(A_n)`\n }
${n}	\n ${tex`${piToTex(series[j][i - j] ? series[j][i - j].An : \"\")}`}\n

`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":163,"value":"range = (start, end) => [...Array(end - start + 1).keys()].map((n) => n + start)","pinned":false,"hidden":false,"mode":"ojs"},{"id":962,"value":"function piToTex(val) {\n if (typeof val === \"undefined\" || val === \"\") return \"\";\n const exact = Math.PI.toString();\n const obs = val.toFixed(15).padEnd(17, \"0\");\n const good = [];\n const bad = [];\n for (let i = 0; i < obs.length; i++) {\n const isGood = !bad.length && obs[i] === exact[i];\n (isGood ? good : bad).push(obs[i]);\n if (bad.length >= 2 && good.length + bad.length >= 3) break;\n }\n return `${good.join(\"\")}${bad.length ? `\\\\color{red}{${bad.join(\"\")}}` : \"\"}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1066,"value":"","pinned":false,"hidden":false,"mode":"html"}],"files":"$W33","stdlib":"1","titling":"auto"}