Mercurial > repos > public > sbplib

annotate +time/+expint/expm_Arnoldi.m @ 513:bc39bb984d88 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added arnoldi krylov subspace approximation
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 26 Jun 2017 20:15:54 +0200
parents
children 8434063ed162
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
513
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
1 function y = expm_Arnoldi(A,v,t,toler,m)
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
2 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
3 % y = expm_Arnoldi(A,v,t,toler,m)
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
4 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
5 % computes $y = \exp(-t A) v$
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
6 % input: A (n x n)-matrix, v n-vector, t>0 time interval,
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
7 % toler>0 tolerance, m maximal Krylov dimension
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
8 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
9 % Copyright (c) 2012 by M.A. Botchev
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
10 % Permission to copy all or part of this work is granted,
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
11 % provided that the copies are not made or distributed
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
12 % for resale, and that the copyright notice and this
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
13 % notice are retained.
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
14 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
15 % THIS WORK IS PROVIDED ON AN "AS IS" BASIS. THE AUTHOR
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
16 % PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED,
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
17 % REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
18 % MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE.
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
19 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
20 n = size (v,1);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
21 V = zeros(n ,m+1);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
22 H = zeros(m+1,m);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
23
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
24 beta = norm(v);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
25 V(:,1) = v/beta;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 resnorm = inf;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
27
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
28 j=0;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
29
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
30 while resnorm > toler
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
31 j = j+1;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
32 w = A*V(:,j);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
33 for i=1:j
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
34 H(i,j) = w'*V(:,i);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
35 w = w - H(i,j)*V(:,i);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
36 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
37 H(j+1,j) = norm(w);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
38 e1 = zeros(j,1); e1(1) = 1;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
39 ej = zeros(j,1); ej(j) = 1;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
40 s = [0.01, 1/3, 2/3, 1]*t;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
41 for q=1:length(s)
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
42 u = expm(-s(q)*H(1:j,1:j))*e1;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
43 beta_j(q) = -H(j+1,j)* (ej'*u);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
44 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
45 resnorm = norm(beta_j,'inf');
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
46 % fprintf('j = %d, resnorm = %.2e\n',j,resnorm);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
47 if resnorm<=toler
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
48 break
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
49 elseif j==m
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
50 disp('warning: no convergence within m steps');
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
51 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
52 V(:,j+1) = w/H(j+1,j);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
53 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
54 y = V(:,1:j)*(beta*u);