30:T3a98,fsolve = { /** * @file Root finding for scalar (n=1) and vector (n>1) functions. * Always accepts an array initial guess. For length 1, takes a scalar fast path. */ /** * Find roots of a nonlinear system F(x) = 0 where x is an array. * - If `x0.length === 1`: scalar mode. Returns a number (root) unless `fullOutput` is true. * - If `x0.length > 1` : vector mode. Returns an array of roots unless `fullOutput` is true. * * Method selection: * - If `options.fprime` provided: * * n=1: Newton's method with derivative. * * n>1: Newton's method with Jacobian. * - Else: * * n=1: Secant method (derivative-free). * * n>1: Newton's method with numerically approximated Jacobian (forward differences). * * Convergence (residual-based): * residual <= max(atol, rtol * residual_initial) * - Scalar: residual = |f(x)|, residual_initial = |f(x0)| * - Vector: residual = ||F(x)||_inf, residual_initial = ||F(x0)||_inf * * @param {function(number[]): number[] } func * Function that maps x (array) -> residual array F(x) of same length. * For n=1, it should return a single-element array: [f(x)]. * @param {number[]} x0 * Initial guess as an array (length >= 1). * @param {Object} [options] * @param {function(number[]): (number[][]|number[]|number)} [options.fprime=null] * Jacobian provider. For n>1: J(x) must be n x n matrix (number[][]). * For n=1: may return [[df/dx]], [df/dx], or df/dx (number). * @param {number} [options.xtol=1e-8] * Backward-compatible absolute tolerance used when `atol` is not provided. * @param {number} [options.atol] * Absolute tolerance on the residual. Defaults to `xtol` if omitted. * @param {number} [options.rtol=0] * Relative tolerance on the residual, scaled by the initial residual. * @param {number} [options.maxfev=200] * Maximum number of function evaluations of `func`. (Numeric Jacobian & backtracking * consume additional evaluations.) * @param {boolean} [options.fullOutput=false] * If true, returns an info object: { root, iterations, funcEvals, converged }. * * @returns {number|number[]|{root:(number|number[]), iterations:number, funcEvals:number, converged:boolean}} * - n=1 & !fullOutput: number * - n>1 & !fullOutput: number[] * - fullOutput: diagnostics object with corresponding root type * * @throws {TypeError} If types/dimensions are invalid or `func(x0)` is malformed. * * @example * // --- Scalar (n=1) --- * // f(x) = x^2 - 2, with derivative. Always pass x0 as array. * const info1 = fsolve( * x => [x[0]*x[0] - 2], * [1], * { fprime: x => [[2*x[0]]], atol: 1e-12, fullOutput: true } * ); * console.log(info1.root); // ~1.41421356237 * * // Without derivative (secant): * const r1 = fsolve(x => [Math.cos(x[0]) - x[0]], [0.7]); * * // --- Vector (n=2) --- * // Solve: x + y = 3, x*y = 2 * const F = ([x, y]) => [x + y - 3, x*y - 2]; * const root2 = fsolve(F, [0.5, 2.5], { atol: 1e-12, rtol: 1e-10 }); * // -> near [1, 2] */ function fsolve(func, x0, options = {}) { if (typeof func !== "function") { throw new TypeError("fsolve: `func` must be a function."); } if (!Array.isArray(x0) || x0.length < 1) { throw new TypeError("fsolve: `x0` must be a non-empty array of numbers."); } if (!x0.every(Number.isFinite)) { throw new TypeError( "fsolve: all entries in `x0` must be finite numbers." ); } const { fprime = null, xtol = 1e-8, atol = undefined, rtol = 0, maxfev = 200, fullOutput = false } = options; const absTol = atol === undefined ? xtol : atol; const relTol = rtol === undefined ? 0 : rtol; const n = x0.length; if (n === 1) { // -------- Scalar fast path (but API is array-in) -------- const x0s = x0[0]; // Wrap array-form func into scalar form const f1 = (xs) => { const out = func([xs]); if ( !Array.isArray(out) || out.length !== 1 || !Number.isFinite(out[0]) ) { throw new TypeError( "fsolve: scalar mode expects `func([x])` -> [finite number]." ); } return out[0]; }; // Optional derivative wrapper (accepts number | [number] | [[number]]) let df1 = null; if (typeof fprime === "function") { df1 = (xs) => { const val = fprime([xs]); return scalarize(val); }; } const f0 = f1(x0s); const resRef = Math.abs(f0); const result = typeof df1 === "function" ? newtonScalar(f1, df1, x0s, absTol, relTol, resRef, maxfev) : secantScalar(f1, x0s, absTol, relTol, resRef, maxfev); return fullOutput ? result : result.root; } // -------- Vector/system path (n > 1) -------- const Fx0 = func(x0); if (!isNumericArray(Fx0) || Fx0.length !== n) { throw new TypeError( `fsolve: vector mode requires func(x0) to return an array of length ${n}.` ); } const resRef = normInf(Fx0); const result = typeof fprime === "function" ? newtonSystem(func, fprime, x0.slice(), absTol, relTol, resRef, maxfev) : newtonSystemNumericJac( func, x0.slice(), absTol, relTol, resRef, maxfev ); return fullOutput ? result : result.root; } /* ===================== Scalar solvers ===================== */ function newtonScalar(f, fp, x0, atol, rtol, resRef, maxfev) { let x = x0; let fev = 0; let iter = 0; let converged = false; const thr = convThreshold(atol, rtol, resRef); while (fev < maxfev) { const fx = f(x); fev++; if (!Number.isFinite(fx)) break; if (Math.abs(fx) <= thr) { converged = true; break; } const fpx = fp(x); if (!Number.isFinite(fpx) || fpx === 0) break; const xNext = x - fx / fpx; if (!Number.isFinite(xNext)) break; x = xNext; iter++; } return { root: x, iterations: iter, funcEvals: fev, converged }; } function secantScalar(f, x0, atol, rtol, resRef, maxfev) { const baseDelta = 1e-4; const delta = Math.abs(x0) > 1 ? Math.abs(x0) * baseDelta : baseDelta; let xPrev = x0; let x = x0 + delta; let fPrev = f(xPrev); let fx = f(x); let fev = 2; let iter = 0; const thr = convThreshold(atol, rtol, resRef); let converged = (Number.isFinite(fPrev) && Math.abs(fPrev) <= thr) || (Number.isFinite(fx) && Math.abs(fx) <= thr); while (!converged && fev < maxfev) { const denom = fx - fPrev; if (!Number.isFinite(denom) || Math.abs(denom) <= Number.EPSILON) { // Perturb slightly x = x + delta; fx = f(x); fev++; iter++; converged = Number.isFinite(fx) && Math.abs(fx) <= thr; if (!Number.isFinite(x) || fev >= maxfev) break; continue; } const xNext = x - (fx * (x - xPrev)) / denom; if (!Number.isFinite(xNext)) break; xPrev = x; fPrev = fx; x = xNext; fx = f(x); fev++; iter++; converged = Number.isFinite(fx) && Math.abs(fx) <= thr; } return { root: x, iterations: iter, funcEvals: fev, converged }; } /* ===================== Vector/system solvers ===================== */ function newtonSystem(F, J, x, atol, rtol, resRef, maxfev) { let fev = 0; let iter = 0; let converged = false; const thr = convThreshold(atol, rtol, resRef); while (fev < maxfev) { const Fx = F(x); fev++; if (!isNumericArray(Fx) || Fx.length !== x.length) break; const res = normInf(Fx); if (res <= thr) { converged = true; break; } const Jx = J(x); if (!isNumericMatrix(Jx, x.length)) break; // Solve Jx * s = Fx, update xNext = x - s const s = solveLinearSystem(Jx, Fx); if ( !isNumericArray(s) || s.length !== x.length || !s.every(Number.isFinite) ) break; let xNext = sub(x, s); // Backtracking if residual increases const FxNorm = res; let FxNext = F(xNext); fev++; let FxNextNorm = normInf(FxNext); let alpha = 1.0; const maxBacktracks = 8; let bt = 0; while (FxNextNorm > FxNorm && bt < maxBacktracks && fev < maxfev) { alpha *= 0.5; xNext = sub(x, scale(s, alpha)); FxNext = F(xNext); fev++; FxNextNorm = normInf(FxNext); bt++; } x = xNext; iter++; if (FxNextNorm <= thr) { converged = true; break; } } return { root: x, iterations: iter, funcEvals: fev, converged }; } function newtonSystemNumericJac(F, x, atol, rtol, resRef, maxfev) { const n = x.length; let fev = 0; let iter = 0; let converged = false; // Evaluate F(x) once let Fx = F(x); fev++; if (!isNumericArray(Fx) || Fx.length !== n) { throw new TypeError( "fsolve: numeric-Jacobian mode requires F(x) to return an array of length n." ); } const thr = convThreshold(atol, rtol, resRef); while (fev < maxfev) { const res = normInf(Fx); if (res <= thr) { converged = true; break; } // Forward-difference Jacobian const Jx = new Array(n); for (let i = 0; i < n; i++) Jx[i] = new Array(n).fill(0); for (let j = 0; j < n; j++) { const h = stepForFD(x[j]); const xPert = x.slice(); xPert[j] = x[j] + h; const Fpert = F(xPert); fev++; if (!isNumericArray(Fpert) || Fpert.length !== n) { throw new TypeError( "fsolve: F(x+h e_j) must return an array of length n." ); } for (let i = 0; i < n; i++) { Jx[i][j] = (Fpert[i] - Fx[i]) / h; } } const s = solveLinearSystem(Jx, Fx); // Jx * s = Fx if (!isNumericArray(s) || s.length !== n || !s.every(Number.isFinite)) break; // Trial step: xNext = x - s let xNext = sub(x, s); // Backtracking if residual increases const FxNorm = normInf(Fx); let FxNext = F(xNext); fev++; let FxNextNorm = normInf(FxNext); let alpha = 1.0; const maxBacktracks = 8; let bt = 0; while (FxNextNorm > FxNorm && bt < maxBacktracks && fev < maxfev) { alpha *= 0.5; xNext = sub(x, scale(s, alpha)); FxNext = F(xNext); fev++; FxNextNorm = normInf(FxNext); bt++; } x = xNext; Fx = FxNext; iter++; } return { root: x, iterations: iter, funcEvals: fev, converged }; } /* ===================== Linear algebra helpers ===================== */ /** * Solve A x = b using Gaussian elimination with partial pivoting. * @param {number[][]} A - n x n matrix * @param {number[]} b - length n vector * @returns {number[]} x */ function solveLinearSystem(A, b) { const n = A.length; // Deep copy const M = A.map((row) => row.slice()); const y = b.slice(); // Forward elimination with partial pivoting for (let k = 0; k < n; k++) { // Find pivot let pivotRow = k; let pivotVal = Math.abs(M[k][k]); for (let i = k + 1; i < n; i++) { const val = Math.abs(M[i][k]); if (val > pivotVal) { pivotVal = val; pivotRow = i; } } if (pivotVal === 0 || !Number.isFinite(pivotVal)) { return new Array(n).fill(NaN); // singular } if (pivotRow !== k) { const tmpRow = M[k]; M[k] = M[pivotRow]; M[pivotRow] = tmpRow; const tmpY = y[k]; y[k] = y[pivotRow]; y[pivotRow] = tmpY; } // Normalize pivot row const pivot = M[k][k]; for (let j = k; j < n; j++) M[k][j] /= pivot; y[k] /= pivot; // Eliminate below for (let i = k + 1; i < n; i++) { const factor = M[i][k]; for (let j = k; j < n; j++) M[i][j] -= factor * M[k][j]; y[i] -= factor * y[k]; } } // Back substitution const x = new Array(n); for (let i = n - 1; i >= 0; i--) { let sum = y[i]; for (let j = i + 1; j < n; j++) sum -= M[i][j] * x[j]; x[i] = sum; // M[i][i] is 1 due to normalization } return x; } /* ===================== Utilities ===================== */ function isNumericArray(arr) { return ( Array.isArray(arr) && arr.every((v) => typeof v === "number" && Number.isFinite(v)) ); } function isNumericMatrix(A, nExpected) { return ( Array.isArray(A) && A.length === nExpected && A.every( (row) => Array.isArray(row) && row.length === nExpected && row.every((v) => typeof v === "number" && Number.isFinite(v)) ) ); } function normInf(v) { let m = 0; for (let i = 0; i < v.length; i++) { const a = Math.abs(v[i]); if (a > m) m = a; } return m; } function sub(a, b) { const n = a.length; const out = new Array(n); for (let i = 0; i < n; i++) out[i] = a[i] - b[i]; return out; } function scale(a, alpha) { const n = a.length; const out = new Array(n); for (let i = 0; i < n; i++) out[i] = alpha * a[i]; return out; } /** * Finite difference step size heuristic: * h = sqrt(eps) * max(1, |x|) */ function stepForFD(xj) { const eps = Number.EPSILON; // ~2.22e-16 const base = Math.sqrt(eps); const mag = Math.max(1, Math.abs(xj)); return base * mag; } /** Residual convergence threshold: max(atol, rtol * resRef) */ function convThreshold(atol, rtol, resRef) { const absTol = Math.max(0, atol); const relTol = Math.max(0, rtol); return Math.max(absTol, relTol * resRef); } /** Convert number | [number] | [[number]] to a scalar number (for n=1 derivative). */ function scalarize(val) { if (typeof val === "number") return val; if (Array.isArray(val)) { if (val.length === 0) throw new TypeError( "fsolve: empty derivative/Jacobian in scalar mode." ); if (Array.isArray(val[0])) { // [[df/dx]] if (val.length !== 1 || val[0].length !== 1) { throw new TypeError("fsolve: 1x1 Jacobian expected in scalar mode."); } return val[0][0]; } // [df/dx] if (val.length !== 1) throw new TypeError( "fsolve: single derivative entry expected in scalar mode." ); return val[0]; } throw new TypeError( "fsolve: unsupported derivative/Jacobian type in scalar mode." ); } return fsolve; }31:[] 2e:{"title":"Newton's method (root finding)","cells":[{"id":0,"value":"# Newton's method (root finding)\nFind the root of a function (the input for which the function equals to 0)\n\n${tex.block`{\\displaystyle f'(x_{n})={\\dfrac {f(x_{n})-0}{x_{n}-x_{n+1}}}}`}\n${tex.block`{\\displaystyle x_{n+1}=x_{n}-{\\frac {f(x_{n})}{f'(x_{n})}}}`}","pinned":false,"hidden":false,"mode":"md"},{"id":190,"value":"# TODO\n\n- [ ] naming convention\n- [ ] solve an array of numbers","pinned":false,"hidden":false,"mode":"md"},{"id":160,"value":"## References\n- [Wikipedia: Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)\n- [Wikipedia: Secant Method](https://en.wikipedia.org/wiki/Secant_method)\n- [Scipy: fsolve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html)","pinned":false,"hidden":false,"mode":"md"},{"id":126,"value":"function f(x) {\n return Math.exp(x) - 2;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":130,"value":"viewof guess = Inputs.range([0, 10], { label: \"Initial Guess\", step: 0.1 })","pinned":true,"hidden":false,"mode":"ojs"},{"id":133,"value":"viewof tol = Inputs.range([0.0001, 0.1], { label: \"Tolerance\", step: 1e-4 })","pinned":true,"hidden":false,"mode":"ojs"},{"id":19,"value":"root = fsolve((x) => [f(x[0])], [guess], { tol })","pinned":true,"hidden":false,"mode":"ojs"},{"id":115,"value":"f(root)","pinned":true,"hidden":false,"mode":"ojs"},{"id":155,"value":"{\n let f = (x) => [x * x[0] - 2];\n let fPrime = (x) => [2 * x[0]];\n return fsolve(f, [1], { fPrime, tol: 0.0001 });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":164,"value":"$30","pinned":true,"hidden":false,"mode":"ojs"},{"id":167,"value":"{\n let f = (x) => [x[0] * x[0] - 2];\n let fPrime = (x) => [2 * x[0]];\n\n return fsolve(f, [1], { fPrime });\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W31","stdlib":"1","titling":"auto"}