Mercurial > repos > public > sbplib

annotate +time/+expint/expm_Arnoldi.m @ 706:a95c89436916 feature/optim

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
add arnoldi
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 06 Nov 2017 11:42:15 +0100
parents 8434063ed162
children
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)
697
8434063ed162 1 objective optimization finished
Ylva Rydin <ylva.rydin@telia.com>
parents: 513
diff changeset
2
8434063ed162 1 objective optimization finished
Ylva Rydin <ylva.rydin@telia.com>
parents: 513
diff changeset
3 global iter
8434063ed162 1 objective optimization finished
Ylva Rydin <ylva.rydin@telia.com>
parents: 513
diff changeset
4
513
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
5 % y = expm_Arnoldi(A,v,t,toler,m)
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
6 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
7 % computes $y = \exp(-t A) v$
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
8 % 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
9 % toler>0 tolerance, m maximal Krylov dimension
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
10 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
11 % Copyright (c) 2012 by M.A. Botchev
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
12 % 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
13 % provided that the copies are not made or distributed
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
14 % for resale, and that the copyright notice and this
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
15 % notice are retained.
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
16 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
17 % 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
18 % PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED,
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
19 % REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
20 % MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE.
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
21 %
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
22 n = size (v,1);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
23 V = zeros(n ,m+1);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
24 H = zeros(m+1,m);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
25
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
26 beta = norm(v);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
27 V(:,1) = v/beta;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
28 resnorm = inf;
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 j=0;
697
8434063ed162 1 objective optimization finished
Ylva Rydin <ylva.rydin@telia.com>
parents: 513
diff changeset
31 iter = 0;
8434063ed162 1 objective optimization finished
Ylva Rydin <ylva.rydin@telia.com>
parents: 513
diff changeset
32
513
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
33
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
34 while resnorm > toler
697
8434063ed162 1 objective optimization finished
Ylva Rydin <ylva.rydin@telia.com>
parents: 513
diff changeset
35 iter = iter +1;
513
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
36 j = j+1;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
37 w = A*V(:,j);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
38 for i=1:j
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
39 H(i,j) = w'*V(:,i);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
40 w = w - H(i,j)*V(:,i);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
41 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
42 H(j+1,j) = norm(w);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
43 e1 = zeros(j,1); e1(1) = 1;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
44 ej = zeros(j,1); ej(j) = 1;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
45 s = [0.01, 1/3, 2/3, 1]*t;
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
46 for q=1:length(s)
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
47 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
48 beta_j(q) = -H(j+1,j)* (ej'*u);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
49 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
50 resnorm = norm(beta_j,'inf');
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
51 % fprintf('j = %d, resnorm = %.2e\n',j,resnorm);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
52 if resnorm<=toler
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
53 break
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
54 elseif j==m
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
55 disp('warning: no convergence within m steps');
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
56 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
57 V(:,j+1) = w/H(j+1,j);
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
58 end
bc39bb984d88 Added arnoldi krylov subspace approximation
Ylva Rydin <ylva.rydin@telia.com>
parents:
diff changeset
59 y = V(:,1:j)*(beta*u);