Hi Lennart,
It seems there are a few things that are not correct with the last part of the code:
1. Instead of li(qq), it should be li(qq,qq).
3. The code below gives you the desired results:
for qq = 1:length(L1)
right1(:, qq) = (A - L1(qq,qq) * B) * X1(:, qq);
left1(qq, :) = Y1(:, qq)' * (A - L1(qq,qq) * B);
right2(:, qq) = (A - L21(qq,qq) * B) * X2(:, qq);
left2(qq, :) = Y2(:, qq)' * (A - L22(qq,qq)' * B);
right3(:, qq) = (A - L3(qq,qq) * B) * X3(:, qq);
left3(qq, :) = Y3(:, qq)' * (A - L3(qq,qq) * B);
right4(:, qq) = (A - L41(qq,qq) * B) * X4(:, qq);
left4(qq, :) = Y4(:, qq)' * (A - L42(qq,qq)' * B);
right5(:, qq) = (A - L5(qq,qq) * B) * X5(:, qq);
left5(qq, :) = Y5(:, qq)' * (A - L5(qq,qq) * B);
right6(:, qq) = (A - L61(qq,qq) * B) * X6(:, qq);
left6(qq, :) = Y6(:, qq)' * (A - L62(qq,qq)' * B);
end
[left1 right1;
left2 right2;
left3 right3;
left4 right4;
left5 right5;
left6 right6]
Please note that I have changed Li2(qq,qq) to Li2(qq,qq)' for the computation of left2, left4 and left6 in the 2nd, 4th and 6th case. Here is why: As per the documentation of eig, [Y2, L22] = eigs(A', B');
ensures that the solution satisfies the following:
A'*Y2 = B'*Y2*L22
Y2'*A = L22'*Y2'*B
Y2'*A - L22'*Y2'*B = 0
Y2'*(A-L22'*B) = 0
(since values of L22 are taken one by one in your calculations, so it can behave as a scalar and moving its position does not change the result). However, taking its (complex conjugate) transpose is crucial since it is a complex number. Instead of computing the residuals inside a loop, it can be written as one single statement (treating L22 as a matrix):
left2 = Y2'*A - L22'*Y2'*B;
The same reasoning applies for cases 4 and 6. Please note this reasoning does not apply to cases 1, 3 and 5 and the original form of equation stays valid.
I hope the above helps.
Nalini
Best Answer