|
706
|
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);
|
