diff --git a/AUTHORS b/AUTHORS index 16fa2c7dd..ae1ceb58b 100644 --- a/AUTHORS +++ b/AUTHORS @@ -199,6 +199,7 @@ FelixSelter <55546882+FelixSelter@users.noreply.github.com> Hansuku <1556207795@qq.com> Windrill Carl Osterwisch +Lucas Egidio Isaac B <66892358+isaacbyr@users.noreply.github.com> Ajinkya Chikhale <86607732+ajinkyac03@users.noreply.github.com> Abdullah Adeel <64294045+mabdullahadeel@users.noreply.github.com> diff --git a/src/expression/embeddedDocs/embeddedDocs.js b/src/expression/embeddedDocs/embeddedDocs.js index 19779e38f..bd620959b 100644 --- a/src/expression/embeddedDocs/embeddedDocs.js +++ b/src/expression/embeddedDocs/embeddedDocs.js @@ -78,6 +78,9 @@ import { roundDocs } from './function/arithmetic/round.js' import { signDocs } from './function/arithmetic/sign.js' import { sqrtDocs } from './function/arithmetic/sqrt.js' import { sqrtmDocs } from './function/arithmetic/sqrtm.js' +import { sylvesterDocs } from './function/matrix/sylvester.js' +import { schurDocs } from './function/matrix/schur.js' +import { lyapDocs } from './function/matrix/lyap.js' import { squareDocs } from './function/arithmetic/square.js' import { subtractDocs } from './function/arithmetic/subtract.js' import { unaryMinusDocs } from './function/arithmetic/unaryMinus.js' @@ -467,6 +470,9 @@ export const embeddedDocs = { zeros: zerosDocs, fft: fftDocs, ifft: ifftDocs, + sylvester: sylvesterDocs, + schur: schurDocs, + lyap: lyapDocs, // functions - probability combinations: combinationsDocs, diff --git a/src/expression/embeddedDocs/function/matrix/lyap.js b/src/expression/embeddedDocs/function/matrix/lyap.js new file mode 100644 index 000000000..8c25d4298 --- /dev/null +++ b/src/expression/embeddedDocs/function/matrix/lyap.js @@ -0,0 +1,15 @@ +export const lyapDocs = { + name: 'lyap', + category: 'Matrix', + syntax: [ + 'lyap(A,Q)' + ], + description: 'Solves the Continuous-time Lyapunov equation AP+PA\'+Q=0 for P', + examples: [ + 'lyap([[-2, 0], [1, -4]], [[3, 1], [1, 3]])', + 'lyap(A,Q)' + ], + seealso: [ + 'schur', 'sylvester' + ] +} diff --git a/src/expression/embeddedDocs/function/matrix/schur.js b/src/expression/embeddedDocs/function/matrix/schur.js new file mode 100644 index 000000000..9256a4923 --- /dev/null +++ b/src/expression/embeddedDocs/function/matrix/schur.js @@ -0,0 +1,15 @@ +export const schurDocs = { + name: 'schur', + category: 'Matrix', + syntax: [ + 'schur(A)' + ], + description: 'Performs a real Schur decomposition of the real matrix A = UTU\'', + examples: [ + 'schur([[1, 0], [-4, 3]])', + 'schur(A)' + ], + seealso: [ + 'lyap', 'sylvester' + ] +} diff --git a/src/expression/embeddedDocs/function/matrix/sylvester.js b/src/expression/embeddedDocs/function/matrix/sylvester.js new file mode 100644 index 000000000..913ccaac9 --- /dev/null +++ b/src/expression/embeddedDocs/function/matrix/sylvester.js @@ -0,0 +1,15 @@ +export const sylvesterDocs = { + name: 'sylvester', + category: 'Matrix', + syntax: [ + 'sylvester(A,B,C)' + ], + description: 'Solves the real-valued Sylvester equation AX+XB=C for X', + examples: [ + 'sylvester([[-1, -2], [1, 1]], [[-2, 1], [-1, 2]], [[-3, 2], [3, 0]])', + 'sylvester(A,B,C)' + ], + seealso: [ + 'schur', 'lyap' + ] +} diff --git a/src/factoriesAny.js b/src/factoriesAny.js index 2ffc291c2..72491a06b 100644 --- a/src/factoriesAny.js +++ b/src/factoriesAny.js @@ -220,6 +220,9 @@ export { createPinv } from './function/matrix/pinv.js' export { createEigs } from './function/matrix/eigs.js' export { createExpm } from './function/matrix/expm.js' export { createSqrtm } from './function/matrix/sqrtm.js' +export { createSylvester } from './function/algebra/sylvester.js' +export { createSchur } from './function/algebra/decomposition/schur.js' +export { createLyap } from './function/algebra/lyap.js' export { createDivide } from './function/arithmetic/divide.js' export { createDistance } from './function/geometry/distance.js' export { createIntersect } from './function/geometry/intersect.js' diff --git a/src/function/algebra/decomposition/schur.js b/src/function/algebra/decomposition/schur.js new file mode 100644 index 000000000..1beebd567 --- /dev/null +++ b/src/function/algebra/decomposition/schur.js @@ -0,0 +1,77 @@ +import { factory } from '../../../utils/factory.js' + +const name = 'schur' +const dependencies = [ + 'typed', + 'matrix', + 'identity', + 'multiply', + 'qr', + 'norm', + 'subtract' +] + +export const createSchur = /* #__PURE__ */ factory(name, dependencies, ( + { + typed, + matrix, + identity, + multiply, + qr, + norm, + subtract + } +) => { + /** + * + * Performs a real Schur decomposition of the real matrix A = UTU' where U is orthogonal + * and T is upper quasi-triangular. + * https://en.wikipedia.org/wiki/Schur_decomposition + * + * Syntax: + * + * math.schur(A) + * + * Examples: + * + * const A = [[1, 0], [-4, 3]] + * math.schur(A) // returns {T: [[3, 4], [0, 1]], R: [[0, 1], [-1, 0]]} + * + * See also: + * + * sylvester, lyap, qr + * + * @param {Array | Matrix} A Matrix A + * @return {{U: Array | Matrix, T: Array | Matrix}} Object containing both matrix U and T of the Schur Decomposition A=UTU' + */ + return typed(name, { + Array: function (X) { + const r = _schur(matrix(X)) + return { + U: r.U.valueOf(), + T: r.T.valueOf() + } + }, + + Matrix: function (X) { + return _schur(X) + } + }) + function _schur (X) { + const n = X.size()[0] + let A = X + let U = identity(n) + let k = 0 + let A0 + do { + A0 = A + const QR = qr(A) + const Q = QR.Q + const R = QR.R + A = multiply(R, Q) + U = multiply(U, Q) + if ((k++) > 100) { break } + } while (norm(subtract(A, A0)) > 1e-4) + return { U, T: A } + } +}) diff --git a/src/function/algebra/lyap.js b/src/function/algebra/lyap.js new file mode 100644 index 000000000..34cb0bf6b --- /dev/null +++ b/src/function/algebra/lyap.js @@ -0,0 +1,61 @@ +import { factory } from '../../utils/factory.js' + +const name = 'lyap' +const dependencies = [ + 'typed', + 'matrix', + 'sylvester', + 'multiply', + 'transpose' +] + +export const createLyap = /* #__PURE__ */ factory(name, dependencies, ( + { + typed, + matrix, + sylvester, + multiply, + transpose + } +) => { + /** + * + * Solves the Continuous-time Lyapunov equation AP+PA'+Q=0 for P, where + * Q is an input matrix. When Q is symmetric, P is also symmetric. Notice + * that different equivalent definitions exist for the Continuous-time + * Lyapunov equation. + * https://en.wikipedia.org/wiki/Lyapunov_equation + * + * Syntax: + * + * math.lyap(A, Q) + * + * Examples: + * + * const A = [[-2, 0], [1, -4]] + * const Q = [[3, 1], [1, 3]] + * const P = math.lyap(A, Q) + * + * See also: + * + * sylvester, schur + * + * @param {Matrix | Array} A Matrix A + * @param {Matrix | Array} Q Matrix Q + * @return {Matrix | Array} Matrix P solution to the Continuous-time Lyapunov equation AP+PA'=Q + */ + return typed(name, { + 'Matrix, Matrix': function (A, Q) { + return sylvester(A, transpose(A), multiply(-1, Q)) + }, + 'Array, Matrix': function (A, Q) { + return sylvester(matrix(A), transpose(matrix(A)), multiply(-1, Q)) + }, + 'Matrix, Array': function (A, Q) { + return sylvester(A, transpose(matrix(A)), matrix(multiply(-1, Q))) + }, + 'Array, Array': function (A, Q) { + return sylvester(matrix(A), transpose(matrix(A)), matrix(multiply(-1, Q))).toArray() + } + }) +}) diff --git a/src/function/algebra/sylvester.js b/src/function/algebra/sylvester.js new file mode 100644 index 000000000..ab882374a --- /dev/null +++ b/src/function/algebra/sylvester.js @@ -0,0 +1,145 @@ +import { factory } from '../../utils/factory.js' + +const name = 'sylvester' +const dependencies = [ + 'typed', + 'schur', + 'matrixFromColumns', + 'matrix', + 'multiply', + 'range', + 'concat', + 'transpose', + 'index', + 'subset', + 'add', + 'subtract', + 'identity', + 'lusolve', + 'abs' +] + +export const createSylvester = /* #__PURE__ */ factory(name, dependencies, ( + { + typed, + schur, + matrixFromColumns, + matrix, + multiply, + range, + concat, + transpose, + index, + subset, + add, + subtract, + identity, + lusolve, + abs + } +) => { + /** + * + * Solves the real-valued Sylvester equation AX+XB=C for X, where A, B and C are + * matrices of appropriate dimensions, being A and B squared. Notice that other + * equivalent definitions for the Sylvester equation exist and this function + * assumes the one presented in the original publication of the the Bartels- + * Stewart algorithm, which is implemented by this function. + * https://en.wikipedia.org/wiki/Sylvester_equation + * + * Syntax: + * + * math.sylvester(A, B, C) + * + * Examples: + * + * const A = [[-1, -2], [1, 1]] + * const B = [[2, -1], [1, -2]] + * const C = [[-3, 2], [3, 0]] + * math.sylvester(A, B, C) // returns DenseMatrix [[-0.25, 0.25], [1.5, -1.25]] + * + * See also: + * + * schur, lyap + * + * @param {Matrix | Array} A Matrix A + * @param {Matrix | Array} B Matrix B + * @param {Matrix | Array} C Matrix C + * @return {Matrix | Array} Matrix X, solving the Sylvester equation + */ + return typed(name, { + 'Matrix, Matrix, Matrix': _sylvester, + 'Array, Matrix, Matrix': function (A, B, C) { + return _sylvester(matrix(A), B, C) + }, + 'Array, Array, Matrix': function (A, B, C) { + return _sylvester(matrix(A), matrix(B), C) + }, + 'Array, Matrix, Array': function (A, B, C) { + return _sylvester(matrix(A), B, matrix(C)) + }, + 'Matrix, Array, Matrix': function (A, B, C) { + return _sylvester(A, matrix(B), C) + }, + 'Matrix, Array, Array': function (A, B, C) { + return _sylvester(A, matrix(B), matrix(C)) + }, + 'Matrix, Matrix, Array': function (A, B, C) { + return _sylvester(A, B, matrix(C)) + }, + 'Array, Array, Array': function (A, B, C) { + return _sylvester(matrix(A), matrix(B), matrix(C)).toArray() + } + }) + function _sylvester (A, B, C) { + const n = B.size()[0] + const m = A.size()[0] + + const sA = schur(A) + const F = sA.T + const U = sA.U + const sB = schur(multiply(-1, B)) + const G = sB.T + const V = sB.U + const D = multiply(multiply(transpose(U), C), V) + const all = range(0, m) + const y = [] + + const hc = (a, b) => concat(a, b, 1) + const vc = (a, b) => concat(a, b, 0) + + for (let k = 0; k < n; k++) { + if (k < (n - 1) && abs(subset(G, index(k + 1, k))) > 1e-5) { + let RHS = vc(subset(D, index(all, k)), subset(D, index(all, k + 1))) + for (let j = 0; j < k; j++) { + RHS = add(RHS, + vc(multiply(y[j], subset(G, index(j, k))), multiply(y[j], subset(G, index(j, k + 1)))) + ) + } + const gkk = multiply(identity(m), multiply(-1, subset(G, index(k, k)))) + const gmk = multiply(identity(m), multiply(-1, subset(G, index(k + 1, k)))) + const gkm = multiply(identity(m), multiply(-1, subset(G, index(k, k + 1)))) + const gmm = multiply(identity(m), multiply(-1, subset(G, index(k + 1, k + 1)))) + const LHS = vc( + hc(add(F, gkk), gmk), + hc(gkm, add(F, gmm)) + ) + const yAux = lusolve(LHS, RHS) + y[k] = yAux.subset(index(range(0, m), 0)) + y[k + 1] = yAux.subset(index(range(m, 2 * m), 0)) + k++ + } else { + let RHS = subset(D, index(all, k)) + for (let j = 0; j < k; j++) { RHS = add(RHS, multiply(y[j], subset(G, index(j, k)))) } + const gkk = subset(G, index(k, k)) + const LHS = subtract(F, multiply(gkk, identity(m))) + + y[k] = lusolve(LHS, RHS) + } + } + const Y = matrix(matrixFromColumns(...y)) + const X = multiply(U, multiply(Y, transpose(V))) + + return X + } +}) diff --git a/test/node-tests/doc.test.js b/test/node-tests/doc.test.js index 0174e389a..7c0533a2e 100644 --- a/test/node-tests/doc.test.js +++ b/test/node-tests/doc.test.js @@ -92,7 +92,7 @@ const knownProblems = new Set([ 'mod', 'invmod', 'floor', 'fix', 'expm1', 'exp', 'dotPow', 'dotMultiply', 'dotDivide', 'divide', 'ceil', 'cbrt', 'add', 'usolveAll', 'usolve', 'slu', 'rationalize', 'qr', 'lusolve', 'lup', 'lsolveAll', 'lsolve', 'derivative', - 'symbolicEqual', 'map' + 'symbolicEqual', 'map', 'schur', 'sylvester' ]) function maybeCheckExpectation (name, expected, expectedFrom, got, gotFrom) { diff --git a/test/unit-tests/function/algebra/decomposition/schur.test.js b/test/unit-tests/function/algebra/decomposition/schur.test.js new file mode 100644 index 000000000..790a5fb82 --- /dev/null +++ b/test/unit-tests/function/algebra/decomposition/schur.test.js @@ -0,0 +1,63 @@ +// test schur decomposition +import assert from 'assert' + +import math from '../../../../../src/defaultInstance.js' + +describe('schur', function () { + it('should calculate schur decomposition of order 5 Array with numbers', function () { + assert.ok(math.norm(math.subtract(math.schur([ + [-5.3, -1.4, -0.2, 0.7, 1.0], + [-0.4, -1.0, -0.1, -1.2, 0.7], + [0.3, 0.7, -2.5, 0.7, -0.3], + [3.6, -0.1, 1.4, -2.4, 0.3], + [2.8, 0.7, 1.4, 0.5, -4.8] + ]).T, [ + [-6.97747746169558, 0.5046853036888738, -0.5269551218982134, 2.9902479419087253, -2.2914719859941908], + [7.686296479877504e-28, -3.667202530573731, -1.3776362163231233, 0.4680921120934126, -0.374760141366345], + [-4.92421882171775e-28, 0.0859443307361262, -3.798852922420336, 0.6595326982269121, 0.38704916773017245], + [-1.8971150504442668e-71, -1.3481560449785515e-44, -8.960727665382592e-44, -1.3867791434503591, 0.5989746088924175], + [5.3038070596554604e-164, -7.015716167009369e-138, 1.0415178925287816e-136, 3.3306222609902884e-93, -0.16968794185997693] + ])) < 1e-3) + assert.ok(math.norm(math.subtract(math.schur([ + [-5.3, -1.4, -0.2, 0.7, 1.0], + [-0.4, -1.0, -0.1, -1.2, 0.7], + [0.3, 0.7, -2.5, 0.7, -0.3], + [3.6, -0.1, 1.4, -2.4, 0.3], + [2.8, 0.7, 1.4, 0.5, -4.8] + ]).U, [ + [0.6039524392527362, -0.11955248228665324, 0.5309978859071411, -0.3239619623824945, 0.48377530651243317], + [0.034196874004165004, -0.3407725193032822, -0.05878741847605931, -0.7490924342670748, -0.5640117270489927], + [-0.024973926752867345, 0.6355774530247909, -0.45835474572889245, -0.5151114453511857, 0.3463938944685342], + [-0.42238909820980175, -0.6350073886392834, -0.26854304887535935, -0.13970317163522103, 0.5715948922294664], + [-0.6745633977804205, 0.24977632323107185, 0.6575567219332444, -0.22147775183161147, 0.03380493473031572] + ])) < 1e-3) + }) + it('should calculate schur decomposition of order 5 Matrix with numbers', function () { + assert.ok(math.norm(math.subtract(math.schur(math.matrix([ + [-5.3, -1.4, -0.2, 0.7, 1.0], + [-0.4, -1.0, -0.1, -1.2, 0.7], + [0.3, 0.7, -2.5, 0.7, -0.3], + [3.6, -0.1, 1.4, -2.4, 0.3], + [2.8, 0.7, 1.4, 0.5, -4.8] + ])).T, math.matrix([ + [-6.97747746169558, 0.5046853036888738, -0.5269551218982134, 2.9902479419087253, -2.2914719859941908], + [7.686296479877504e-28, -3.667202530573731, -1.3776362163231233, 0.4680921120934126, -0.374760141366345], + [-4.92421882171775e-28, 0.0859443307361262, -3.798852922420336, 0.6595326982269121, 0.38704916773017245], + [-1.8971150504442668e-71, -1.3481560449785515e-44, -8.960727665382592e-44, -1.3867791434503591, 0.5989746088924175], + [5.3038070596554604e-164, -7.015716167009369e-138, 1.0415178925287816e-136, 3.3306222609902884e-93, -0.16968794185997693] + ]))) < 1e-3) + assert.ok(math.norm(math.subtract(math.schur(math.matrix([ + [-5.3, -1.4, -0.2, 0.7, 1.0], + [-0.4, -1.0, -0.1, -1.2, 0.7], + [0.3, 0.7, -2.5, 0.7, -0.3], + [3.6, -0.1, 1.4, -2.4, 0.3], + [2.8, 0.7, 1.4, 0.5, -4.8] + ])).U, math.matrix([ + [0.6039524392527362, -0.11955248228665324, 0.5309978859071411, -0.3239619623824945, 0.48377530651243317], + [0.034196874004165004, -0.3407725193032822, -0.05878741847605931, -0.7490924342670748, -0.5640117270489927], + [-0.024973926752867345, 0.6355774530247909, -0.45835474572889245, -0.5151114453511857, 0.3463938944685342], + [-0.42238909820980175, -0.6350073886392834, -0.26854304887535935, -0.13970317163522103, 0.5715948922294664], + [-0.6745633977804205, 0.24977632323107185, 0.6575567219332444, -0.22147775183161147, 0.03380493473031572] + ]))) < 1e-3) + }) +}) diff --git a/test/unit-tests/function/algebra/lyap.test.js b/test/unit-tests/function/algebra/lyap.test.js new file mode 100644 index 000000000..4ac64937e --- /dev/null +++ b/test/unit-tests/function/algebra/lyap.test.js @@ -0,0 +1,49 @@ +// test lyapunov equation solver +import assert from 'assert' + +import math from '../../../../src/defaultInstance.js' + +describe('lyap', function () { + it('should solve lyapunov equation of order 5 with Matrices', function () { + assert.ok(math.norm(math.subtract(math.lyap(math.matrix([ + [-5.3, -1.4, -0.2, 0.7, 1.0], + [-0.4, -1.0, -0.1, -1.2, 0.7], + [0.3, 0.7, -2.5, 0.7, -0.3], + [3.6, -0.1, 1.4, -2.4, 0.3], + [2.8, 0.7, 1.4, 0.5, -4.8] + ]), math.matrix([ + [-2, 0, 0, 0, 0], + [0, -2, 0, 0, 0], + [0, 0, -2, 0, 0], + [0, 0, 0, -2, 0], + [0, 0, 0, 0, -2] + ])), math.matrix([ + [-0.8567032819332775, 1.6675327139319873, 0.05193348368861031, -1.2336381488442063, -0.3320481938129234], + [1.6675327139319864, -4.577292046354884, -0.24491693525451969, 2.7566635155267862, 0.533179377606574], + [0.05193348368861139, -0.24491693525452285, -0.4786127989547805, -0.10398529391846956, -0.10706505975853312], + [-1.233638148844205, 2.7566635155267827, -0.10398529391847175, -2.5199807906596963, -0.6187026623519729], + [-0.33204819381292294, 0.5331793776065724, -0.10706505975853345, -0.6187026623519725, -0.4199482902478158] + ])), 'fro') < 1e-3) + }) + it('should solve lyapunov equation of order 5 with Arrays', function () { + assert.ok(math.norm(math.subtract(math.lyap([ + [-5.3, -1.4, -0.2, 0.7, 1.0], + [-0.4, -1.0, -0.1, -1.2, 0.7], + [0.3, 0.7, -2.5, 0.7, -0.3], + [3.6, -0.1, 1.4, -2.4, 0.3], + [2.8, 0.7, 1.4, 0.5, -4.8] + ], [ + [-2, 0, 0, 0, 0], + [0, -2, 0, 0, 0], + [0, 0, -2, 0, 0], + [0, 0, 0, -2, 0], + [0, 0, 0, 0, -2] + ]), [ + [-0.8567032819332775, 1.6675327139319873, 0.05193348368861031, -1.2336381488442063, -0.3320481938129234], + [1.6675327139319864, -4.577292046354884, -0.24491693525451969, 2.7566635155267862, 0.533179377606574], + [0.05193348368861139, -0.24491693525452285, -0.4786127989547805, -0.10398529391846956, -0.10706505975853312], + [-1.233638148844205, 2.7566635155267827, -0.10398529391847175, -2.5199807906596963, -0.6187026623519729], + [-0.33204819381292294, 0.5331793776065724, -0.10706505975853345, -0.6187026623519725, -0.4199482902478158] + ]), 'fro') < 1e-3) + }) +}) diff --git a/test/unit-tests/function/algebra/sylvester.test.js b/test/unit-tests/function/algebra/sylvester.test.js new file mode 100644 index 000000000..aac909f63 --- /dev/null +++ b/test/unit-tests/function/algebra/sylvester.test.js @@ -0,0 +1,55 @@ +// test schur decomposition +import assert from 'assert' + +import math from '../../../../src/defaultInstance.js' + +describe('sylvester', function () { + it('should solve sylvester equation of order 5 with Arrays', function () { + assert.ok(math.norm(math.subtract(math.sylvester([ + [-5.3, -1.4, -0.2, 0.7], + [-0.4, -1.0, -0.1, -1.2], + [0.3, 0.7, -2.5, 0.7], + [3.6, -0.1, 1.4, -2.4] + ], [ + [1.1, -0.3, -0.9, 0.8, -2.5], + [-0.6, 2.6, 0.2, 0.4, 1.3], + [0.4, -0.5, -0.2, 0.2, -0.1], + [-0.4, -1.9, -0.2, 0.5, 1.4], + [0.4, -1.0, -0.1, -0.8, -1.3] + ], [ + [1.4, 1.1, -1.9, 0.1, 1.2], + [-1.7, 0.1, -0.4, 2.1, 0.5], + [1.9, 2.3, -0.8, -0.7, 1.0], + [-1.1, 2.8, -0.8, -0.3, 0.9] + ]), [ + [-0.19393862606643053, -0.17101629636521865, 2.709348263225366, -0.0963000767188319, 0.5244718194343121], + [0.38421326955977486, -0.21159588555260944, -6.544262021555474, -0.15113424769761136, -2.312533293658291], + [-2.2708235174374747, 4.498279916441834, 1.4553799673144823, -0.9300926971755248, 2.5508111398452353], + [-1.0935231387905004, 4.113817086842746, 5.747671819196675, -0.9408309030864932, 2.967655969930743] + ]), 'fro') < 1e-3) + }) + it('should solve sylvester equation of order 5 with Matrices', function () { + assert.ok(math.norm(math.subtract(math.sylvester(math.matrix([ + [-5.3, -1.4, -0.2, 0.7], + [-0.4, -1.0, -0.1, -1.2], + [0.3, 0.7, -2.5, 0.7], + [3.6, -0.1, 1.4, -2.4] + ]), math.matrix([ + [1.1, -0.3, -0.9, 0.8, -2.5], + [-0.6, 2.6, 0.2, 0.4, 1.3], + [0.4, -0.5, -0.2, 0.2, -0.1], + [-0.4, -1.9, -0.2, 0.5, 1.4], + [0.4, -1.0, -0.1, -0.8, -1.3] + ]), math.matrix([ + [1.4, 1.1, -1.9, 0.1, 1.2], + [-1.7, 0.1, -0.4, 2.1, 0.5], + [1.9, 2.3, -0.8, -0.7, 1.0], + [-1.1, 2.8, -0.8, -0.3, 0.9] + ])), math.matrix([ + [-0.19393862606643053, -0.17101629636521865, 2.709348263225366, -0.0963000767188319, 0.5244718194343121], + [0.38421326955977486, -0.21159588555260944, -6.544262021555474, -0.15113424769761136, -2.312533293658291], + [-2.2708235174374747, 4.498279916441834, 1.4553799673144823, -0.9300926971755248, 2.5508111398452353], + [-1.0935231387905004, 4.113817086842746, 5.747671819196675, -0.9408309030864932, 2.967655969930743] + ])), 'fro') < 1e-3) + }) +}) diff --git a/types/index.d.ts b/types/index.d.ts index e14873378..06a54f96b 100644 --- a/types/index.d.ts +++ b/types/index.d.ts @@ -465,6 +465,11 @@ declare namespace math { R: MathCollection } + interface SchurDecomposition { + U: MathCollection + T: MathCollection + } + interface FractionDefinition { a: number b: number @@ -1807,6 +1812,41 @@ declare namespace math { */ expm(x: Matrix): Matrix + /** + * Solves the real-valued Sylvester equation AX-XB=C for X, where A, B and C are + * matrices of appropriate dimensions, being A and B squared. The method used is + * the Bartels-Stewart algorithm. + * https://en.wikipedia.org/wiki/Sylvester_equation + * @param A Matrix A + * @param B Matrix B + * @param C Matrix C + * @returns Matrix X, solving the Sylvester equation + */ + sylvester( + A: Matrix | MathArray, + B: Matrix | MathArray, + C: Matrix | MathArray + ): Matrix | MathArray + + /** + * Performs a real Schur decomposition of the real matrix A = UTU' where U is orthogonal + * and T is upper quasi-triangular. + * https://en.wikipedia.org/wiki/Schur_decomposition + * @param A Matrix A + * @returns Object containing both matrix U and T of the Schur Decomposition A=UTU' + */ + schur(A: Matrix | MathArray): SchurDecomposition + + /** + * Solves the Continuous-time Lyapunov equation AP+PA'=Q for P, where Q is a positive semidefinite + * matrix. + * https://en.wikipedia.org/wiki/Lyapunov_equation + * @param A Matrix A + * @param Q Matrix Q + * @returns Matrix P solution to the Continuous-time Lyapunov equation AP+PA'=Q + */ + lyap(A: Matrix | MathArray, Q: Matrix | MathArray): Matrix | MathArray + /** * Create a 2-dimensional identity matrix with size m x n or n x n. The * matrix has ones on the diagonal and zeros elsewhere. @@ -3572,6 +3612,9 @@ declare namespace math { invDependencies: FactoryFunctionMap expmDependencies: FactoryFunctionMap sqrtmDependencies: FactoryFunctionMap + sylvesterDependencies: FactoryFunctionMap + schurDependencies: FactoryFunctionMap + lyapDependencies: FactoryFunctionMap divideDependencies: FactoryFunctionMap distanceDependencies: FactoryFunctionMap intersectDependencies: FactoryFunctionMap @@ -5272,6 +5315,26 @@ declare namespace math { expm(this: MathJsChain): MathJsChain + /** + * Performs a real Schur decomposition of the real matrix A = UTU' where U is orthogonal + * and T is upper quasi-triangular. + * https://en.wikipedia.org/wiki/Schur_decomposition + * @returns Object containing both matrix U and T of the Schur Decomposition A=UTU' + */ + schur(this: MathJsChain): SchurDecomposition + + /** + * Solves the Continuous-time Lyapunov equation AP+PA'=Q for P, where Q is a positive semidefinite + * matrix. + * https://en.wikipedia.org/wiki/Lyapunov_equation + * @param Q Matrix Q + * @returns Matrix P solution to the Continuous-time Lyapunov equation AP+PA'=Q + */ + lyap( + this: MathJsChain, + Q: Matrix | MathArray + ): MathJsChain + /** * Create a 2-dimensional identity matrix with size m x n or n x n. The * matrix has ones on the diagonal and zeros elsewhere.