Mercurial > matrix-functions

annotate 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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
1 function R = expm_frechet_quad(A,E,theta,rule,k)
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
2 %EXPM_FRECHET_QUAD Frechet derivative of matrix exponential via quadrature.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
3 % L = EXPM_FRECHET_QUAD(A,E,THETA,RULE) is an approximation to the
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
4 % Frechet derivative of the matrix exponential at A in the direction E
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
5 % intended to have norm of the correct order of magnitude.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
6 % It is obtained from the repeated trapezium rule (RULE = 'T'),
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
7 % the repeated Simpson rule (RULE = 'S', default),
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
8 % or the repeated midpoint rule (RULE = 'M').
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
9 % L = EXPM_FRECHET_QUAD(A,E,THETA,RULE,k) uses either matrix
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
10 % exponentials (k = 0, the default) or repeated squaring (k = 1)
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
11 % in the final phase of the algorithm.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
12
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
13 % A is scaled so that norm(A/2^s) <= THETA. Defalt: THETA = 1/2.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
14
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
15 if nargin < 3 || isempty(theta), theta = 1/2; end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
16 if nargin < 4 || isempty(rule), rule = 'S'; end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
17 if nargin < 5, k = 0; end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
18
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
19 s = ceil( log2(norm(A,1)/theta) );
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
20 As = A/2^s;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
21
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
22 X = expm(As);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
23
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
24 switch upper(rule)
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
25
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
26 case 'T'
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
27 R = 2^(-s) * (X*E + E*X)/2 ;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
28
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
29 case 'S'
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
30 Xmid = expm(As/2);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
31 R = 2^(-s) * (X*E + 4*Xmid*E*Xmid + E*X)/6;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
32
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
33 case 'M'
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
34 Xmid = expm(As/2);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
35 R = 2^(-s) * Xmid*E*Xmid;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
36
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
37 otherwise
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
38 error('Illegal value of RULE.')
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
39
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
40 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
41
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
42 for i = s:-1:1
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
43 if i < s
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
44 if k == 0
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
45 X = expm(2^(-i)*A);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
46 else
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
47 X = X^2;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
48 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
49 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
50 R = X*R + R*X;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
51 end