1 | 'use strict'
|
2 |
|
3 | const format = require('../../utils/string').format
|
4 |
|
5 | function factory (type, config, load, typed) {
|
6 | const abs = load(require('../arithmetic/abs'))
|
7 | const add = load(require('../arithmetic/add'))
|
8 | const identity = load(require('./identity'))
|
9 | const inv = load(require('./inv'))
|
10 | const multiply = load(require('../arithmetic/multiply'))
|
11 |
|
12 | const SparseMatrix = type.SparseMatrix
|
13 |
|
14 | |
15 |
|
16 |
|
17 |
|
18 |
|
19 |
|
20 |
|
21 |
|
22 |
|
23 |
|
24 |
|
25 |
|
26 |
|
27 |
|
28 |
|
29 |
|
30 |
|
31 |
|
32 |
|
33 |
|
34 |
|
35 |
|
36 |
|
37 |
|
38 |
|
39 | const expm = typed('expm', {
|
40 |
|
41 | 'Matrix': function (A) {
|
42 |
|
43 | const size = A.size()
|
44 |
|
45 | if (size.length !== 2 || size[0] !== size[1]) {
|
46 | throw new RangeError('Matrix must be square ' +
|
47 | '(size: ' + format(size) + ')')
|
48 | }
|
49 |
|
50 | const n = size[0]
|
51 |
|
52 |
|
53 |
|
54 | const eps = 1e-15
|
55 |
|
56 |
|
57 |
|
58 |
|
59 |
|
60 |
|
61 |
|
62 | const infNorm = infinityNorm(A)
|
63 |
|
64 |
|
65 |
|
66 | const params = findParams(infNorm, eps)
|
67 | const q = params.q
|
68 | const j = params.j
|
69 |
|
70 |
|
71 |
|
72 |
|
73 |
|
74 |
|
75 |
|
76 |
|
77 | const Apos = multiply(A, Math.pow(2, -j))
|
78 |
|
79 |
|
80 | let N = identity(n)
|
81 | let D = identity(n)
|
82 |
|
83 |
|
84 | let factor = 1
|
85 |
|
86 |
|
87 | let AposToI = Apos
|
88 | let alternate = -1
|
89 |
|
90 | for (let i = 1; i <= q; i++) {
|
91 | if (i > 1) {
|
92 | AposToI = multiply(AposToI, Apos)
|
93 | alternate = -alternate
|
94 | }
|
95 | factor = factor * (q - i + 1) / ((2 * q - i + 1) * i)
|
96 |
|
97 | N = add(N, multiply(factor, AposToI))
|
98 | D = add(D, multiply(factor * alternate, AposToI))
|
99 | }
|
100 |
|
101 | let R = multiply(inv(D), N)
|
102 |
|
103 |
|
104 | for (let i = 0; i < j; i++) {
|
105 | R = multiply(R, R)
|
106 | }
|
107 |
|
108 | return type.isSparseMatrix(A)
|
109 | ? new SparseMatrix(R)
|
110 | : R
|
111 | }
|
112 |
|
113 | })
|
114 |
|
115 | function infinityNorm (A) {
|
116 | const n = A.size()[0]
|
117 | let infNorm = 0
|
118 | for (let i = 0; i < n; i++) {
|
119 | let rowSum = 0
|
120 | for (let j = 0; j < n; j++) {
|
121 | rowSum += abs(A.get([i, j]))
|
122 | }
|
123 | infNorm = Math.max(rowSum, infNorm)
|
124 | }
|
125 | return infNorm
|
126 | }
|
127 |
|
128 | |
129 |
|
130 |
|
131 |
|
132 |
|
133 | function findParams (infNorm, eps) {
|
134 | const maxSearchSize = 30
|
135 | for (let k = 0; k < maxSearchSize; k++) {
|
136 | for (let q = 0; q <= k; q++) {
|
137 | const j = k - q
|
138 | if (errorEstimate(infNorm, q, j) < eps) {
|
139 | return { q: q, j: j }
|
140 | }
|
141 | }
|
142 | }
|
143 | throw new Error('Could not find acceptable parameters to compute the matrix exponential (try increasing maxSearchSize in expm.js)')
|
144 | }
|
145 |
|
146 | |
147 |
|
148 |
|
149 |
|
150 | function errorEstimate (infNorm, q, j) {
|
151 | let qfac = 1
|
152 | for (let i = 2; i <= q; i++) {
|
153 | qfac *= i
|
154 | }
|
155 | let twoqfac = qfac
|
156 | for (let i = q + 1; i <= 2 * q; i++) {
|
157 | twoqfac *= i
|
158 | }
|
159 | const twoqp1fac = twoqfac * (2 * q + 1)
|
160 |
|
161 | return 8.0 *
|
162 | Math.pow(infNorm / Math.pow(2, j), 2 * q) *
|
163 | qfac * qfac / (twoqfac * twoqp1fac)
|
164 | }
|
165 |
|
166 | expm.toTex = { 1: `\\exp\\left(\${args[0]}\\right)` }
|
167 |
|
168 | return expm
|
169 | }
|
170 |
|
171 | exports.name = 'expm'
|
172 | exports.factory = factory
|