Mercurial > matrix-functions

comparison mftoolbox/expm_frechet_quad.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function R = expm_frechet_quad(A,E,theta,rule,k)
2 %EXPM_FRECHET_QUAD Frechet derivative of matrix exponential via quadrature.
3 % L = EXPM_FRECHET_QUAD(A,E,THETA,RULE) is an approximation to the
4 % Frechet derivative of the matrix exponential at A in the direction E
5 % intended to have norm of the correct order of magnitude.
6 % It is obtained from the repeated trapezium rule (RULE = 'T'),
7 % the repeated Simpson rule (RULE = 'S', default),
8 % or the repeated midpoint rule (RULE = 'M').
9 % L = EXPM_FRECHET_QUAD(A,E,THETA,RULE,k) uses either matrix
10 % exponentials (k = 0, the default) or repeated squaring (k = 1)
11 % in the final phase of the algorithm.
12
13 % A is scaled so that norm(A/2^s) <= THETA. Defalt: THETA = 1/2.
14
15 if nargin < 3 || isempty(theta), theta = 1/2; end
16 if nargin < 4 || isempty(rule), rule = 'S'; end
17 if nargin < 5, k = 0; end
18
19 s = ceil( log2(norm(A,1)/theta) );
20 As = A/2^s;
21
22 X = expm(As);
23
24 switch upper(rule)
25
26 case 'T'
27 R = 2^(-s) * (X*E + E*X)/2 ;
28
29 case 'S'
30 Xmid = expm(As/2);
31 R = 2^(-s) * (X*E + 4*Xmid*E*Xmid + E*X)/6;
32
33 case 'M'
34 Xmid = expm(As/2);
35 R = 2^(-s) * Xmid*E*Xmid;
36
37 otherwise
38 error('Illegal value of RULE.')
39
40 end
41
42 for i = s:-1:1
43 if i < s
44 if k == 0
45 X = expm(2^(-i)*A);
46 else
47 X = X^2;
48 end
49 end
50 R = X*R + R*X;
51 end