diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 0000000..123bfff --- /dev/null +++ b/.editorconfig @@ -0,0 +1,7 @@ +root = true + +[*] +end_of_line = crlf +insert_final_newline = true +indent_style = space +indent_size = 2 \ No newline at end of file diff --git a/jake.yml b/jake.yml index efb78f9..acc11b5 100644 --- a/jake.yml +++ b/jake.yml @@ -14,6 +14,7 @@ builds: packages: sylvester: - sylvester + - space - vector - matrix - line diff --git a/src/matrix.js b/src/matrix.js index 4fb9c2e..c7c57ed 100644 --- a/src/matrix.js +++ b/src/matrix.js @@ -436,7 +436,98 @@ Sylvester.Matrix.prototype = { this.elements.push([elements[i]]); } return this; + }, + + eigenvalues: function() { + if (!this.elements.length === 2) { return null; } + if (!this.isSquare()) { return null; } + var elems = this.elements, + aa = elems[0][0], + ab = elems[0][1], + ba = elems[1][0], + bb = elems[1][1]; + var firsteigenvalue = 1.0 / 2.0 * (aa + bb + + Math.sqrt((aa - bb) * (aa - bb) + 4.0 * ab * ba)), + secondeigenvalue = 1.0 / 2.0 * (aa + bb + - Math.sqrt((aa - bb) * (aa - bb) + 4.0 * ab * ba)), + sortedvalues = [firsteigenvalue, secondeigenvalue].sort(); + return Sylvester.Vector.create(sortedvalues); + }, + + // For 2x2 matrices only. Assumes real eigenvalues exist. + // Norms to absolute value 1. + eigenvectors: function() { + if (!this.elements.length === 2) { return null; } + if (!this.isSquare()) { return null; } + + // if real eigenvalues exist + var elems = this.elements, + EV = this.eigenvalues().elements, + EPSILON = Sylvester.precision, + vectors; + + function computeRealEigenvector(m, ev) { + var aa = m[0][0], ab = m[0][1], ba = m[1][0], bb = m[1][1], + x, y, norm; + if (Math.abs(aa - ev) < EPSILON && Math.abs(ab) < EPSILON) { + // return (l-d, c), where l is the eigenvalue + if (ba < 0.0) { + x = -ev + bb; + y = -ba; + } else { + x = ev - bb; + y = ba; + } + } else { + // return (b, l-a), where l is the eigenvalue + if (ev - aa < 0.0) { + x = -ab; + y = -ev + aa; + } else { + x = ab; + y = ev - aa; + } + } + norm = Math.sqrt(x * x + y * y); + return [x / norm, y / norm]; + } + + if (Math.abs(EV[0] - EV[1]) < EPSILON) { + // the same real eigenvalues + if (Math.abs(elems[0][0] - EV[0]) < EPSILON + && Math.abs(elems[0][1]) < EPSILON && Math.abs(elems[1][0]) < EPSILON + && Math.abs(elems[1][1] - EV[0]) < EPSILON) { + // two eigenvectors which are a basis of IR^2 (2-dim eigenspace) + vectors = [1.0, 0.0, 0.0, 1.0]; + } else { + vectors = [computeRealEigenvector(elems, EV[0])[0], + computeRealEigenvector(elems, EV[0])[1], NaN, NaN]; + } + // only one eigenvector + } else { + // two different real eigenvalues + var firstEigenvector = computeRealEigenvector(elems, EV[0]), + secondEigenvector = computeRealEigenvector(elems, EV[1]); + if (firstEigenvector[0] * secondEigenvector[1] + - firstEigenvector[1] * secondEigenvector[0] > 0) { + vectors = [firstEigenvector[0], firstEigenvector[1], + secondEigenvector[0], secondEigenvector[1]]; + } else { + vectors = [firstEigenvector[0], firstEigenvector[1], + -secondEigenvector[0], -secondEigenvector[1]]; + } + } + + if (vectors.length === 4) { + return Sylvester.Matrix.create([ + [vectors[0], vectors[1]], + [vectors[2], vectors[3]] + ]); + } else { + return Sylvester.Vector.create([vectors[0], vectors[1]]); + } } + }; Sylvester.Matrix.prototype.toUpperTriangular = Sylvester.Matrix.prototype.toRightTriangular; diff --git a/src/space.js b/src/space.js new file mode 100644 index 0000000..077e541 --- /dev/null +++ b/src/space.js @@ -0,0 +1,28 @@ +Sylvester.Space = function(){}; +Sylvester.Space.prototype.norm=function(v){ + return Math.sqrt(this.dot(v,v)); +} +Sylvester.Space.prototype.metric=function(v1,v2){ + return this.norm(v2.map(function(e,i){return e-v1[i];})); +} +Sylvester.Space.prototype.dot=function(v1,v2){} +Sylvester.Spaces={}; + +Sylvester.Spaces.Eucledian=function(){}; +Sylvester.Spaces.Eucledian.prototype=new Sylvester.Space(); +Sylvester.Spaces.Eucledian.prototype.dot=function(v1,v2) { + var product=0; + for(var i=0;i 1) { theta = 1; } return Math.acos(theta); }, + + oAngleFrom: function(vector) { + if (this.elements.length !== vector.elements.length)return null; + var a,b; + if (this.elements.length == 2){ + a=this.to3D(); + b=vector.to3D(); + }else{ + a=this; + b=vector; + } + var cross=a.cross(b); + var dot=a.dot(b); + var sign=1; + cross.each(function(e){sign*=(e>=0?1:-1);}); + return sign*Math.atan2(a.cross(b).modulus(), a.dot(b)); + }, isParallelTo: function(vector) { var angle = this.angleFrom(vector); @@ -103,25 +115,30 @@ Sylvester.Vector.prototype = { add: function(vector) { var V = vector.elements || vector; if (this.elements.length !== V.length) { return null; } - return this.map(function(x, i) { return x + V[i-1]; }); + return this.map(function(x, i) { return x + V[i]; }); }, subtract: function(vector) { var V = vector.elements || vector; if (this.elements.length !== V.length) { return null; } - return this.map(function(x, i) { return x - V[i-1]; }); + return this.map(function(x, i) { return x - V[i]; }); }, multiply: function(k) { return this.map(function(x) { return x*k; }); }, - dot: function(vector) { + dot: function(vector,space) { + space=space||Sylvester.Vector.Space; var V = vector.elements || vector; - var i, product = 0, n = this.elements.length; - if (n !== V.length) { return null; } - while (n--) { product += this.elements[n] * V[n]; } - return product; + return space.dot(this.elements,V); + }, + + wedge: function(vector) { + var B = vector.elements || vector; + if (this.elements.length !== 2 || B.length !== 2) { return null; } + var A = this.elements; + return (A[0] * B[1]) - (A[1] * B[0]); }, cross: function(vector) { @@ -167,16 +184,12 @@ Sylvester.Vector.prototype = { }); }, - distanceFrom: function(obj) { + distanceFrom: function(obj,space) { if (obj.anchor || (obj.start && obj.end)) { return obj.distanceFrom(this); } var V = obj.elements || obj; if (V.length !== this.elements.length) { return null; } - var sum = 0, part; - this.each(function(x, i) { - part = x - V[i-1]; - sum += part * part; - }); - return Math.sqrt(sum); + space=space||Sylvester.Vector.Space; + return space.metric(V,this.elements); }, liesOn: function(line) { @@ -188,6 +201,7 @@ Sylvester.Vector.prototype = { }, rotate: function(t, obj) { + obj=obj||Sylvester.Vector.Origin; var V, R = null, x, y, z; if (t.determinant) { R = t.elements; } switch (this.elements.length) { @@ -230,7 +244,7 @@ Sylvester.Vector.prototype = { // obj is a point var Q = obj.elements || obj; if (this.elements.length !== Q.length) { return null; } - return this.map(function(x, i) { return Q[i-1] + (Q[i-1] - x); }); + return this.map(function(x, i) { return Q[i] + (Q[i] - x); }); } }, @@ -260,3 +274,5 @@ Sylvester.Vector.prototype.each = Sylvester.Vector.prototype.forEach; Sylvester.Vector.i = Sylvester.Vector.create([1,0,0]); Sylvester.Vector.j = Sylvester.Vector.create([0,1,0]); Sylvester.Vector.k = Sylvester.Vector.create([0,0,1]); +Sylvester.Vector.Origin = Sylvester.Vector.Zero(); +Sylvester.Vector.Space=new Sylvester.Spaces.Eucledian(); \ No newline at end of file diff --git a/test/specs/matrix_spec.js b/test/specs/matrix_spec.js index d93036f..cae21cb 100644 --- a/test/specs/matrix_spec.js +++ b/test/specs/matrix_spec.js @@ -334,4 +334,27 @@ JS.ENV.MatrixSpec = JS.Test.describe("Matrix", function() { with(this) { [0,0,0,7] ])) }}) + + test("Eigenvalues", function() { with(this) { + assert($M([ + [1, -2], + [-2, 0] + ]).eigenvalues() + .eql([(1-Math.sqrt(17))/2, (1+Math.sqrt(17))/2])) + }}) + + test("Eigenvectors", function() { with(this) { + assert( + $M([ + [0.8, 0.3], + [0.2, 0.7] + ]) + .eigenvectors() + .eql([ + [-0.7071067811865475, 0.7071067811865476], + [-0.8320502943378437, -0.554700196225229] + ]) + ) + }}) + }})