Arnoldi method to find eigenvalues

I'm trying to implement the Arnoldi method with the inverse power method to find eigenvalues of large pencil matrix.

I have followed the practical implementation in Saad's book and I started with a small matrix to check if the code work well. As I understand H is a square matrix and has size of the number of the iterations but the resulted H is of size 3x2 and V is 4x3. I'm not sure if this is correct and I do'nt know how I can find the eigenvalues of H and the corresponding eigenvectors. Here is my attempt, and I really appreciate any help..

A=[1 0 0 0;0 17 0 0;0 0 -10 0 ;0 0 0 5]
n=length(A)
B=eye(n)
a=0  %interval [a,b] 
b=10
m=2;  %iterations
x=ones(n,1); %constant vector
shift=0.5;
V = zeros(n,m);  %orthonormal basis of Krylov space
V(:,1) = x/norm(x);
H = zeros(n,m);  %upper Hessenberg matrix
C=A-shift*B; 
[L, U] = lu(C, 'vector');
for j=1:m     
    w=U\(L\(B*V(:,j)));
    for i=1:j    %Gram-schmidt 
        H(i,j) = V(:,i)'*w;
        w = w - V(:,i)*H(i,j);
    end
        H(j+1,j) = norm(w,2);
       V(:,j+1) = w/H(j+1,j); 
end 
   

Answer 1

[Pavl M.]

Hi, happy to help:

Use next code:

 clc
 clear all
 close all

%non-square matrix:
%A=[1 0 0 0;0 17 0 0;0 0 -10 0 ;0 0 0 5;0 0 0 0]

A=[1 0 0 0;0 17 0 0;0 0 -10 0 ;0 0 0 5]

A = 4×4

1 0 0 0 0 17 0 0 0 0 -10 0 0 0 0 5

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%DMD code:
%version 1:

n1=size(A,1);
n2=size(A,2);

B=eye(n1,n2);

a=0  %interval [a,b]

a = 0

b=10

b = 10

m=n2;  %iterations
x=ones(n2,1); %constant vector

shift=0.5;

V = zeros(n2,m);  %orthonormal basis of Krylov space
V(:,1) = x/norm(x);

H = zeros(min(m+1,m),n1);  %upper Hessenberg matrix

W = zeros(n1,m);
Z = zeros(n1,m);

C=A-shift*B;

[L, U] = lu(C, 'vector');

for j=1:m

    w=U\(L\(B*V(:,j)))
    z = A*V(:,j)
    W(:,j) = w;
    Z(:,j) = z;
    
    for i=1:j    %Gram-schmidt
        H(i,j) = V(:,i)'*w;
        w = w - V(:,i)*H(i,j);
    end

    if j < n1
        H(j+1,j) = norm(w,2);
        if H(j+1,j) == 0, return, end
        V(:,j+1) = w/H(j+1,j);
    end

end

w = 4×1

1.0000 0.0303 -0.0476 0.1111

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

z = 4×1

0.5000 8.5000 -5.0000 2.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

w = 4×1

1.7168 -0.0174 0.0361 -0.0426

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

z = 4×1

0.8584 -4.8835 3.7932 -0.9590

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

w = 4×1

0.2294 -0.0042 -0.0645 -0.1607

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

z = 4×1

0.1147 -1.1711 -6.7746 -3.6164

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

w = 4×1

0.0110 0.0493 0.0365 -0.0969

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

z = 4×1

0.0055 13.8395 3.8362 -2.1799

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

H

H = 4×4

0.5469 0.8464 -0.0000 0.0000 0.8464 1.4731 0.2534 0.0000 0 0.2534 0.0991 0.0927 0 0 0.0927 0.0684

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

V

V = 4×4

0.5000 0.8584 0.1147 0.0055 0.5000 -0.2873 -0.0689 0.8141 0.5000 -0.3793 0.6775 -0.3836 0.5000 -0.1918 -0.7233 -0.4360

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

W

W = 4×4

1.0000 1.7168 0.2294 0.0110 0.0303 -0.0174 -0.0042 0.0493 -0.0476 0.0361 -0.0645 0.0365 0.1111 -0.0426 -0.1607 -0.0969

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Z

Z = 4×4

0.5000 0.8584 0.1147 0.0055 8.5000 -4.8835 -1.1711 13.8395 -5.0000 3.7932 -6.7746 3.8362 2.5000 -0.9590 -3.6164 -2.1799

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%find dynamical modes:

%vector of mode amplitudes and ? estimation:
v1 = H*A

v1 = 4×4

0.5469 14.3892 0.0000 0.0000 0.8464 25.0426 -2.5343 0.0000 0 4.3084 -0.9915 0.4634 0 0 -0.9269 0.3422

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

v2 = A*diag(Z)

v2 = 4×1

0.5000 -83.0188 67.7460 -10.8993

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

v23 = diag(z)*A

v23 = 4×4

0.0055 0 0 0 0 235.2707 0 0 0 0 -38.3617 0 0 0 0 -10.8993

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

v11 = W*A

v11 = 4×4

1.0000 29.1848 -2.2943 0.0550 0.0303 -0.2960 0.0417 0.2467 -0.0476 0.6141 0.6452 0.1827 0.1111 -0.7245 1.6073 -0.4844

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

v22 = Z*A

v22 = 4×4

0.5000 14.5924 -1.1471 0.0275 8.5000 -83.0188 11.7107 69.1973 -5.0000 64.4848 67.7460 19.1809 2.5000 -16.3023 36.1643 -10.8993

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

display('Accuracy measure 1:')

Accuracy measure 1:

acc = mean(mean(abs(V-Z)))

acc = 3.4670