Mercurial > repos > public > sbplib

comparison +time/+expint/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
children
comparison
equal deleted inserted replaced
705:e6fbdc9ccfc4 706:a95c89436916
1 function y = Arnoldi(A,v,t,tol,m)
2
3 n = size (v,1);
4 V = zeros(n ,m+1);
5 H = zeros(m+1,m);
6
7 beta = norm(v);
8 V(:,1) = v/beta;
9 resnorm = inf;
10
11 j=0;
12
13 while resnorm > tol
14 j = j+1;
15 z = A*V(:,j);
16 for i=1:j
17 H(i,j) = z'*V(:,i);
18 z = z - H(i,j)*V(:,i);
19 end
20 H(j+1,j) = norm(z);
21 e1 = zeros(j,1); e1(1) = 1;
22 ej = zeros(j,1); ej(j) = 1;
23 s = [0.01, 1/3, 2/3, 1]*t;
24 for q=1:length(s)
25 u = expm(-s(q)*H(1:j,1:j))*e1;
26 beta_j(q) = -H(j+1,j)* (ej'*u);
27 end
28 resnorm = norm(beta_j,'inf');
29 fprintf('j = %d, resnorm = %.2e\n',j,resnorm);
30 if resnorm<=toler
31 break
32 elseif j==m
33 disp('warning: no convergence within m steps');
34 end
35 V(:,j+1) = z/H(j+1,j);
36 end
37 y = V(:,1:j)*(beta*u);