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