Instantly share code, notes, and snippets.

# Fil/newton.js

Last active August 29, 2015 13:56
Show Gist options
• Save Fil/8903368 to your computer and use it in GitHub Desktop.
newton function to solve f(V)=0
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 /* * Newton's method for finding roots * * code adapted from D.V. Fedorov, * “Introduction to Numerical Methods with examples in Javascript” * http://owww.phys.au.dk/~fedorov/nucltheo/Numeric/11/book.pdf * (licensed under the GPL) * by Philippe Riviere March 2014 * modified for compatibility with Chrome/Safari * added a max iterations parameter * * Usage: Newton.Solve(Math.exp, 2)); => ~ log(2) * Newton.Solve(d3.geo.chamberlin(), [200,240]) */ (function() { Newton = {version: "1.0.0"}; // semver Newton.Norm = function (v) { return Math.sqrt(v.reduce(function (s, e) { return s + e * e }, 0)); } Newton.Dot = function (a, b) { var s = 0; for (var i in a) s += a[i] * b[i]; return s; } // QR - decomposition A=QR of matrix A Newton.QRDec = function(A){ var m = A.length, R = []; for (var j in A) { R[j] = []; for (var i in A) R[j][i] = 0; } var Q = []; for (var i in A) { Q[i] = []; for (var j in A[0]) Q[i][j] = A[i][j]; } // Q is a copy of A for (var i = 0; i < m; i++) { var e = Q[i], r = Math.sqrt(Newton.Dot(e, e)); if (r == 0) throw "Newton.QRDec: singular matrix" R[i][i] = r; for (var k in e) e[k] /= r; // normalization for (var j = i + 1; j < m; j++) { var q = Q[j], s = Newton.Dot(e, q); for (var k in q) q[k] -= s * e[k]; // orthogonalization R[j][i] = s; } } return [Q, R]; }; // QR - backsubstitution // input: matrices Q,R, array b; output: array x such that QRx=b Newton.QRBack = function(Q, R, b) { var m = Q.length, c = new Array(m), x = new Array(m); for (var i in Q) { // c = QˆT b c[i] = 0; for (var k in b) c[i] += Q[i][k] * b[k]; } for (var i = m - 1; i >= 0; i--) { // back substitution for (var s = 0, k = i + 1; k < m; k++) s += R[k][i] * x[k]; x[i] = (c[i] - s) / R[i][i]; } return x; } // calculates inverse of matrix A Newton.Inverse = function(A){ var t = Newton.QRDec(A), Q = t[0], R = t[1]; var m = []; for (i in A) { n = []; for (k in A) { n[k] = (k == i ? 1 : 0) } m[i] = Newton.QRBack(Q, R, n); } return m; }; Newton.Zero = function (fs, x, opts={} /* acc, dx, max */) { // Newton's root-finding method if (opts.acc == undefined) opts.acc = 1e-6 if (opts.dx == undefined) opts.dx = 1e-3 if (opts.max == undefined) opts.max = 40 // max iterations var J = []; for (i in x) { J[i] = []; for (j in x) J[i][j] = 0; } var minusfx = []; var v = fs(x); if (v == null) throw "unable to compute fs at "+JSON.stringify(x) for (i in x) minusfx[i] = -v[i]; do { if (opts.max-- < 0) return null; for (i in x) for (k in x) { // calculate Jacobian x[k] += opts.dx v = fs(x); if (v == null) throw "unable to compute fs at "+JSON.stringify(x) J[k][i] = (v[i] + minusfx[i]) / opts.dx x[k] -= opts.dx } var t = Newton.QRDec(J), Q = t[0], R = t[1], Dx = Newton.QRBack(Q, R, minusfx) // Newton's step var s = 2 do { // simple backtracking line search s = s / 2; var z = []; for (i in x) { z[i] = x[i] + s * Dx[i]; } var minusfz = []; v = fs(z); if (v == null) throw "unable to compute fs at "+JSON.stringify(z) for (i in x) { minusfz[i] = -v[i]; } } while (Newton.Norm(minusfz) > (1 - s / 2) * Newton.Norm(minusfx) && s > 1. / 128) minusfx = minusfz; x = z; // step done } while (Newton.Norm(minusfx) > opts.acc) return x; } Newton.Solve = function(fs, res, opts={}){ if (typeof res != 'object') res = [ typeof res == 'number' ? + res : 0 ]; var _fs = fs, fs = function(x) { var r = _fs(x); if (typeof r == 'number') r = [ r ]; for (var i in r) r[i] -= res[i]; return r; } var start = []; if (opts.start) { start = opts.start; } else { for (var i in res) start[i]=0; } var n = Newton.Zero(fs, start, opts); if (n && n.length == 1) n = n[0]; return n; }; })();

### Fil commented Feb 9, 2014

il manque un paramètre : nb max d'itérations avant de déclarer forfait, en cas de non convergence

max=10 ajouté

### Fil commented Feb 18, 2014

TODO: works only on Firefox, throws syntax error on Chrome and Safari

### Fil commented Feb 19, 2014

compatiblilty with Webkit solved

todo: proper wrapping; f(V) should be a vector result of one function instead of a vetor of functions.