Mercurial > repos > public > sbplib

comparison +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
comparison
equal deleted inserted replaced
512:4ef2d2a493f1 513:bc39bb984d88
1 function y = expm_Arnoldi(A,v,t,toler,m)
2 %
3 % y = expm_Arnoldi(A,v,t,toler,m)
4 %
5 % computes $y = \exp(-t A) v$
6 % input: A (n x n)-matrix, v n-vector, t>0 time interval,
7 % toler>0 tolerance, m maximal Krylov dimension
8 %
9 % Copyright (c) 2012 by M.A. Botchev
10 % Permission to copy all or part of this work is granted,
11 % provided that the copies are not made or distributed
12 % for resale, and that the copyright notice and this
13 % notice are retained.
14 %
15 % THIS WORK IS PROVIDED ON AN "AS IS" BASIS. THE AUTHOR
16 % PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED,
17 % REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS
18 % MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE.
19 %
20 n = size (v,1);
21 V = zeros(n ,m+1);
22 H = zeros(m+1,m);
23
24 beta = norm(v);
25 V(:,1) = v/beta;
26 resnorm = inf;
27
28 j=0;
29
30 while resnorm > toler
31 j = j+1;
32 w = A*V(:,j);
33 for i=1:j
34 H(i,j) = w'*V(:,i);
35 w = w - H(i,j)*V(:,i);
36 end
37 H(j+1,j) = norm(w);
38 e1 = zeros(j,1); e1(1) = 1;
39 ej = zeros(j,1); ej(j) = 1;
40 s = [0.01, 1/3, 2/3, 1]*t;
41 for q=1:length(s)
42 u = expm(-s(q)*H(1:j,1:j))*e1;
43 beta_j(q) = -H(j+1,j)* (ej'*u);
44 end
45 resnorm = norm(beta_j,'inf');
46 % fprintf('j = %d, resnorm = %.2e\n',j,resnorm);
47 if resnorm<=toler
48 break
49 elseif j==m
50 disp('warning: no convergence within m steps');
51 end
52 V(:,j+1) = w/H(j+1,j);
53 end
54 y = V(:,1:j)*(beta*u);