(replay, draw({ f: (x) => Math.pow(x, 3) - 2 * x - 4, a: -3, b: 4 }))
f = x => Math.exp(-1 / x ** 3) - .4
chandrupatla(f)
Math.pow(-Math.log(.4), -1 / 3)
Ê = 0.61489
import { Sidi_method } from "0d9882008faed17b"
tests = [
{ f: x => Math.pow(x, 3) - 2 * x - 5, a: -2.8, b: 3.8 },
{ f: x => 1 - Math.pow(x, -2) },
{ f: x => Math.pow(x - 3, 3) },
{ f: x => 6 * Math.pow(x - 2, 5) },
{ f: x => Math.pow(x, 9), a: -1, b: 2 },
{ f: x => Math.pow(x, 19), a: -1, b: 2 },
{
f: x => (Math.abs(x) < 3.8e-4 ? 0 : x * Math.exp(-Math.pow(x, -2))),
a: -1,
b: 2
},
{
f: x =>
(-3062 * (1 - Ê) * Math.exp(-x)) / (Ê + (1 - Ê) * Math.exp(-x)) -
1013 +
1628 / x,
a: 1,
b: 2
},
{
f: x => Math.exp(x) - 2 - 0.01 / (x * x) + 0.000002 / Math.pow(x, 3),
a: 0.1,
b: 2
}
]
function draw(test) {
const f = test.f,
a = test.a || 0,
b = test.b || 0,
VV = chandrupatla(test.f, a, b),
points = VV.points,
xx = d3.range(...d3.extent(points, (d) => d[0]), 1e-2),
fx = xx.map(f);
const height = 500,
x = d3
.scaleLinear()
.domain(d3.extent(points, (d) => d[0]))
.range([50, width - 30])
.nice(),
y = d3
.scaleLinear()
.domain(d3.extent(fx.filter((d) => Math.abs(d) < 100)))
.range([height - 10, 10])
.nice(),
svg = d3.create("svg").attr("viewBox", [0, 0, width, height]);
svg
.append("g")
.call(d3.axisBottom(x))
.attr("transform", `translate(0,${y(0)})`)
.call((g) => g.select("text").remove());
svg
.append("g")
.call(d3.axisLeft(y))
.attr("transform", `translate(${x(a)},0)`);
svg
.append("path")
.style("fill", "none")
.style("stroke", "#000")
.attr("d", d3.line()(d3.zip(xx.map(x), fx.map(y))));
svg
.append("path")
.style("fill", "none")
.style("stroke", "red")
.attr("d", d3.line()(points.map((d) => [x(d[0]), y(d[1])])));
svg
.append("g")
.style("fill", "none")
.style("stroke", "red")
.selectAll("circle")
.data(points)
.join("circle")
.attr("r", (d, i) => (points.length - i) * 2)
.attr("cx", (d) => x(d[0]))
.attr("cy", (d) => y(d[1]))
.style("opacity", 0)
.transition()
.delay((d, i) => i * 300)
.style("opacity", 1);
return svg.node();
}
tests.map(d => chandrupatla(d.f, d.a, d.b))
tests.map(d => Sidi_method(d.f, d.a, d.b))
function chandrupatla(
f,
b = 0,
a = 1,
verbose = true,
eps_m = null,
eps_a = null
) {
// adapted from Chandrupatla's algorithm as described in Scherer
// https://books.google.com/books?id=cC-8BAAAQBAJ&pg=PA95
// ported to JavaScript by Philippe Rivière based on Jason Sachs’ Python implementation
// https://www.embeddedrelated.com/showarticle/855.php
// https://github.com/jason-s
// Initialization
const { abs, max, min } = Math,
log = verbose ? [] : null,
points = verbose ? [] : null;
let count = 0;
// fil: if the default values are the same, open the interval
if (a === b) b = a + 1;
let c = a,
fa = f(a),
fb = f(b),
fc = fa,
Nmax = 100;
count += 2;
points.push([b, fb], [a, fa]);
// fil: instead of throwing, try to find a bracket by exagerating the secant method
while (fa * fb > 0) {
b += 1.2 * ((fb - fa) / (b - a));
fb = f(b);
count++;
if (verbose) log.push(`BRA b = ${b}`), points.push([b, fb]);
if (Nmax-- === 0) throw "not bracketing";
}
let t = 0.5;
let do_iqi = false;
// jms: some guesses for default values of the eps_m and eps_a settings
// based on machine precision... not sure exactly what to do here
const eps = 1e-16;
if (!eps_m) eps_m = eps;
if (!eps_a) eps_a = 2 * eps;
let xm, fm;
do {
// use t to linearly interpolate between a and b,
// and evaluate this function as our newest estimate xt
const xt = a + t * (b - a),
ft = f(xt);
count++;
if (verbose)
log.push(
`${
do_iqi ? "IQI" : "BIS"
} t=${t} xt=${xt} ft=${ft} a=${a} b=${b} c=${c}`
),
points.push([xt, ft]);
// update our history of the last few points so that
// - a is the newest estimate (we're going to update it from xt)
// - c and b get the preceding two estimates
// - a and b maintain opposite signs for f(a) and f(b)
if (ft * fa >= 0) {
c = a;
fc = fa;
} else {
c = b;
b = a;
fc = fb;
fb = fa;
}
a = xt;
fa = ft;
// set xm so that f(xm) is the minimum magnitude of f(a) and f(b)
if (abs(fa) < abs(fb)) {
xm = a;
fm = fa;
} else {
xm = b;
fm = fb;
}
if (fm == 0) return { x: xm, log, fm, points, count };
// Figure out values xi and phi
// to determine which method we should use next
let tol = 2 * eps_m * abs(xm) + eps_a,
tlim = tol / abs(b - c);
if (tlim > 0.5) return { x: xm, log, fm, points, count };
let xi = (a - b) / (c - b),
phi = (fa - fb) / (fc - fb);
do_iqi = phi ** 2 < xi && (1 - phi) ** 2 < 1 - xi;
// inverse quadratic interpolation
if (do_iqi)
t =
((fa / (fb - fa)) * fc) / (fb - fc) +
(((((c - a) / (b - a)) * fa) / (fc - fa)) * fb) / (fc - fb);
// bisection
else t = 0.5;
// limit to the range (tlim, 1-tlim)
t = min(1 - tlim, max(tlim, t));
} while (Nmax-- > 0);
return { x: xm, log, fm, points, count };
}
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
