mirror of
https://github.com/josdejong/mathjs.git
synced 2026-01-18 14:59:29 +00:00
Added functions: schur, sylvester, and lyap (#2646)
* Added matrix functions `schur`, `sylvester` and `lyap` * Added docs for `schur`, `sylvester` and `lyap` * Added unit tests for `schur`, `sylvester` and `lyap` * Fixed lint and tests errors * fixed typescript and added Matrix + Array tests * lint fixed * fixed example in sylvester.js doc * Fixed docs and ci * fixed definition of lyap and sylvester * remark on diff defs for lyap * rm a trailing space Co-authored-by: Jos de Jong <wjosdejong@gmail.com>
This commit is contained in:
parent
b9678ca029
commit
a2e85ca698
1
AUTHORS
1
AUTHORS
@ -199,6 +199,7 @@ FelixSelter <55546882+FelixSelter@users.noreply.github.com>
|
||||
Hansuku <1556207795@qq.com>
|
||||
Windrill <windrill@outlook.com>
|
||||
Carl Osterwisch <costerwi@gmail.com>
|
||||
Lucas Egidio <lucasegidio1@gmail.com>
|
||||
Isaac B <66892358+isaacbyr@users.noreply.github.com>
|
||||
Ajinkya Chikhale <86607732+ajinkyac03@users.noreply.github.com>
|
||||
Abdullah Adeel <64294045+mabdullahadeel@users.noreply.github.com>
|
||||
|
||||
@ -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,
|
||||
|
||||
15
src/expression/embeddedDocs/function/matrix/lyap.js
Normal file
15
src/expression/embeddedDocs/function/matrix/lyap.js
Normal file
@ -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'
|
||||
]
|
||||
}
|
||||
15
src/expression/embeddedDocs/function/matrix/schur.js
Normal file
15
src/expression/embeddedDocs/function/matrix/schur.js
Normal file
@ -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'
|
||||
]
|
||||
}
|
||||
15
src/expression/embeddedDocs/function/matrix/sylvester.js
Normal file
15
src/expression/embeddedDocs/function/matrix/sylvester.js
Normal file
@ -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'
|
||||
]
|
||||
}
|
||||
@ -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'
|
||||
|
||||
77
src/function/algebra/decomposition/schur.js
Normal file
77
src/function/algebra/decomposition/schur.js
Normal file
@ -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 }
|
||||
}
|
||||
})
|
||||
61
src/function/algebra/lyap.js
Normal file
61
src/function/algebra/lyap.js
Normal file
@ -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()
|
||||
}
|
||||
})
|
||||
})
|
||||
145
src/function/algebra/sylvester.js
Normal file
145
src/function/algebra/sylvester.js
Normal file
@ -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
|
||||
}
|
||||
})
|
||||
@ -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) {
|
||||
|
||||
63
test/unit-tests/function/algebra/decomposition/schur.test.js
Normal file
63
test/unit-tests/function/algebra/decomposition/schur.test.js
Normal file
@ -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)
|
||||
})
|
||||
})
|
||||
49
test/unit-tests/function/algebra/lyap.test.js
Normal file
49
test/unit-tests/function/algebra/lyap.test.js
Normal file
@ -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)
|
||||
})
|
||||
})
|
||||
55
test/unit-tests/function/algebra/sylvester.test.js
Normal file
55
test/unit-tests/function/algebra/sylvester.test.js
Normal file
@ -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)
|
||||
})
|
||||
})
|
||||
63
types/index.d.ts
vendored
63
types/index.d.ts
vendored
@ -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<Matrix>): MathJsChain<Matrix>
|
||||
|
||||
/**
|
||||
* 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<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 Q Matrix Q
|
||||
* @returns Matrix P solution to the Continuous-time Lyapunov equation AP+PA'=Q
|
||||
*/
|
||||
lyap(
|
||||
this: MathJsChain<Matrix | MathArray>,
|
||||
Q: Matrix | MathArray
|
||||
): MathJsChain<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.
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user